Not bad! But: if I want to compute a pow(), in which cases it might be ok to get an error up to 25%? Could you provide us some examples?

In Java use BigInteger and BigDecimal:
http://javadoc.ankerl.com/results.html?cx=006156709672261707051%3Atrd_lbkwu2y&cof=FORID%3A9&q=bigdecimal&sa=Search#1356

I forgot to mention that we are using C++ and the g++ compiler.

Клёво, мне понравилось!

Hi,

Your aproximation seems to be really usefor for slow embeddede devices.Besides,I have no RAM free space enough to load math library. However,i tried to used your aproximation for pow in an embedded micro(JN5139) 32-bit RISC CPU using Big-Endian format and using gcc compiler and it doesn’t work. Any idea of what kind of changes can be done to make it work ?
Thanks.

Well, code looks OK, but I’ve Cut and Paste it into J2ME (on Windows XP, SP3),
run it like pow(1.026f, 0.077f) an the result is: 0.97423… !?
The result is [b]less[/b] than 1.0 !.

Here is some code I’ve written to complement what I have available in my J2ME libraries…..related to pow, log, etc.. There is some stuff commented out here, so don’t worry about it. The rest is good. Note I didn’t include the concept of bittians for fast trig function lookups to 8 bits accuracy (great for graphics). Good luck with these.

   /**
     * Fast integer to the power of another integer. Note that if the exponent is < 1
     * then the function returns 0, as you'd expect for an integer function. This
     * function relies on the JVM to catch overflow, so if the int operations just
     * wrap around, then this function will not return the correct value for very large
     * values.
     *
     * @param   x   int that is to be raised to a power.
     * @param   e   int that is the exponent.
     * @return      x^e calculated as an int.
     */
    public static int pow(int x, int e)
    {
        if (e >= 1;    // next bit
        }
        return r;
    }   // pow(int,int)

    /**
     * Another pow function that raises a float to the power of an int. This
     * is a common operation. This uses the same algorithm as the int,int function
     * but with a float as the value to be raised. As with the int algorithm, the
     * exponent must be positive. If e < 0, it will return 0.
     *
     * @param   x   float value to be raised to a power.
     * @param   e   int exponent. e < 0 returns 0.
     */
    public static float pow(float x, int e)
    {
        if (e >= 1;    // next bit
        }
        return r;
    }   // pow(float,int)

    /* TODO: Later debug...not working right now.
     * Fast floating point (single precision) power function. x^e, where both
     * x and e are floats. This is based on the principle that log(x^e) = elog(x).
     * It doesn't matter what base the log is in. If we use 2 then, log2 can be very
     * fast as float is already stored as M*2^E, where M and E are the mantissa and
     * exponent respectively. So, using some math.
     *
     *          a       =   M*2^E
     *  =>      log2(a) =   log2(M) + E
     *
     *  Since M is stored as 1.xxxx (IEEE format), then log2(M) ~ 0
     *
     *  =>      log2(a) =   E approx, with log2(M) as error factor
     *
     *  In J2ME, the Float.floatToIntBits and Float.intBitsToFloat
     *  functions can be used to get floats as IEEE binary format and vice-versa.
     *  So, this allows the exponent and Mantissa parts to be extracted.
     *  As mentioned above, a = 2^Eold * Mold
     *  =>  log2(a) = Eold + log2(Mold)         [Equation 2]
     *
     *  Repeatedly squaring the mantissa, Mold, we get the sequence…
     *      Xnew = Mold^N, where 1.0 <= Xnew   Mold^N = 2^Enew * Mnew
     *  =>  N*log2(Mold) = Enew + log2(Mnew)
     *  =>  log2(Mold) = Enew/N + log2(Mnew)/N
     *
     *  This last equation means that we need to follow the following steps.
     *  1. Extract Eold and Mold from the input value x
     *  2. Put Eold aside and the repeatedly square Mold say 7 times (N = 2^7 = 128) == Xnew
     *  3. Break Xnew into Enew and Mnew
     *  4. Concatenate Enew and the high 9 bits of Mnew (including the invisible first bit).
     *     indicated as [Enew:1.Mnew] and we only use the top 2 bits of Mnew rounded up from 3.
     *     This is explained in http://focus.ti.com/lit/an/spra218/spra218.pdf
     *  5. Normalize this into a float, using the value in step 4 as a new mantissa = log2(Mold)
     *  6. Use Equation 2 above to calculate log2(Xold) = Eold + log2(Mold) {from step 6}
     *  7. This is the answer to return.
     *
     *
     * The squaring can only be 7 times for single precision numbers due to the exponent size,
     * also explained in that article. Note that this is accurate to about 5 decimal places
     * after the decimal point.
     *
     * Note that it might be faster to add the original Eold to Mold as an estimate to a
     * converging algorithm. Note also that other log bases can be calculated by this equality.
     * logn(x) = log2(x)/log2(n)
     *
     * It might also be possible to make this algorithm faster using the float to int
     * techniques used in flog2. Requires experimentation at a future date. Chebychev
     * Polynomials are also worth looking at.
     *
     * @param   f   float to get log to the base 2 of (log2(f)).
     * @return      float returned that holds the value of log2(f).
     * @see         flog2
     *
     */
    /*
    public static float fplog2(float f)
    {
        int Xold = Float.floatToIntBits(f);    // Convert to IEEE binary

        // Extract Eold and Mold from Xold
        int Eold = ((Xold >> 23) & 0xff) - 128;     // Extract the exponent.
        int i = Xold & 0x807fffff;                  // Change the exponent to 0 (127)
        i |= 0x3f800000;                            // and calculate the mantissa, Mold
        float Mold = Float.intBitsToFloat(i);  // Store Mold as a float
        // Mold = ((-1.0/3 * Mold + 2) * Mold = 2.0/3.0;
        // For low precision faster results, can use Eold + Mold here to return
        // not a bad result (3 decimal places or so).

        // Square Mold 7 times
        int N = 128;                                    // N = 2^7 = 128 (7-bit exponent)
        float temp = pow(Mold, N);                  // Xnew = Mold^N, N = 128
        int Xnew = Float.floatToIntBits(temp); // ….in binary format

        //
        //  seeeeeeeemmmmmmmmmmmmmmmmmmmmmmm
        //  33222222222211111111110000000000
        //  10987654321098765432109876543210
        //  12345678EEEEEEE1MMM0100010001000
        //

        // Extract Enew and Mnew from Xnew
        int Enew = ((Xnew >> 23) & 0xff) - 128;     // Extract Exponent.
        i = Xnew & 0x807fffff;                      // Clear exponent part
        i = (i >> 7) | 0x00010000;                  // Place mantissa in [16:0] including hidden bit
        i += (Enew <>= 2;                                    // worst case is if bit 24 got set due to rounding
        while((i & 0x00800000) != 0)                // shift until this bit is set
        {
            i <<= 1;    // shift left
            ++e;        // increment exponent
        }
        i &= 0xff7fffff;                            // Clear bit 23, the hidden mantissa bit
        i |= e << 23;                               // Put the exponent in there.
        temp = Float.intBitsToFloat(i);             // Convert to log2(Mold) value
        temp += (float)Eold;                        // Calculate Eold + log2(Mold) = Xnew = log2(f)

        return (temp);                              // return log2(f)
    };  // fplog2(float
    */

    /**
     * Log to the base 2 of a single precision floating point number. This is a
     * fast algorithm exchanging speed for precision. Is accurate to about 3 digits.
     *
     * This function converts the IEEE floating point input to bits as a 32-bit
     * integer. IEEE FP is sign bit, 8 exponent bits and 23 mantissa bits. The
     * exponent is e + 127 so that negative exponents can be supported. Finally,
     * the mantissa is of the format 1.mmmmmmm where the 1 is assumed and the
     * mmmm’s are the low 23 bits. After every floating point operation, the
     * number is normalized by shifting and adjusting the exponent until the
     * desired 1.mmmmmm mantissa is achieved.
     *
     * So, the input FP is converted to an int, is divided by 2^23 to put the exp to the LHS
     * of the decimal and the mantissa to the RHS of the decimal. 127 is subtracted from it
     * to give the real exponent. Next, the mantissa  is extracted and passed through an
     * error polynomial to converge on the result and added to the origninal ee.mmmmmmm
     * to give a fairly accurate estimate.
     *
     * @param   f       float input value.
     * @return          float result of log2(f)
     * @see             fpow2
     */
    public static float flog2(float f)
    {
        float Xnorm = (float)Float.floatToIntBits(f); // Convert to integer
        Xnorm = (Xnorm * 1.1920929e-7f) - 127;      // eemmmmmmm *(1/2^23)-127 = ee.mmmmm
        float frac = modf(Xnorm);                   // Get RHS of decimal = mantissa
        frac = 0.346607f * frac * (1 - frac);       // Polynomial error calculation
        return Xnorm + frac;                        // New value
    }   // flog2

    /**
     * This is a high performance to a few decimal places for 2^x, where x
     * is a float.
     *
     * It’s reasonable to assume that the input value to be raised to the
     * power is no greater than 128. At too high a value, the power will be
     * larger than the largest number representable by a single precision float
     * and since the exponent is from -127 to 128, then it’s reasonable to
     * assume the input argument is limited to this. This means that the
     * input argument is always in the format iii.fffffff and that due to
     * floating point normalization in IEEE numbers, we can represent this
     * closer to eee.mmmmmmm, where eee is the exponent and mmmmmmm is the
     * mantissa. This allows us to greatly speed up our algorithm.
     *
     * @param   x   float input value to calculate alog2 of (alog2(x) == pow2(x)).
     * @return      float result = 2^x.
     * @see         flog2
     */
    public static float fpow2(float x)
    {
        float frac = modf(x);                       // Get fractional part
        frac = 0.33971f * frac * (1 - frac);        // Polynomial error calculation
        float Xnorm = x + 127 - frac;               // Offset exponent and frac adjust
        Xnorm *= 8388608.0f;                               // * 2^23 (shift left) = eemmmmmmm
        Xnorm = Float.intBitsToFloat((int)Xnorm);   // Convert to proper float
        return Xnorm;
    }   // fpow2

    /**
     * Returns the log to the base 2 of an integer. This is done by converting to an IEEE
     * single precision floating point number, casting this to an int to get the bits and
     * shifting it right 23 bits to seperate the exponent from the mantissa. Then the
     * exponent is normalized by subtracting 128 from it. Finally, if the MS bit after the 1.
     * in the mantissa is a 1, then the exponent is rounded up by 1. The idea is to return the
     * nearest value to the log rounded to nearest.
     *
     * @param   x       int to get the log to the base 2 of.
     * @return  int log2(x) rounded to the nearest int.
     */
    public static int log2(int x)
    {
        int y = Float.floatToIntBits((float)x);
        //int y = (int)Float.intBitsToFloat(x);
        y = (((y & 0x00400000) != 0) ? ((y >> 23) + 1) : (y >> 23)) - 127;
        return y;
    }

    /**
     * Raise input int, x, to the power of 2. Used integer pow function.
     *
     * @param   x       Int to raise to power of 2.
     * @return  2^x returned as int.
     * @see             pow
     */
    public static int pow2(int x)
    {
        return pow(2, x);
    }

    // inv sqrt of float has a trick. It’s not magic.
    //
    // Algorithm fast 1/sqrt(x)
    // f -> int
    //
    // int >> 2 divides exponent by 2
    // if exponent is even mantissa 1.mmmm -> 1.0mmmm
    // if exponent is odd mantissa 1.mmmm -> 1.1mmmm
    //
    // Last step is very clever. sqrt of exp does result id divide by 2.
    // sqrt of 1.x is 1.0x or 1.1x approximately depending on odd or even exponent. APPROX!!!!!!
    // Thing sqrt(2) is similar to sqrt(200), sqrt(20) is similar to sqrt(2000), so this is true.
    //
    // The subtract is handling the 1/x operation.
    //
    // The eqt is 1/sqrt(x) = a - sqrt(x) This is a subtract not in fp numbers but in ints don’t
    // forget…
    //
    // 0x5f3759df -> + exp = BE = 190 - 128 = 62 decimal. Mantissa = 1.01101110101100111011111
    // = 1.375 x 2^62 approx.
    //
    // If exp of sqrt is 0 -> 0x80…..turns exp from BE to 3E or -66 decimal. The mantissa also gets
    // affected. The 1.0 part remains the same (since it’s hidden). The RHS of the decimal point gets
    // subtracted.
    //
    // So finally, convert int back to fp number.
    // Run one pass of newtons sqrt.
/*
float InvSqrt (float x)
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}
*/
    //
    //
    /**
     * Get square root of int, x, and return the result as an int. This is based on
     * an algorithm for the microchip pic 8-bit micro extended to twice the width
     * and optimized for java rather than assembler.
     *
     * @param   x   Value to get square root of. Must be > 0. There is no error check
     *              for negative x.
     * @result      int sqrt(x).
     *
     */
    public static int sqrt(int x)
    {
        int bx = 0xc000;        // Used to shift or set a bit.
        int v = 0x8000;         // sqrt value initial estimate
        int sq; // TODO: was unsigned int (doesn’t exist in java)
        while(bx > 0)           // Loops 16 times through
        {
            if ((sq = v*v) == x)    // check for exact match.
                break;
            v = (sq >= 1;
        }
        return v;               // Return the closest estimate
    }   // sqrt