Hi everybody!
I'm writing a C program for a PIC18F microcontroller.
I need to calculate a power function, in which both base and exponent
are fixed point numbers (ex: 3.15^1.13).
Using pow() function is too expensive...
Is there another way to do that?
It doesn't seem obvious to me. I guess you would want a break down
like:
two_pow_fromIM ( y * two_log_toIM ( x ) );
The idea would be that two_log_toIM and two_pow_fromIM could be
implemented as a scaling (normalize to the range 1 <= x < 2) then
either a post or pre-shift along with a table look up if the
resolution was small enough (and possibly perform interpolations). To
_fromIM and _toIM reflect the fact you might like to convert it to a
temporarily higher resolution intermediate value, or range corrected
for the particular input values.
I am not aware of any really good approximations to log() or 2exp()
except for taylor series or rational function approximations, which
will end up doing no better than using pow() directly. This table
based stuff would obviously compromise accuracy/resolution.
The logarithm can also be calculated effectively using the AGM.
Carlson's (1972) algorithm is described in Math. Comp. 26(118):
543--549.
It may be a good choice since very limited precision is needed, so I
suspect that only a few iterations would be required.
There is a thread on it from 2002 in and
the title of which is:
"Speaking of natural log and AGM... I've gone either blind or mad or
both."
Given an effective logarithm calculation, it can be inverted to create
exp() via something like this idea:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
long double f( long double x, long double k )
{
return expl( x ) - k;
}
int main(int argc, char **argv)
{
long double x;
long double k;
long double lnx;
long double y;
char pszString[256] = "2";
if ( argc < 2 )
{
printf( "Enter a number to find exponential of:" );
fgets( pszString, sizeof( pszString ), stdin );
puts( pszString );
}
else
strncpy( pszString, argv[1], sizeof( pszString ) );
k = atof( pszString );
x = expl( k ) * 1.00001; /* form estimate correct to 5 digits */
lnx = logl(x);
printf("answer = %30.20LeL\nx = %30.20LeL\n", expl( k ), x - x*lnx + x*k );
y = lnx - k;
x -= 3. * x * ( 2 + 3 * y + 2 * y * y ) * y /
( 6 + 12 * y + 11 * y * y + 6 * y * y * y );
printf("answer = %30.20LeL\nx = %30.20LeL\n", expl( k ), x );
return 0;
}
Another possiblity is to use cubic splines. That's why I asked about
memory.
I briefly looked at Carlson's Algorithm as you suggested, and it does
indeed look interesting. Its always good to have an algorithm whose
accuracy can simply be increased *incrementally* by running it through
more and more iterations. This means making good multi-precision
implementations become viable. Thanks for the pointer!
But in general using splines or other forms of interpolation provide
maybe a few more bits of accuracy at most. So its not that useful
except in improving initial guess before running it through
iterations, or in low-precision situations (which may be suitable for
the OP.)
My notion on the spline (strictly for the OP) was to use a spline to
create an approximate exp() function. It looked from his examples as
though he only needs 2 digits of precision + some integral size (not
sure what), and so a spline could be used (potentially) but it may
also be necessary to do range reduction (which is why I asked about
the extreme values that are possible).
It's the same notion and the Numerical Recipes' using a spline to
approximate a functions (I modified xsplint to include exp()):
E:\nr\c\ansi\recipes>cl xsplint.c splint.c spline.c nrutil.c
Microsoft (R) 32-bit C/C++ Optimizing Compiler Version 14.00.50727.762
for 80x86
Copyright (C) Microsoft Corporation. All rights reserved.
xsplint.c
splint.c
spline.c
nrutil.c
Generating Code...
Microsoft (R) Incremental Linker Version 8.00.50727.762
Copyright (C) Microsoft Corporation. All rights reserved.
/out:xsplint.exe
xsplint.obj
splint.obj
spline.obj
nrutil.obj
E:\nr\c\ansi\recipes>xsplint
sine function from 0 to pi
x f(x) interpolation
0.157080 0.156434 0.156351
0.471239 0.453990 0.453981
0.785398 0.707107 0.707088
1.099557 0.891007 0.890984
1.413717 0.987688 0.987663
1.727876 0.987688 0.987663
2.042035 0.891007 0.890983
2.356194 0.707107 0.707089
2.670354 0.453991 0.453978
2.984513 0.156435 0.156433
***********************************
Press RETURN
exponential function from 0 to 1
x f(x) interpolation
0.050000 1.051271 1.051268
0.150000 1.161834 1.161834
0.250000 1.284025 1.284025
0.350000 1.419068 1.419067
0.450000 1.568312 1.568312
0.550000 1.733253 1.733253
0.650000 1.915541 1.915540
0.750000 2.117000 2.116999
0.850000 2.339647 2.339646
0.950000 2.585710 2.585709
***********************************
Press RETURN
The use of a cubic spline means that only 3 multiplies are needed to
estimate the exp() function.