/*
This program is distributed under the terms of the 'MIT license'. The text
of this licence follows...

Copyright (c) 2004-2005 J.D.Medhurst (a.k.a. Tixy)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/

/**
@file

@brief 32-bit fixed-point maths routines
*/

#include "common.h"
#include "fix.h"


/**
If this macro is defined, the implementation of the multiply methods
will use 64bit arithametic. This produces much more efficient code
on some architecture/compiler combinations, like ARM+GCC.
*/
#define FIX_USE_64BIT_MUL

/**
If this macro is defined, the implementation of the Fix::Div will
use an unrolled loop. This may produce faster code.
*/
#define FIX_UNROLL_DIV_LOOP


/**
Helper function for interpolatating lookup tables used by transcendental functions.

If the lookup table represents the function <i>f(x)</i> over the range
<i>0</i> to <i>N<<shift</i>, then it must contain entries for:
<i>{ f(-1<<shift), f(0), f(1<<shift), f(2<<shift), ... f((N-1)<<shift)), f(N<<shift), f((N+1)<<shift)) }</i>.
This function then returns <i>f(value)</i>.

Notes:

- Each entry in the table must have a value greater than the preceding one.
- This function produces undefined results due to arithmatic overflow
when the difference between any adjacent table entries is close to the maximum value
of an <code>int</code> shifted right by \a shift. (The exact conditions which cause
overflow are more complex than this, therefore it is advised that any
table is tested for correctness using all posible values of \a value.)

@param table	The lookup table.
@param value	The value to lookup in \a table.
@param shift	The number of least significant bits in \a value which are
				discarded when indexing into \a table.

@return The interpolated result of looking up \a value in \a table.
*/
static inline int Interpolate(const int* table, int value, int shift)
	{
	// get values from table (required value lies between b and c)
	int i = value>>shift;
	int a = table[i+0];
	int b = table[i+1];
	int c = table[i+2];
	int d = table[i+3];

	// interpolate
	int f = value&((1<<shift)-1); // f = factional part of index
	int cadb = ( (c-a) - (d-b) ) >> 2;
	int r = (c-b) + cadb - (cadb*f>>shift);
	r = (unsigned)(r*f)>>shift;  // this cast to unsigned assumes table only contains increasing values (and gains us 1 more bit headroom)
	r += b;

	return r;
	}


/*
Members of class Fix
*/


fix Fix::Add(fix a,fix b)
	{
	fix r = a+b;
	if((~(a^b) & (a^r)) < 0)
		r = (r>>31)^0x80000000;  // produce saturated result if overflow
	return r;
	}


fix Fix::Sub(fix a,fix b)
	{
	fix r = a-b;
	if(((a^b) & (a^r)) < 0)
		r = (r>>31)^0x80000000;  // produce saturated result if overflow
	return r;
	}


fix Fix::Mul(fix a,fix b)
	{
#ifdef FIX_USE_64BIT_MUL

	int64_t r = ((int64_t)a*(int64_t)b);
	r += 0x8000;
	int32_t hi = r>>32;
	r >>= 16;
	if(hi>=0x8000)
		return 0x7fffffff;
	if(hi<-0x8000)
		return 0x80000000;
	return (int32_t)r;

#else

	// calculate sign result
	int sign = a^b;

	// calculate absolute values
	if(a<0) a=-a;
	if(b<0) b=-b;

	// do the multiply in four parts
	uint32_t al = a&0xFFFF;
	uint32_t bl = b&0xFFFF;
	uint32_t c = al*bl;
	c += 0x8000;
	c >>= 16;
	uint32_t c1 = bl*((uint32_t)a>>16);
	c += c1;
	if(c>=c1) // No carry from addition
		{
		uint32_t c2 = al*((uint32_t)b>>16);
		c += c2;
		if(c>=c2) // No carry from addition
			{
			uint32_t d = ((uint32_t)a>>16)*((uint32_t)b>>16);
			if(d<0x10000) // No overflow from multiply
				{
				uint32_t dl = (uint32_t)d<<16;
				c += dl;
				if(c>=dl) // No carry from addition
					{
					if(sign<0)
						{
						if(c<=0x80000000)
							return -(int32_t)c;
						}
					else
						{
						if(c<=0x7fffffff)
							return c;
						}
					}
				}
			}
		}

	// produce saturated result
	return (sign<0) ? 0x80000000 : 0x7fffffff;

#endif
	}


fix Fix::MulNS(fix a,fix b)
	{
#ifdef FIX_USE_64BIT_MUL
	int64_t r = ((int64_t)a*(int64_t)b);
	r += 0x8000;
	r >>= 16;
	return (int32_t)r;
#else
	uint32_t al = a&0xFFFF;
	uint32_t bl = b&0xFFFF;
	uint32_t r = al*bl;
	r += 0x8000;
	r >>= 16;
	r += bl*(a>>16);
	r += al*(b>>16);
	r += ((a>>16)*(b>>16))<<16;
	return r;
#endif
	}


fix Fix::Div(fix a,fix b)
	{
	// calculate sign bit of result
	int32_t r = a^b;
	r &= (1<<31);

	// calculate absolute values
	if(a<0) a = -a;
	if(b<0) b = -b;

	// calculate the number of integer bits is the result
	int32_t intBits = 0;
	uint32_t q = a;
	if((uint32_t)b<=((uint32_t)q>>16))
		intBits += 16, q >>= 16;
	if((uint32_t)b<=((uint32_t)q>>8))
		intBits += 8,  q >>= 8;
	if((uint32_t)b<=((uint32_t)q>>4))
		intBits += 4,  q >>= 4;
	if((uint32_t)b<=((uint32_t)q>>2))
		intBits += 2,  q >>= 2;
	if((uint32_t)b<=((uint32_t)q>>1))
		intBits += 1,  q >>= 1;

#ifdef FIX_UNROLL_DIV_LOOP

	// calculate the integer part of result (bits 31 to 16)
	switch(intBits)
		{
		case 14:
			if((uint32_t)a>=((uint32_t)b<<14))
				a -= (b<<14), r += 1<<(14+16);
		case 13:
			if((uint32_t)a>=((uint32_t)b<<13))
				a -= (b<<13), r += 1<<(13+16);
		case 12:
			if((uint32_t)a>=((uint32_t)b<<12))
				a -= (b<<12), r += 1<<(12+16);
		case 11:
			if((uint32_t)a>=((uint32_t)b<<11))
				a -= (b<<11), r += 1<<(11+16);
		case 10:
			if((uint32_t)a>=((uint32_t)b<<10))
				a -= (b<<10), r += 1<<(10+16);
		case 9:
			if((uint32_t)a>=((uint32_t)b<<9))
				a -= (b<<9), r += 1<<(9+16);
		case 8:
			if((uint32_t)a>=((uint32_t)b<<8))
				a -= (b<<8), r += 1<<(8+16);
		case 7:
			if((uint32_t)a>=((uint32_t)b<<7))
				a -= (b<<7), r += 1<<(7+16);
		case 6:
			if((uint32_t)a>=((uint32_t)b<<6))
				a -= (b<<6), r += 1<<(6+16);
		case 5:
			if((uint32_t)a>=((uint32_t)b<<5))
				a -= (b<<5), r += 1<<(5+16);
		case 4:
			if((uint32_t)a>=((uint32_t)b<<4))
				a -= (b<<4), r += 1<<(4+16);
		case 3:
			if((uint32_t)a>=((uint32_t)b<<3))
				a -= (b<<3), r += 1<<(3+16);
		case 2:
			if((uint32_t)a>=((uint32_t)b<<2))
				a -= (b<<2), r += 1<<(2+16);
		case 1:
			if((uint32_t)a>=((uint32_t)b<<1))
				a -= (b<<1), r += 1<<(1+16);
		case 0:
			if((uint32_t)a>=((uint32_t)b<<0))
				a -= (b<<0), r += 1<<(0+16);
			break;
	default:
		// produce saturated result
		return (r<0) ? 0x80000000 : 0x7fffffff; // saturated result
		}

	// calculate the fractional part of result (bits 15 to 0)
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<15;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<14;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<13;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<12;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<11;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<10;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<9;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<8;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<7;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<6;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<5;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<4;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<3;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<2;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<1;
	a <<= 1;
	if((uint32_t)a>=(uint32_t)b)
		a -= b, r += 1<<0;

#else

	if(intBits>14)
		return (r<0) ? 0x80000000 : 0x7fffffff; // saturated result

	// calculate the integer part of result (bits 31 to 16)
	b <<= intBits;
	int32_t bit = 0x10000<<intBits;
	while(bit>0x10000)
		{
		if((uint32_t)a>=(uint32_t)b)
			{
			a -= b;
			r += bit;
			}
		b >>= 1;
		bit >>= 1;
		}
	if((uint32_t)a>=(uint32_t)b)
		{
		a -= b;
		r += bit;
		}
	bit >>= 1;

	// calculate the fractional part of result (bits 15 to 0)
	do
		{
		a <<= 1;
		if((uint32_t)a>=(uint32_t)b)
			{
			a -= b;
			r += bit;
			}
		bit >>= 1;
		}
	while(bit);

#endif

	return (r<0) ? 0x80000000-r : r;

	}


fix Fix::Sqrt(ufix a)
	{
	ufix r,nr,m;

	// calculate integer part (bits 31 to 16)
	r = 0;
	m = 0x40000000;
	do
		{
		nr = r+m;
		if(nr<=a)
			{
			a -= nr;
			r = nr+m;
			}
		r >>= 1;
		m >>= 2;
		}
	while(m);

	// calculate bits 15 to 8
	r <<= 8;
	a <<= 8;
	m = 0x40;
	do
		{
		nr = r+m;
		if(nr<=a)
			{
			a -= nr;
			r = nr+m;
			}
		r >>= 1;
		m >>= 2;
		}
	while(m);

	// calculate bits 7 to 0
	r <<= 8;
	a <<= 8;
	m = 0x40;
	do
		{
		nr = r+m;
		if(nr<=a)
			{
			a -= nr;
			r = nr+m;
			}
		r >>= 1;
		m >>= 2;
		}
	while(m);

	// round result
	if(r<a)
		r++;

	return r;
	}


fix Fix::Log2(ufix a)
	{
	// trap 0
	if(a==0)
		return (fix)-0x80000000;

	// calculate integer part of result in i
	// and set n = normalised value of a
	int32_t i=15*0x10000;
	uint32_t n=a;
	if(n<(uint32_t)(1<<(32-16)))
		n <<= 16, i -= 16*0x10000;
	if(n<(uint32_t)(1<<(32-8)))
		n <<= 8, i -= 8*0x10000;
	if(n<(uint32_t)(1<<(32-4)))
		n <<= 4, i -= 4*0x10000;
	if(n<(uint32_t)(1<<(32-2)))
		n <<= 2, i -= 2*0x10000;
	if(n<(uint32_t)(1<<(32-1)))
		n <<= 1, i -= 1*0x10000;

	// reduce n to the 23 most significant bits and clear
	// the most significant bit, leaving a 22 bit value in n
	n = (n-0x80000000+(1<<8))>>9;

	// calculate fractional part of result (in f) by interpolation of lookup table values
	static const int32_t LogTable[] =
		{
		0xffff45e1,
		0x00000000,0x0000b73d,0x00016bad,0x00021d67,0x0002cc7f,0x00037908,0x00042316,0x0004caba,
		0x00057007,0x0006130b,0x0006b3d8,0x0007527c,0x0007ef06,0x00088984,0x00092204,0x0009b892,
		0x000a4d3c,0x000ae00d,0x000b7111,0x000c0053,0x000c8ddd,0x000d19bb,0x000da3f6,0x000e2c98,
		0x000eb3aa,0x000f3935,0x000fbd43,0x00103fdb,0x0010c105,0x001140ca,0x0011bf31,0x00123c42,
		0x0012b803,0x0013327c,0x0013abb4,0x001423b0,0x00149a78,0x00151012,0x00158482,0x0015f7d0,
		0x00166a01,0x0016db19,0x00174b20,0x0017ba19,0x0018280a,0x001894f7,0x001900e6,0x00196bdb,
		0x0019d5da,0x001a3ee8,0x001aa709,0x001b0e41,0x001b7495,0x001bda07,0x001c3e9d,0x001ca259,
		0x001d053f,0x001d6754,0x001dc89a,0x001e2914,0x001e88c7,0x001ee7b4,0x001f45e1,0x001fa34e,
		0x00200000,0x00205bf9,0x0020b73d
		};
	int32_t f = Interpolate(LogTable,n,16);

	// scale result to 16 bits (from the 22 bit precision used in lookup table)
	f = (f+(1<<(5-1)))>>5;

	// return sum of integer and fractional part
	return i+f;
	}


ufix Fix::Exp2(fix a)
	{
	if(a>=0x00100000)
		return 0xffffffffu; // result will be too big so result max value

	// special cases for small values
	if(a<-0x000ead96)
		{
		if(a<-0x00110000)
			return 0;
		if(a<-0x000f6a3f)
			return 1;
		return 2;
		}

	// table of values for ((2^n)-1)<<32 for n in range 0x0000.0000 to 0x0000.ff00
	static const int32_t Exp2TableFF00[] =
		{
		0x00000000,0x00b1afa5,0x0163da9f,0x02168143,0x02c9a3e7,0x037d42e1,0x04315e86,0x04e5f72f,
		0x059b0d31,0x0650a0e3,0x0706b29d,0x07bd42b7,0x08745187,0x092bdf66,0x09e3ecac,0x0a9c79b1,
		0x0b5586cf,0x0c0f145e,0x0cc922b7,0x0d83b233,0x0e3ec32d,0x0efa55fd,0x0fb66aff,0x1073028d,
		0x11301d01,0x11edbab5,0x12abdc06,0x136a814f,0x1429aaea,0x14e95934,0x15a98c8a,0x166a4547,
		0x172b83c7,0x17ed4869,0x18af9388,0x19726583,0x1a35beb6,0x1af99f81,0x1bbe0840,0x1c82f952,
		0x1d487316,0x1e0e75eb,0x1ed5022f,0x1f9c1843,0x2063b886,0x212be357,0x21f49917,0x22bdda27,
		0x2387a6e7,0x2451ffb8,0x251ce4fb,0x25e85711,0x26b4565e,0x2780e341,0x284dfe1f,0x291ba759,
		0x29e9df51,0x2ab8a66d,0x2b87fd0d,0x2c57e397,0x2d285a6e,0x2df961f6,0x2ecafa93,0x2f9d24ab,
		0x306fe0a3,0x31432ede,0x32170fc4,0x32eb83ba,0x33c08b26,0x3496266e,0x356c55f9,0x36431a2d,
		0x371a7373,0x37f26231,0x38cae6d0,0x39a401b7,0x3a7db34e,0x3b57fbfe,0x3c32dc31,0x3d0e544e,
		0x3dea64c1,0x3ec70df1,0x3fa4504a,0x40822c36,0x4160a21f,0x423fb270,0x431f5d95,0x43ffa3f8,
		0x44e08606,0x45c2042a,0x46a41ed1,0x4786d668,0x486a2b5c,0x494e1e19,0x4a32af0d,0x4b17dea6,
		0x4bfdad53,0x4ce41b81,0x4dcb299f,0x4eb2d81d,0x4f9b2769,0x508417f4,0x516daa2c,0x5257de83,
		0x5342b569,0x542e2f4f,0x551a4ca5,0x56070dde,0x56f4736b,0x57e27dbe,0x58d12d49,0x59c0827f,
		0x5ab07dd4,0x5ba11fba,0x5c9268a5,0x5d845909,0x5e76f15a,0x5f6a320d,0x605e1b97,0x6152ae6c,
		0x6247eb03,0x633dd1d1,0x6434634c,0x652b9feb,0x66238825,0x671c1c70,0x68155d44,0x690f4b19,
		0x6a09e667,0x6b052fa7,0x6c012750,0x6cfdcddd,0x6dfb23c6,0x6ef92985,0x6ff7df95,0x70f7466f,
		0x71f75e8e,0x72f8286e,0x73f9a48a,0x74fbd35d,0x75feb564,0x77024b1a,0x780694fd,0x790b938a,
		0x7a11473e,0x7b17b097,0x7c1ed013,0x7d26a62f,0x7e2f336c,0x7f387849,0x80427543,0x814d2add,
		0x82589994,0x8364c1eb,0x8471a462,0x857f4179,0x868d99b4,0x879cad93,0x88ac7d98,0x89bd0a47,
		0x8ace5422,0x8be05bad,0x8cf3216b,0x8e06a5e0,0x8f1ae991,0x902fed02,0x9145b0b9,0x925c353a,
		0x93737b0c,0x948b82b5,0x95a44cbc,0x96bdd9a7,0x97d829fd,0x98f33e47,0x9a0f170c,0x9b2bb4d5,
		0x9c49182a,0x9d674194,0x9e86319e,0x9fa5e8d0,0xa0c667b5,0xa1e7aed8,0xa309bec4,0xa42c9804,
		0xa5503b23,0xa674a8af,0xa799e133,0xa8bfe53c,0xa9e6b557,0xab0e5213,0xac36bbfd,0xad5ff3a3,
		0xae89f995,0xafb4ce62,0xb0e07298,0xb20ce6c9,0xb33a2b84,0xb468415b,0xb59728de,0xb6c6e29f,
		0xb7f76f2f,0xb928cf22,0xba5b030a,0xbb8e0b79,0xbcc1e904,0xbdf69c3f,0xbf2c25bd,0xc0628614,
		0xc199bdd8,0xc2d1cd9f,0xc40ab5ff,0xc544778f,0xc67f12e5,0xc7ba8898,0xc8f6d940,0xca340575,
		0xcb720dce,0xccb0f2e6,0xcdf0b555,0xcf3155b5,0xd072d4a0,0xd1b532b0,0xd2f87080,0xd43c8eac,
		0xd5818dcf,0xd6c76e86,0xd80e316c,0xd955d71f,0xda9e603d,0xdbe7cd63,0xdd321f30,0xde7d5641,
		0xdfc97337,0xe11676b1,0xe264614f,0xe3b333b1,0xe502ee78,0xe6539246,0xe7a51fbc,0xe8f7977c,
		0xea4afa2a,0xeb9f4867,0xecf482d8,0xee4aaa21,0xefa1bee6,0xf0f9c1cb,0xf252b376,0xf3ac948d,
		0xf50765b6,0xf6632798,0xf7bfdad9,0xf91d8022,0xfa7c1819,0xfbdba369,0xfd3c22b8,0xfe9d96b2
		};

	// table of values for ((2^n)-1)<<40 for n in range 0x0000.0000 to 0x0000.00ff
	static const int32_t Exp2Table00FF[] =
		{
		0x00000000,0x00b17255,0x0162e525,0x02145871,0x02c5cc37,0x03774079,0x0428b535,0x04da2a6d,
		0x058ba01f,0x063d164d,0x06ee8cf5,0x07a00419,0x08517bb7,0x0902f3d1,0x09b46c65,0x0a65e575,
		0x0b175eff,0x0bc8d905,0x0c7a5386,0x0d2bce81,0x0ddd49f8,0x0e8ec5e9,0x0f404256,0x0ff1bf3e,
		0x10a33ca1,0x1154ba7e,0x120638d7,0x12b7b7ab,0x136936fa,0x141ab6c4,0x14cc3709,0x157db7c9,
		0x162f3904,0x16e0baba,0x17923ceb,0x1843bf97,0x18f542be,0x19a6c660,0x1a584a7d,0x1b09cf16,
		0x1bbb5429,0x1c6cd9b7,0x1d1e5fc1,0x1dcfe645,0x1e816d45,0x1f32f4bf,0x1fe47cb5,0x20960526,
		0x21478e11,0x21f91778,0x22aaa15a,0x235c2bb7,0x240db68f,0x24bf41e2,0x2570cdb0,0x262259f9,
		0x26d3e6bd,0x278573fd,0x283701b7,0x28e88fed,0x299a1e9d,0x2a4badc9,0x2afd3d6f,0x2baecd91,
		0x2c605e2e,0x2d11ef46,0x2dc380d9,0x2e7512e7,0x2f26a570,0x2fd83874,0x3089cbf4,0x313b5fee,
		0x31ecf464,0x329e8954,0x33501ec0,0x3401b4a7,0x34b34b09,0x3564e1e6,0x3616793e,0x36c81111,
		0x3779a95f,0x382b4228,0x38dcdb6d,0x398e752c,0x3a400f67,0x3af1aa1d,0x3ba3454e,0x3c54e0fa,
		0x3d067d21,0x3db819c3,0x3e69b6e0,0x3f1b5479,0x3fccf28c,0x407e911b,0x41303025,0x41e1cfaa,
		0x42936faa,0x43451025,0x43f6b11b,0x44a8528d,0x4559f479,0x460b96e1,0x46bd39c3,0x476edd21,
		0x482080fa,0x48d2254e,0x4983ca1e,0x4a356f68,0x4ae7152e,0x4b98bb6e,0x4c4a622a,0x4cfc0961,
		0x4dadb113,0x4e5f5941,0x4f1101e9,0x4fc2ab0d,0x507454ab,0x5125fec5,0x51d7a95a,0x5289546a,
		0x533afff5,0x53ecabfc,0x549e587d,0x5550057a,0x5601b2f2,0x56b360e5,0x57650f53,0x5816be3d,
		0x58c86da1,0x597a1d81,0x5a2bcddc,0x5add7eb2,0x5b8f3003,0x5c40e1cf,0x5cf29417,0x5da446da,
		0x5e55fa17,0x5f07add1,0x5fb96205,0x606b16b4,0x611ccbdf,0x61ce8184,0x628037a5,0x6331ee41,
		0x63e3a559,0x64955ceb,0x654714f9,0x65f8cd82,0x66aa8686,0x675c4005,0x680df9ff,0x68bfb475,
		0x69716f66,0x6a232ad2,0x6ad4e6b9,0x6b86a31b,0x6c385ff9,0x6cea1d52,0x6d9bdb26,0x6e4d9975,
		0x6eff583f,0x6fb11785,0x7062d746,0x71149782,0x71c65839,0x7278196b,0x7329db19,0x73db9d42,
		0x748d5fe6,0x753f2305,0x75f0e6a0,0x76a2aab5,0x77546f46,0x78063452,0x78b7f9da,0x7969bfdc,
		0x7a1b865a,0x7acd4d53,0x7b7f14c7,0x7c30dcb7,0x7ce2a522,0x7d946e07,0x7e463769,0x7ef80145,
		0x7fa9cb9d,0x805b9670,0x810d61be,0x81bf2d87,0x8270f9cc,0x8322c68c,0x83d493c7,0x8486617d,
		0x85382fae,0x85e9fe5b,0x869bcd83,0x874d9d27,0x87ff6d45,0x88b13ddf,0x89630ef4,0x8a14e084,
		0x8ac6b290,0x8b788517,0x8c2a5819,0x8cdc2b96,0x8d8dff8f,0x8e3fd403,0x8ef1a8f2,0x8fa37e5c,
		0x90555442,0x91072aa3,0x91b9017f,0x926ad8d6,0x931cb0a9,0x93ce88f7,0x948061c0,0x95323b05,
		0x95e414c5,0x9695ef00,0x9747c9b6,0x97f9a4e8,0x98ab8095,0x995d5cbd,0x9a0f3961,0x9ac1167f,
		0x9b72f41a,0x9c24d22f,0x9cd6b0c0,0x9d888fcc,0x9e3a6f53,0x9eec4f55,0x9f9e2fd3,0xa05010cc,
		0xa101f241,0xa1b3d430,0xa265b69b,0xa3179982,0xa3c97ce3,0xa47b60c0,0xa52d4519,0xa5df29ec,
		0xa6910f3b,0xa742f505,0xa7f4db4b,0xa8a6c20b,0xa958a947,0xaa0a90ff,0xaabc7932,0xab6e61e0,
		0xac204b09,0xacd234ae,0xad841ece,0xae360969,0xaee7f480,0xaf99e011,0xb04bcc1f,0xb0fdb8a7
		};

	// Treat 'a' as i+j+k where i is integer part of value, j is bits 15 to 8 of
	// the fractional part and k is bits 7 to 0 of the fractional part
	// We calculate 2^a as 2^(i+j+k) which is (2^i)*(2^j)*(2^k)

	int i = a>>16; // i = integer part of value
	uint32_t x=Exp2TableFF00[(a>>8)&0xff];	// x = (2^j-1) << 32
	uint32_t y=Exp2Table00FF[a&0xff]; 		// y = (2^k-1) << 40

	// calculate p = (x*y)>>32 = (2^j-1)(2^k-1)<<40
	uint32_t xl = x&0xffff;
	uint32_t yl = y&0xffff;
	uint32_t xh = x>>16;
	uint32_t yh = y>>16;
	uint32_t pl = (xl*yl+(1<<15))>>16;
	uint32_t pm1 = xh*yl+pl;
	uint32_t pm2 = xl*yh;
	unsigned p = xh*yh;
	uint32_t pm = pm1+pm2;
	if(pm<pm2) p += (1<<16); // check for carry
	if(pm&(1<<15)) p += 1; // round up if required
	p += pm>>16;

	// calculate n = p+x+y = x*y+x+y = (x+1)(y+1)-1 = ((2^j)(2^k)-1)<<32
	unsigned round = ((p&0xff)+(y&0xff)+0x80)>>8;
	unsigned n = x+(y>>8)+round;
	n += p>>8;

	// n = n>>(16-i) = 2^(i-16)*n = ((2^i)*((2^j)*(2^k)-1)<<16
	n >>= 15-i;
	n = (n+1)>>1;

	// finally add 2^i to get result of 2^a
	n += 0x80000000>>(15-i);

	return (ufix)n;
	}


fix Fix::Cos(fixangle angle)
	{
	return Fix::Sin(angle+0x4000);
	}


fix Fix::Sin(fixangle angle)
	{
	// reduce angle to first quadrant
	unsigned n = angle&0x3FFF;
	if(angle&(1<<14))
		n = 0x4000-n;

	// calc Sin through table lookup
	static const int32_t SinTable[] =
		{
		-0x000c8fb3,
		0x00000000,0x000c8fb3,0x001917a7,0x00259021,0x0031f170,0x003e33f3,0x004a5019,0x00563e6a,
		0x0061f78b,0x006d7440,0x0078ad75,0x00839c3d,0x008e39da,0x00987fc0,0x00a26799,0x00abeb4a,
		0x00b504f3,0x00bdaef9,0x00c5e403,0x00cd9f02,0x00d4db31,0x00db941a,0x00e1c598,0x00e76bd8,
		0x00ec835e,0x00f10908,0x00f4fa0b,0x00f853f8,0x00fb14be,0x00fd3aac,0x00fec46d,0x00ffb10f,
		0x01000000,0x00ffb10f,0x00ffb10f
		};
	int r = Interpolate(SinTable,n,9);

	// scale result to 16 bits (from the 24 bit precision used in lookup table)
	r = (r+(1<<(8-1)))>>8;

	// produce correct sign for result
	if(angle&(1<<15))
		r = -r;

	return r;
	}


fix Fix::Tan(fixangle angle)
	{
	// n = angle in first quadrant
	int n = angle&0x7fff;
	bool neg = n>0x4000;
	if(neg) n = 0x8000-n;

	// calculate tangent by interpolation of lookup table values
	unsigned r;
	if(n<=0x2000)
		{
		static const int32_t TanTable0000[] =
			{
			0xffe6dcba,
			0x00000000,0x00192346,0x00324e4f,0x004b88e7,0x0064daee,0x007e4c61,0x0097e564,0x00b1ae4d,
			0x00cbafaf,0x00e5f267,0x01007fa7,0x011b6104,0x0136a083,0x015248ae,0x016e649f,0x018b0019,
			0x01a8279a,0x01c5e872,0x01e450e1,0x02037033,0x022356e3,0x024416be,0x0265c311,0x028870db,
			0x02ac3705,0x02d12ea3,0x02f7733e,0x031f232f,0x03486005,0x03734ef8,0x03a01979,0x03ceedd6,
			0x04000000,0x04338a74
			};
		r = Interpolate(TanTable0000,n-0x0000,8);
		// scale result to 16 bits (from the 26 bit precision used in lookup table)
		r = (r+(1<<(10-1)))>>10;
		}
	else if(n<=0x2d80)
		{
		static const int32_t TanTable2000[] =
			{
			0x01f395da,
			0x02000000,0x020cb920,0x0219c53a,0x02272890,0x0234e7ab,0x02430762,0x02518ce0,0x02607dac,
			0x026fdfb2,0x027fb948,0x0290113f,0x02a0eee8,0x02b25a22,0x02c45b6c,0x02d6fbf1,0x02ea4598,
			0x02fe431c,0x03130022,0x0328894f,0x033eec66,0x0356386c,0x036e7dc8,0x0387ce70,0x03a23e1b,
			0x03bde277,0x03dad368,0x03f92b57,0x04190786,0x043a8876
			};
		r = Interpolate(TanTable2000,n-0x2000,7);
		// scale result to 16 bits (from the 25 bit precision used in lookup table)
		r = (r+(1<<(9-1)))>>9;
		}
	else if(n<=0x35c0)
		{
		static const int32_t TanTable2D80[] =
			{
			0x0204737a,0x020c83c3,0x0214c89d,
			0x021d443b,0x0225f8ef,0x022ee92e,0x02381791,0x024186d9,0x024b39ef,0x025533ea,0x025f7813,
			0x026a09e6,0x0274ed19,0x0280259c,0x028bb7a6,0x0297a7b1,0x02a3fa88,0x02b0b54a,0x02bddd70,
			0x02cb78da,0x02d98dd2,0x02e8231d,0x02f73fff,0x0306ec4e,0x0317307b,0x032815a6,0x0339a5ac,
			0x034beb3d,0x035ef1f2,0x0372c665,0x03877650,0x039d10ad,0x03b3a5d7,0x03cb47bd,0x03e40a0a,
			0x03fe0261
			};
		r = Interpolate(TanTable2D80,n-0x2d80,6);
		// scale result to 16 bits (from the 24 bit precision used in lookup table)
		r = (r+(1<<(8-1)))>>8;
		}
	else if(n<=0x39e0)
		{
		static const int32_t TanTable35C0[] =
			{
			0x01ebc1c4,0x01f20505,0x01f86f02,
			0x01ff0130,0x0205bd16,0x020ca44f,0x0213b88a,0x021afb90,0x02226f3e,0x022a158f,0x0231f097,
			0x023a0288,0x02424db5,0x024ad491,0x025399b6,0x025c9fe3,0x0265ea01,0x026f7b26,0x0279569c,
			0x02837fdc,0x028dfa9d,0x0298cace,0x02a3f4a5,0x02af7c9d,0x02bb6780,0x02c7ba6a,0x02d47ad7,
			0x02e1aea3,0x02ef5c1b,0x02fd8a00,0x030c3f99,0x031b84b8,0x032b61d0,0x033be000,0x034d0923,
			0x035ee7ea
			};
		r = Interpolate(TanTable35C0,n-0x35c0,5);
		// scale result to 16 bits (from the 23 bit precision used in lookup table)
		r = (r+(1<<(7-1)))>>7;
		}
	else if(n<=0x3c70)
		{
		static const int32_t TanTable39E0[] =
			{
			0x01a22f46,0x01a68491,0x01aaf08f,
			0x01af73f5,0x01b40f80,0x01b8c3f6,0x01bd9224,0x01c27ae2,0x01c77f0f,0x01cc9f96,0x01d1dd6a,
			0x01d7398d,0x01dcb509,0x01e250f7,0x01e80e7b,0x01edeec8,0x01f3f321,0x01fa1cd8,0x02006d4d,
			0x0206e5f6,0x020d885a,0x02145612,0x021b50d0,0x02227a5a,0x0229d48f,0x02316169,0x023922fc,
			0x02411b7b,0x02494d38,0x0251baa7,0x025a6660,0x02635323,0x026c83da,0x0275fb9a,0x027fbdac,
			0x0289cd8b,0x02942eec,0x029ee5c0,0x02a9f63c,0x02b564da,0x02c13664,0x02cd6ff9,0x02da1711,
			0x02e7318b
			};
		r = Interpolate(TanTable39E0,n-0x39e0,4);
		// scale result to 16 bits (from the 22 bit precision used in lookup table)
		r = (r+(1<<(6-1)))>>6;
		}
	else if(n<=0x3df0)
		{
		static const int32_t TanTable3C70[] =
			{
			0x0169dabd,0x016d0b89,0x01704ac0,
			0x017398c5,0x0176f600,0x017a62d9,0x017ddfbf,0x01816d24,0x01850b7f,0x0188bb49,0x018c7d04,
			0x01905133,0x01943860,0x0198331a,0x019c41f6,0x01a0658d,0x01a49e82,0x01a8ed7c,0x01ad5328,
			0x01b1d03c,0x01b66576,0x01bb139b,0x01bfdb79,0x01c4bde6,0x01c9bbc2,0x01ced5f9,0x01d40d7d,
			0x01d96350,0x01ded87c,0x01e46e1a,0x01ea254f,0x01efff4d,0x01f5fd57,0x01fc20be,0x02026ae5,
			0x0208dd40,0x020f7955,0x021640bf,0x021d352e,0x0224586a,0x022bac51,0x023332de,0x023aee23,
			0x0242e055,0x024b0bc4,0x025372e6,0x025c1852,0x0264fec8,0x026e2932,0x02779aa6,0x0281566b
			};
		r = Interpolate(TanTable3C70,n-0x3c70,3);
		// scale result to 16 bits (from the 21 bit precision used in lookup table)
		r = (r+(1<<(5-1)))>>5;
		}
	else if(n<=0x3ebc)
		{
		static const int32_t TanTable3DF0[] =
			{
			0x01396c6c,0x013bcd53,0x013e3784,0x0140ab35,0x014328a0,
			0x0145afff,0x0148418d,0x014add89,0x014d8433,0x015035cd,0x0152f29c,0x0155bae5,0x01588ef2,
			0x015b6f0e,0x015e5b87,0x016154ad,0x01645ad4,0x01676e52,0x016a8f7f,0x016dbeb8,0x0170fc5c,
			0x017448cf,0x0177a477,0x017b0fbd,0x017e8b10,0x018216e2,0x0185b3aa,0x018961e2,0x018d220a,
			0x0190f4a6,0x0194da41,0x0198d368,0x019ce0b1,0x01a102b6,0x01a53a18,0x01a9877e,0x01adeb98,
			0x01b26719,0x01b6fac1,0x01bba753,0x01c06d9e,0x01c54e79,0x01ca4ac3,0x01cf6366,0x01d49958,
			0x01d9ed98,0x01df6132,0x01e4f53d,0x01eaaadf,0x01f0834b,0x01f67fc3,0x01fca197,0x0202ea2c,
			0x02095af4
			};
		r = Interpolate(TanTable3DF0,n-0x3df0,2);
		// scale result to 16 bits (from the 20 bit precision used in lookup table)
		r = (r+(1<<(4-1)))>>4;
		}
	else if(n<=0x3f90)
		{
		static const int32_t TanTable3EBC[] =
			{
			0x00ffe079,0x01017516,0x01030eb9,
			0x0104ad7a,0x01065172,0x0107fabb,0x0109a96e,0x010b5da7,0x010d1780,0x010ed715,0x01109c84,
			0x011267ea,0x01143965,0x01161115,0x0117ef18,0x0119d391,0x011bbea1,0x011db06b,0x011fa912,
			0x0121a8bb,0x0123af8b,0x0125bda9,0x0127d33e,0x0129f071,0x012c156d,0x012e425e,0x0130776f,
			0x0132b4cf,0x0134faae,0x0137493b,0x0139a0a9,0x013c012c,0x013e6af8,0x0140de45,0x01435b4b,
			0x0145e245,0x0148736f,0x014b0f06,0x014db54c,0x01506681,0x015322eb,0x0155ead0,0x0158be78,
			0x015b9e30,0x015e8a44,0x01618306,0x016488c8,0x01679be1,0x016abcaa,0x016deb7e,0x017128be,
			0x017474cc,0x0177d00f,0x017b3af1,0x017eb5e0,0x0182414d,0x0185ddb0,0x01898b84,0x018d4b47,
			0x01911d7f,0x019502b5,0x0198fb77,0x019d085b,0x01a129fc,0x01a560f9,0x01a9adfb,0x01ae11b0,
			0x01b28ccd,0x01b72010,0x01bbcc3e,0x01c09224,0x01c5729a,0x01ca6e80,0x01cf86bf,0x01d4bc4c,
			0x01da1028,0x01df835d,0x01e51704,0x01eacc41,0x01f0a448,0x01f6a05b,0x01fcc1cc,0x020309fc,
			0x02097a5f,0x0210147d,0x0216d9f0,0x021dcc68,0x0224edad,0x022c3f9d,0x0233c433,0x023b7d81,
			0x02436dbc,0x024b9734,0x0253fc5f,0x025c9fd4,0x02658453,0x026eacc6,0x02781c43,0x0281d611,
			0x028bddad,0x029636cb,0x02a0e55c,0x02abed94,0x02b753f0,0x02c31d37,0x02cf4e89,0x02dbed5f,
			0x02e8ff97,0x02f68b7b
			};
		r = Interpolate(TanTable3EBC,n-0x3ebc,1);
		// scale result to 16 bits (from the 19 bit precision used in lookup table)
		r = (r+(1<<(3-1)))>>3;
		}
	else if(n<0x4000)
		{
		static const int32_t TanTable3F90[] =
			{
			0x005d1ff3,0x005df6bd,0x005ed16f,0x005fb025,0x006092fa,0x00617a0c,0x0062657b,0x00635566,
			0x006449ed,0x00654334,0x0066415f,0x00674491,0x00684cf3,0x00695aac,0x006a6de7,0x006b86ce,
			0x006ca58f,0x006dca59,0x006ef55e,0x007026d1,0x00715ee9,0x00729ddc,0x0073e3e5,0x00753142,
			0x00768633,0x0077e2fa,0x007947dd,0x007ab526,0x007c2b22,0x007daa20,0x007f3276,0x0080c47c,
			0x0082608e,0x00840710,0x0085b866,0x008774fe,0x00893d49,0x008b11bf,0x008cf2de,0x008ee12b,
			0x0090dd33,0x0092e78b,0x009500d0,0x009729a7,0x009962c0,0x009bacd6,0x009e08af,0x00a0771d,
			0x00a2f8fd,0x00a58f3f,0x00a83add,0x00aafce4,0x00add675,0x00b0c8c1,0x00b3d50f,0x00b6fcbe,
			0x00ba4145,0x00bda438,0x00c12747,0x00c4cc43,0x00c89521,0x00cc83fe,0x00d09b21,0x00d4dd01,
			0x00d94c4b,0x00ddebe4,0x00e2bef2,0x00e7c8e5,0x00ed0d7a,0x00f290c9,0x00f8574c,0x00fe65ee,
			0x0104c218,0x010b71c1,0x01127b80,0x0119e6a4,0x0121bb49,0x012a027b,0x0132c656,0x013c122e,
			0x0145f2c4,0x0150767c,0x015bada6,0x0167aad4,0x0174833b,0x01824f38,0x01912ae6,0x01a136df,
			0x01b2992c,0x01c57e75,0x01da1b7f,0x01f0af1b,0x020984ae,0x0224f778,0x02437703,0x02658d16,
			0x028be5ec,0x02b75bab,0x02e906ce,0x0322561d,0x036532a4,0x03b43743,0x0413099b,0x0486ee3e,
			0x0517cc0b,0x05d20dc8,0x06ca656d,0x08261355,0x0a2f982f,0x0d94caee,0x145f306a,0x28be60d9,
			0xffffffff
			};
		r = TanTable3F90[n-0x3f90];
		}
	else
		{
		// result is infinity so return saturated result
		return angle<0 ? (int32_t)-0x80000000 : (int32_t)0x7fffffff;
		}

	return neg ? -(int32_t)r: (int32_t)r;
	}


fixangle Fix::ACos(fix value)
	{
	if(value<-0x10000)
		value = -0x10000;
	else if(value>0x10000)
		value = 0x10000;
	return 0x4000-Fix::ASin(value);
	}


fixangle Fix::ASin(fix value)
	{
	unsigned n = value;
	if(value<0)
		n = (unsigned)-(int)n;

	// calculare arc-sine by interpolation of lookup table values
	int r;
	if(n<0xc000)
		{
		static const int32_t ASinTable0000[] =
			{
			0xfffeb9ff,
			0x00000000,0x00014601,0x00028c53,0x0003d349,0x00051b37,0x00066474,0x0007af57,0x0008fc3f,
			0x000a4b8d,0x000b9daa,0x000cf305,0x000e4c19,0x000fa967,0x00110b83,0x0012730c,0x0013e0b9,
			0x00155555,0x0016d1cc,0x0018572d,0x0019e6b6,0x001b81e1,0x001d2a76,0x001ee2a7,0x0020ad31,
			0x00228d9c,0x00248890
			};
		r = Interpolate(ASinTable0000,n,11);
		}
	else if(n<0xf200)
		{
		static const int32_t ASinTableC000[] =
			{
																						 0x00221339,
			0x00228d9c,0x002309a5,0x0023876b,0x00240706,0x00248890,0x00250c26,0x002591e7,0x002619f5,
			0x0026a476,0x00273194,0x0027c17d,0x00285465,0x0028ea85,0x00298420,0x002a217e,0x002ac2f5,
			0x002b68e6,0x002c13c1,0x002cc40b,0x002d7a5e,0x002e3777,0x002efc36,0x002fc9b5,0x0030a152,
			0x003184d3,0x0032768f,0x003379c1
			};
		r = Interpolate(ASinTableC000,n-0xc000,9);
		}
	else if(n<0xfe00)
		{
		static const int32_t ASinTableF200[] =
			{
											 0x003238a2,0x0032768f,0x0032b592,0x0032f5ba,0x00333719,
			0x003379c1,0x0033bdc8,0x00340345,0x00344a52,0x0034930c,0x0034dd94,0x00352a10,0x003578a9,
			0x0035c991,0x00361d01,0x0036733a,0x0036cc8b,0x00372953,0x00378a02,0x0037ef25,0x0038596e,
			0x0038c9be,0x00394144,0x0039c19e,0x003a4d1f,0x003ae75a,0x003b9654
			};
		r = Interpolate(ASinTableF200,n-0xF200,7);
		}
	else if(n<=0xffe0)
		{
		static const int32_t ASinTableFE00[] =
			{
																						 0x003ad319,
			0x003ae75a,0x003afbed,0x003b10d5,0x003b2617,0x003b3bb7,0x003b51ba,0x003b6827,0x003b7f02,
			0x003b9654,0x003bae22,0x003bc677,0x003bdf5a,0x003bf8d6,0x003c12f8,0x003c2dcb,0x003c4960,
			0x003c65c7,0x003c8314,0x003ca160,0x003cc0c5,0x003ce165,0x003d0369,0x003d2702,0x003d4c6e,
			0x003d73ff,0x003d9e1e,0x003dcb5f,0x003dfc94,0x003e3300,0x003e70c5,0x003eba0a,0x003f1984
			};
		r = Interpolate(ASinTableFE00,n-0xFE00,4);
		}
	else if(n<=0x10000)
		{
		static const int32_t ASinTableFFE1[] =
			{
					   0x003ebf2c,0x003ec464,0x003ec9b2,0x003ecf17,0x003ed496,0x003eda2f,0x003edfe4,
			0x003ee5b6,0x003eeba8,0x003ef1bb,0x003ef7f2,0x003efe4f,0x003f04d5,0x003f0b88,0x003f126c,
			0x003f1984,0x003f20d5,0x003f2867,0x003f303e,0x003f3865,0x003f40e5,0x003f49c9,0x003f5323,
			0x003f5d06,0x003f678d,0x003f72dc,0x003f7f28,0x003f8cc2,0x003f9c33,0x003fae83,0x003fc661,
			0x00400000,
			};
		r = ASinTableFFE1[n-0xffe1];
		}
	else
		{
		// value out of range, return PI/2
		r = 0x00400000;
		}

	// scale result to 16 bits (from the 24 bit precision used in lookup table)
	r = (r+(1<<(8-1)))>>8;

	// produce correct sign for result
	if(value<0)
		r = -r;

	return r;
	}

fixangle Fix::ATan(fix value)
	{
	unsigned n = value;
	if(value<0)
		n = (unsigned)-(int)n;

	// calculare arc-tangent by interpolation of lookup table values
	static const int32_t ATanTable[] =
		{
		0xfffeba28,
		0x00000000,0x000145d8,0x00028b0d,0x0003cf00,0x00051112,0x000650ad,0x00078d40,0x0008c643,
		0x0009fb38,0x000b2bac,0x000c5734,0x000d7d75,0x000e9e1d,0x000fb8e7,0x0010cd99,0x0011dc04,
		0x0012e405,0x0013e582,0x0014e06a,0x0015d4b7,0x0016c267,0x0017a982,0x00188a16,0x00196434,
		0x001a37f6,0x001b0575,0x001bccd2,0x001c8e2d,0x001d49ab,0x001dff72,0x001eafa7,0x001f5a74,
		0x00200000,0x0020a074,
		};
	int r;
	if(n<=0x10000)
		{
		r = Interpolate(ATanTable,n,11);
		}
	else
		{
		// For n>1 use the identity ATAN(n) = PI/2-ATAN(1/n)
		n = (unsigned)-(int)(n>>1)/(unsigned)n+1;	// n = 1/n
		r = Interpolate(ATanTable,n,11);
		r = 0x400000-r; // r = PI/2-r
		}

	// scale result to 16 bits (from the 24 bit precision used in lookup table)
	r = (r+(1<<(8-1)))>>8;

	if(value<0)
		r = -r;

	return r;
	}


fix Fix::Random(uint32_t& seed)
	{
	return seed=(seed*69069+1);
	}


ufix Fix::Random(uint32_t& seed,ufix range)
	{
	// generate next random number
	ufix n = seed=(seed*69069+1);

	// multiply range by random number (we won't bother with carry from low order half
	// of 64bit result because we don't really need acuracy for random numbers.)
	ufix lo = range&0xffff;
	ufix hi = range>>16;
	ufix r = hi*(n&0xffff)>>16;
	n >>= 16;
	r += lo*n>>16;
	r += hi*n;

	return r;
	}