vclog.c

/***********************************************************************
*                                                                      *
*               This software is part of the ast package               *
*          Copyright (c) 2003-2011 AT&T Intellectual Property          *
*                      and is licensed under the                       *
*                 Eclipse Public License, Version 1.0                  *
*                    by AT&T Intellectual Property                     *
*                                                                      *
*                A copy of the License is available at                 *
*          http://www.eclipse.org/org/documents/epl-v10.html           *
*         (with md5 checksum b35adb5213ca9657e911e9befb180842)         *
*                                                                      *
*              Information and Software Systems Research               *
*                            AT&T Research                             *
*                           Florham Park NJ                            *
*                                                                      *
*                   Phong Vo <[email protected]>                    *
*                                                                      *
***********************************************************************/
#include	"vchdr.h"

/*	Function to approximate log_base_2(unsigned int v)".
**	Return 0. if v is zero.
**
**	The basic idea is to write v = (2^k)*w + r = (2^k)*w*(1 + r/((2^k)*w))
**	so that log2(v) = k + log2(w) + log2(1 + r/((2^k)*w)).
**	Now find k so that w < 256. Then, k + log2(w) is a pretty good
**	approximation of log2(v) up to about 2 decimal places. Both k and log2(w)
**	can be computed using table look-up with a couple of small tables.
**
**	For further precision, we add the first term of the Taylor expansion of
**	log2(1 + r/((2^k)*w)). This term can be computed quickly with a few
**	bit masks and divisions.
**
**	Compile with -DPRECISE=0 if precision can be fudged.
**
**	Written by Kiem-Phong Vo ([email protected])
*/

#ifndef PRECISE
#define PRECISE	1
#endif

static unsigned char _Vcpow2[256] = {
 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7 
};

static double _Vclog2[256] = {
 0.0000000, 0.0000000, 1.0000000, 1.5849625, 2.0000000, 2.3219281, 2.5849625, 2.8073549,
 3.0000000, 3.1699250, 3.3219281, 3.4594316, 3.5849625, 3.7004397, 3.8073549, 3.9068906,
 4.0000000, 4.0874628, 4.1699250, 4.2479275, 4.3219281, 4.3923174, 4.4594316, 4.5235620,
 4.5849625, 4.6438562, 4.7004397, 4.7548875, 4.8073549, 4.8579810, 4.9068906, 4.9541963,
 5.0000000, 5.0443941, 5.0874628, 5.1292830, 5.1699250, 5.2094534, 5.2479275, 5.2854022,
 5.3219281, 5.3575520, 5.3923174, 5.4262648, 5.4594316, 5.4918531, 5.5235620, 5.5545889,
 5.5849625, 5.6147098, 5.6438562, 5.6724253, 5.7004397, 5.7279205, 5.7548875, 5.7813597,
 5.8073549, 5.8328900, 5.8579810, 5.8826430, 5.9068906, 5.9307373, 5.9541963, 5.9772799,
 6.0000000, 6.0223678, 6.0443941, 6.0660892, 6.0874628, 6.1085245, 6.1292830, 6.1497471,
 6.1699250, 6.1898246, 6.2094534, 6.2288187, 6.2479275, 6.2667865, 6.2854022, 6.3037807,
 6.3219281, 6.3398500, 6.3575520, 6.3750394, 6.3923174, 6.4093909, 6.4262648, 6.4429435,
 6.4594316, 6.4757334, 6.4918531, 6.5077946, 6.5235620, 6.5391588, 6.5545889, 6.5698556,
 6.5849625, 6.5999128, 6.6147098, 6.6293566, 6.6438562, 6.6582115, 6.6724253, 6.6865005,
 6.7004397, 6.7142455, 6.7279205, 6.7414670, 6.7548875, 6.7681843, 6.7813597, 6.7944159,
 6.8073549, 6.8201790, 6.8328900, 6.8454901, 6.8579810, 6.8703647, 6.8826430, 6.8948178,
 6.9068906, 6.9188632, 6.9307373, 6.9425145, 6.9541963, 6.9657843, 6.9772799, 6.9886847,
 7.0000000, 7.0112273, 7.0223678, 7.0334230, 7.0443941, 7.0552824, 7.0660892, 7.0768156,
 7.0874628, 7.0980321, 7.1085245, 7.1189411, 7.1292830, 7.1395514, 7.1497471, 7.1598713,
 7.1699250, 7.1799091, 7.1898246, 7.1996723, 7.2094534, 7.2191685, 7.2288187, 7.2384047,
 7.2479275, 7.2573878, 7.2667865, 7.2761244, 7.2854022, 7.2946207, 7.3037807, 7.3128830,
 7.3219281, 7.3309169, 7.3398500, 7.3487282, 7.3575520, 7.3663222, 7.3750394, 7.3837043,
 7.3923174, 7.4008794, 7.4093909, 7.4178525, 7.4262648, 7.4346282, 7.4429435, 7.4512111,
 7.4594316, 7.4676056, 7.4757334, 7.4838158, 7.4918531, 7.4998459, 7.5077946, 7.5156998,
 7.5235620, 7.5313815, 7.5391588, 7.5468945, 7.5545889, 7.5622424, 7.5698556, 7.5774288,
 7.5849625, 7.5924570, 7.5999128, 7.6073303, 7.6147098, 7.6220518, 7.6293566, 7.6366246,
 7.6438562, 7.6510517, 7.6582115, 7.6653359, 7.6724253, 7.6794801, 7.6865005, 7.6934870,
 7.7004397, 7.7073591, 7.7142455, 7.7210992, 7.7279205, 7.7347096, 7.7414670, 7.7481928,
 7.7548875, 7.7615512, 7.7681843, 7.7747871, 7.7813597, 7.7879026, 7.7944159, 7.8008999,
 7.8073549, 7.8137812, 7.8201790, 7.8265485, 7.8328900, 7.8392038, 7.8454901, 7.8517490,
 7.8579810, 7.8641861, 7.8703647, 7.8765169, 7.8826430, 7.8887432, 7.8948178, 7.9008668,
 7.9068906, 7.9128893, 7.9188632, 7.9248125, 7.9307373, 7.9366379, 7.9425145, 7.9483672,
 7.9541963, 7.9600019, 7.9657843, 7.9715436, 7.9772799, 7.9829936, 7.9886847, 7.9943534 
};

#if __STD_C
size_t vclogi(size_t v)
#else
size_t vclogi(v)
size_t	v;
#endif
{
	register size_t	lg, b;

	switch(sizeof(v) ) /* find the left most non-zero byte */
	{ default: /* in case sizeof(unsigned int) > 4 */
		for(lg = (sizeof(v)-1)*8; lg >= 32; lg -= 8)
			   if((b = v >> lg) != 0 )
				goto got_byte;
	  case 4: lg = 24; if((b = v >> 24) != 0)
				goto got_byte;
	  case 3: lg = 16; if((b = v >> 16) != 0)
				goto got_byte;
	  case 2: lg =  8; if((b = v >>  8) != 0)
				goto got_byte;
	  case 1: lg =  0; b = v;
	}

got_byte: return lg + _Vcpow2[b];
}


#if __STD_C
double vclog(size_t v)
#else
double vclog(v)
size_t	v;
#endif
{
	register int	lg, b;

	switch(sizeof(v) ) /* find the left most non-zero byte */
	{ default: /* in case sizeof(unsigned int) > 4 */
		for(lg = (sizeof(v)-1)*8; lg >= 32; lg -= 8)
			   if((b = v >> lg) != 0 )
				goto got_byte;
	  case 4: lg = 24; if((b = v >> 24) != 0)
				goto got_byte;
	  case 3: lg = 16; if((b = v >> 16) != 0)
				goto got_byte;
	  case 2: lg =  8; if((b = v >>  8) != 0)
				goto got_byte;
	  case 1: return v == 0 ? 0. : _Vclog2[v];
	}

got_byte:
	lg = lg + _Vcpow2[b] - 7;

#if PRECISE
#define LOG2	0.6931472 /* natural log of 2 */
	b = (1 << lg) - 1;
	return lg + _Vclog2[v >> lg] +  (v & b)/((v & ~b)*LOG2);
#else
	return lg + _Vclog2[v >> lg];
#endif
}


#ifdef TABLE /* code to generate the above _Vc* tables */
#include <stdio.h>
#include <math.h>
main()
{
	int	i, k, lg, cnt;

	printf("static unsigned char _Vcpow2[256] = {\n");
	for(i = 2, k = 0, lg = 0, cnt = 0; i <= 256; i <<= 1, lg += 1)
	{	for(; k < i; ++k)
		{	printf(" %d%c", lg, (k < 255 ? ',' : ' ') );
			if((cnt += 1)%16 == 0 )
				printf("\n");
		}
	}
	printf("};\n\n");

	printf("static double _Vclog2[256] = {\n");
	for(i = 0, cnt = 0; i < 256; ++i)
	{	printf(" %9.7lf%c", (i == 0 ? 0 : log((double)i)/log(2.)),
				   (i < 255 ? ',' : ' ') );
		if((cnt += 1)%8 == 0)
			printf("\n");
	}
	printf("};\n\n");
}
#endif