Power by squaring for negative exponents

by

Integer examples are for 32 bit int arithmetics, DWORD is 32bit unsigned int

  1. floating pow(x,y)=x^y

Is usually evaluated like this:

so the fractional exponent can be evaluated: pow(x,y) = exp2(y*log2(x)). This can be done also on fixed point:

  1. integer pow(a,b)=a^b where a>=0 , b>=0

This is easy (you already have that) done by squaring:

        DWORD powuu(DWORD a,DWORD b)
            {   
            int i,bits=32;
            DWORD d=1;
            for (i=0;i<bits;i++)
                {
                d*=d;
                if (DWORD(b&0x80000000)) d*=a;
                b<<=1;
                }
            return d;
            }
  1. integer pow(a,b)=a^b where b>=0

Just add few ifs to handle the negative a

        int powiu(int a,DWORD b)
         {
         int sig=0,c;
         if ((a<0)&&(DWORD(b&1)) { sig=1; a=-a; } // negative output only if a<0 and b is odd
         c=powuu(a,b); if (sig) c=-c;
         return c;
         }
  1. integer pow(a,b)=a^b

So if b<0 then it means 1/powiu(a,-b) As you can see the result is not integer at all so either ignore this case or return floating value or add a multiplier variable (so you can evaluate PI equations on pure Integer arithmetics). This is float result:

        float powfii(int a,int b)
         {
         if (b<0) return 1.0/float(powiu(a,-b));
         else return powiu(a,b);
         }
  1. integer pow(a,b)=a^b where b is fractional

You can do something like this a^(1/bb) where bb is integer. In reality this is rooting so you can use binary search to evaluate:

  • a^(1/2) is square root(a)
  • a^(1/bb) is bb_root(a)

so do a binary search for c from MSB to LSB and evaluate if pow(c,bb)<=a then leave the bit as is else clear it. This is sqrt example:

        int bits(DWORD p) // count how many bits is p
            {
            DWORD m=0x80000000; int b=32;
            for (;m;m>>=1,b--)
             if (p>=m) break;
            return b;
            }
        
        DWORD sqrt(const DWORD &x)
            {
            DWORD m,a;
            m=(bits(x)>>1);
            if (m) m=1<<m; else m=1;
            for (a=0;m;m>>=1) { a|=m; if (a*a>x) a^=m; }
            return a;
            }

so now just change the if (a*a>x) with if (pow(a,bb)>x) where bb=1/b … so b is fractional exponent you looking for and bb is integer. Also m is the number of bits of the result so change m=(bits(x)>>1); to m=(bits(x)/bb);

[edit1] fixed point sqrt example

//---------------------------------------------------------------------------
const int _fx32_fract=16;       // fractional bits count
const int _fx32_one  =1<<_fx32_fract;
DWORD fx32_mul(const DWORD &x,const DWORD &y)   // unsigned fixed point mul
    {
    DWORD a=x,b=y;              // asm has access only to local variables
    asm {                       // compute (a*b)>>_fx32_fract
        mov eax,a               // eax=a
        mov ebx,b               // ebx=b
        mul eax,ebx             // (edx,eax)=eax*ebx
        mov ebx,_fx32_one
        div ebx                 // eax=(edx,eax)>>_fx32_fract
        mov a,eax;
        }
    return a;
    }
DWORD fx32_sqrt(const DWORD &x) // unsigned fixed point sqrt
    {
    DWORD m,a;
    if (!x) return 0;
    m=bits(x);                  // integer bits
    if (m>_fx32_fract) m-=_fx32_fract; else m=0;
    m>>=1;                      // sqrt integer result is half of x integer bits
    m=_fx32_one<<m;             // MSB of result mask
    for (a=0;m;m>>=1)           // test bits from MSB to 0
        {
        a|=m;                   // bit set
        if (fx32_mul(a,a)>x)    // if result is too big
         a^=m;                  // bit clear
        }
    return a;
    }
//---------------------------------------------------------------------------

so this is unsigned fixed point. High 16 bits are integer and low 16 bits are fractional part.

  • this is fp -> fx conversion: DWORD(float(x)*float(_fx32_one))
  • this is fp <- fx conversion: float(DWORD(x))/float(_fx32_one))
  • fx32_mul(x,y) is x*y it uses assembler of 80386+ 32bit architecture (you can rewrite it to karatsuba or whatever else to be platform independent)
  • fx32_sqrt(x) is sqrt(x)

In fixed point you should be aware of the fractional bit shift for multiplication: (a<<16)*(b<<16)=(a*b<<32) you need to shift back by >>16 to get result (a*b<<16). Also the result can overflow 32 bit therefore I use 64 bit result in assembly.

[edit2] 32bit signed fixed point pow C++ example

When you put all the previous steps together you should have something like this:

//---------------------------------------------------------------------------
//--- 32bit signed fixed point format (2os complement)
//---------------------------------------------------------------------------
// |MSB              LSB|
// |integer|.|fractional|
//---------------------------------------------------------------------------
const int _fx32_bits=32;                                // all bits count
const int _fx32_fract_bits=16;                          // fractional bits count
const int _fx32_integ_bits=_fx32_bits-_fx32_fract_bits; // integer bits count
//---------------------------------------------------------------------------
const int _fx32_one       =1<<_fx32_fract_bits;         // constant=1.0 (fixed point)
const float _fx32_onef    =_fx32_one;                   // constant=1.0 (floating point)
const int _fx32_fract_mask=_fx32_one-1;                 // fractional bits mask
const int _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask; // integer bits mask
const int _fx32_sMSB_mask =1<<(_fx32_bits-1);           // max signed bit mask
const int _fx32_uMSB_mask =1<<(_fx32_bits-2);           // max unsigned bit mask
//---------------------------------------------------------------------------
float fx32_get(int   x) { return float(x)/_fx32_onef; }
int   fx32_set(float x) { return int(float(x*_fx32_onef)); }
//---------------------------------------------------------------------------
int fx32_mul(const int &x,const int &y) // x*y
    {
    int a=x,b=y;                // asm has access only to local variables
    asm {                       // compute (a*b)>>_fx32_fract
        mov eax,a
        mov ebx,b
        mul eax,ebx             // (edx,eax)=a*b
        mov ebx,_fx32_one
        div ebx                 // eax=(a*b)>>_fx32_fract
        mov a,eax;
        }
    return a;
    }
//---------------------------------------------------------------------------
int fx32_div(const int &x,const int &y) // x/y
    {
    int a=x,b=y;                // asm has access only to local variables
    asm {                       // compute (a*b)>>_fx32_fract
        mov eax,a
        mov ebx,_fx32_one
        mul eax,ebx             // (edx,eax)=a<<_fx32_fract
        mov ebx,b
        div ebx                 // eax=(a<<_fx32_fract)/b
        mov a,eax;
        }
    return a;
    }
//---------------------------------------------------------------------------
int fx32_abs_sqrt(int x)            // |x|^(0.5)
    {
    int m,a;
    if (!x) return 0;
    if (x<0) x=-x;
    m=bits(x);                  // integer bits
    for (a=x,m=0;a;a>>=1,m++);  // count all bits
    m-=_fx32_fract_bits;        // compute result integer bits (half of x integer bits)
    if (m<0) m=0; m>>=1;
    m=_fx32_one<<m;             // MSB of result mask
    for (a=0;m;m>>=1)           // test bits from MSB to 0
        {
        a|=m;                   // bit set
        if (fx32_mul(a,a)>x)    // if result is too big
         a^=m;                  // bit clear
        }
    return a;
    }
//---------------------------------------------------------------------------
int fx32_pow(int x,int y)       // x^y
    {
    // handle special cases
    if (!y) return _fx32_one;                           // x^0 = 1
    if (!x) return 0;                                   // 0^y = 0  if y!=0
    if (y==-_fx32_one) return fx32_div(_fx32_one,x);    // x^-1 = 1/x
    if (y==+_fx32_one) return x;                        // x^+1 = x
    int m,a,b,_y; int sx,sy;
    // handle the signs
    sx=0; if (x<0) { sx=1; x=-x; }
    sy=0; if (y<0) { sy=1; y=-y; }
    _y=y&_fx32_fract_mask;      // _y fractional part of exponent
     y=y&_fx32_integ_mask;      //  y integer part of exponent
    a=_fx32_one;                // ini result
    // powering by squaring x^y
    if (y)
        {
        for (m=_fx32_uMSB_mask;(m>_fx32_one)&&(m>y);m>>=1);     // find mask of highest bit of exponent
        for (;m>=_fx32_one;m>>=1)
            {
            a=fx32_mul(a,a);
            if (int(y&m)) a=fx32_mul(a,x);
            }
        }
    // powering by rooting x^_y
    if (_y)
        {
        for (b=x,m=_fx32_one>>1;m;m>>=1)                            // use only fractional part
            {
            b=fx32_abs_sqrt(b);
            if (int(_y&m)) a=fx32_mul(a,b);
            }
        }
    // handle signs
    if (sy) { if (a) a=fx32_div(_fx32_one,a); else a=0; /*Error*/ }     // underflow
    if (sx) { if (_y) a=0; /*Error*/ else if(int(y&_fx32_one)) a=-a; }  // negative number ^ non integer exponent, here could add test if 1/_y is integer instead
    return a;
    }
//---------------------------------------------------------------------------

I have tested it like this:

float a,b,c0,c1,d;
int x,y;
for (a=0.0,x=fx32_set(a);a<=10.0;a+=0.1,x=fx32_set(a))
 for (b=-2.5,y=fx32_set(b);b<=2.5;b+=0.1,y=fx32_set(b))
    {
    if (!x) continue; // math pow has problems with this
    if (!y) continue; // math pow has problems with this
    c0=pow(a,b);
    c1=fx32_get(fx32_pow(x,y));
    d=0.0;
    if (fabs(c1)<1e-3) d=c1-c0; else d=(c0/c1)-1.0;
    if (fabs(d)>0.1)
     d=d; // here add breakpoint to check inconsistencies with math pow
    }
  • a,b are floating point
  • x,y are closest fixed point representations of a,b
  • c0 is math pow result
  • c1 is fx32_pow result
  • d is difference

hope did not forget something trivial but it seems like it works properly. Do not forget that fixed point has very limited precision so the results will differ a bit …

P.S. Take a look at this:

More Related Contents:

  1. Recursive C function
  2. The most efficient way to implement an integer based power function pow(int, int)
  3. Explain this snippet which finds the maximum of two integers without using if-else or any other comparison operator?
  4. How do you print the EXACT value of a floating point number?
  5. C – How to implement Set data structure?
  6. Gaussian elimination without result for acceleration
  7. Linked list recursive reverse
  8. signed and unsigned arithmetic implementation on x86
  9. A possible algorithm for determining whether two strings are anagrams of one another? [closed]
  10. Algorithm to pick values from array that sum closest to a target value?
  11. Recursive and non-recursive traversal of three degree tree
  12. why does this recursive program for factorial not work? [closed]
  13. how to show time performance for merge sort? [closed]
  14. Sum of the all the numbers between a and b [closed]
  15. What integer hash function are good that accepts an integer hash key?
  16. All Possible Combinations of a list of Values
  17. Calculating factorial of large numbers in C
  18. Rounding integer division (instead of truncating)
  19. Combine Gyroscope and Accelerometer Data
  20. C How to “draw” a Binary Tree to the console [closed]
  21. How to do unsigned saturating addition in C?
  22. Traverse Matrix in Diagonal strips
  23. unexpected output in C (recursion)
  24. Are floating point operations in C associative?
  25. Why does gdb evaluate sqrt(3) to 0?
  26. Inverting a 4×4 matrix
  27. To the power of in C? [duplicate]
  28. How does this bitwise operation check for a power of 2?
  29. Conditional tests in primality by trial division
  30. QuickSort and Hoare Partition

Leave a Comment