fix.cpp

Go to the documentation of this file.

/*
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.
*/

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


#define FIX_USE_64BIT_MUL

#define FIX_UNROLL_DIV_LOOP


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;
    }