Some things could be added and improved, for example:
A method to translate an input for the base to the necessary magic constant.
The "result_scaler" can use a larger range, it was done here with a bitshift for speed.
A well constructed loop that can be unrolled at compile time, to allow a speed/accuracy tradeoff.
About the magic constant, with a constant of 1 this will approximate e^x. At the same time the errata within the [-.5,.5) domain will require small adjustment of the constant to minimize total error.
To use a base that isn't a power of two requires modulo in the equation.
Code: Select all
#include <stdio.h>
#include <math.h>
#include <stdint.h>
int main()
{
double in = 1.0; // input should be positive, compensate for input domain with offset
// input domain is (-.5,63.5)
double in_rounded = round(in);
double result_scaler = (double)( ((uint64_t)1) << ((uint64_t)in_rounded) );
double h = 0.693038794439242; // magic multiplier for the base (this is base2)
double x = (in - in_rounded) * h;
double dividend = x+16.0;
double divisor = x-16.0;
double result = (dividend / divisor);
result *= result;
result *= result;
result *= result;
result *= result_scaler;
// multiply by the 1/2^offset that compensates the positive input
printf("\nresult = %.15f", result);
return 0;
}