#include <math.h>
#include <stdint.h>
#include <stdio.h>

/* fixed point math definitions */
typedef int32_t fpnum; /* fixed-point number */

/* fixed point number characteristics */
enum {
	FPBITS  = 32, /* number of bits in a fixed point number */
	FPSCALE =  8, /* number of binary digits before the decimal point */
};

/*
 * macro to build a fixed point number from a part before the decimal
 * point and another one after it.  b must be a positive integer.
 */
#define MKFPNUM(a, b) ((fpnum)(a) << FPBITS - FPSCALE | (fpnum)(b))

/*
 * lookup table for fplog2(x).  Each entry i contains the 32 binary places of
 * log2(1+2^(-i-1)) after the decimal point as a 32 bit number.
 */
uint32_t logtab[32] = {
	2512394810, 1382670639,  729822324,  375650043,
	 190671291,   96069025,   48220695,   24157255,
	  12090400,    6048149,    3024812,    1512591,
	    756342,     378182,     189094,      94548,
	     47274,      23637,      11819,       5909,
	      2955,       1477,        739,        369,
	       185,         92,         46,         23,
	        12,          6,          3,          1,
};

/* compute log2(x) */
extern fpnum
fplog2(fpnum x)
{
	fpnum k, y;
	uint32_t k1, x1, y1, a;
	int i;

	/* the logarithm of x < 0 is a complex number, return -1 */
	if (x <= 0)
		return (-1L);

	/*
	 * step 1: scale x so it's in the range 1 <= x < 2 by multiplying
	 * with 2^-k for an appropriate k.  As we only care about the places
	 * after the decimal point here, we shift x into a fixed point form
	 * with 0 places before the decimal point, i.e. the scale is 2^-FPBITS.
	 * The value k takes the range FPBITS - FPSCALE <= k < FPSCALE:
	 */
	k1 = __builtin_clz(x); /* compute k - FPSCALE - 1 */
	x1 = x << k1 + 1; /* x1 is x scaled to range 1 <= x < 2, but just the places after the decimal point */
	k1 = FPSCALE - k1 - 1; /* adjust k to actually be that 2^-k factor */
	k = MKFPNUM(k1, 0); /* and turn it into a fixed point number */

	/*
	 * step 2: approximate x1 as a product of factors of the form (1+2^a).
	 * In y1 (represented like x1 as just places after the decimal point),
	 * we accumulate the logarithm of this approximation.  Note that as
	 * x1 = (1+2^a0)*(1+2^a1)*... we have y1 = log2(1+2^a0)+log2(1+2^a1)+...
	 * each of the summands of y1 is looked up as an entry in logtab.
	 * a is the accumulator used to approximate x1.  It starts out as a=1.0,
	 * but as we only care about places after the decimal point, the 1
	 * before the decimal point is not stored.  Just remember that it's
	 * there.  This approximation is done for up to 32 binary places, but I
	 * am sure just FPBITS - FPSCALE rounds would suffice.
	 */
	a = 0;
	y1 = 0;
	for (i = 1; i <= 32; i++) {
		/* 1 << FPBITS - i accounts for the implicit 1 */
		if (a + (a >> i) + (1 << FPBITS - i) < x1) {
			a += (a >> i) + (1 << FPBITS - i);
			y1 += logtab[i - 1];
		}
	}

	/*
	 * step 3: form the result from y1 and k.
	 */
	y = MKFPNUM(0, y1 >> FPSCALE); /* convert y to a fixed point number */
	y += k;                        /* apply scale */

	return (y);
}

/*
 * Test program: read numbers from console, and then print the logarithm both
 * using our routine as well as the floating point log2 function.
 */
extern int
main(void)
{
	/* infinite loop */
	for (;;) {
		double x, y, yfd;
		fpnum xf, yf;

		if (scanf("%lf", &x) != 1)
			break;

		xf = (fpnum)scalbln(x, FPBITS - FPSCALE);
		y = log2(x);
		yf = fplog2(xf);
		yfd = scalbln((double)yf, FPSCALE - FPBITS);

		printf("log2 %g  fplog2 %g  diff %g\n", y, yfd, y - yfd);
	}

	return (0);
}