How to build math library

For those who are interested on how to build themselves their own math library here are a few ideas that I have been testing for CalcReg.

The exponential

The first functions needed is the exponential. It has to be as precise as possible with regards to the speed requested of the program.

It is possible to calculate apart from the program the functions and then go and search for the data in a big file. But this will need a big amount of memory at the end... but why not.

Other possibility is to find a reasonably fast way to calculate the function.

First possibility is the sum of the Taylor series: see http://en.wikipedia.org/wiki/Exponential_function

Of course you do not take the serie up to the infinite but you could get far enough to get reasonable precision. this solution is time consumming for the processor.

Another solution (back to the memory) is to put in memory the exp(1/128)*exp(1/128)*...*exp(1/128) to build the different exp(x=n*128). You therefore need only the exact value of exp(1/128) to start with. You'll understand that if you considere that exp(2x)=exp(x)*exp(x).

The last solution that I present here is to decompose a float number into a UInt32 (an int on 32 bits, could be 64 as well).

float x; //the number to get the exponential

UInt32 a;

Now you will remarque one thing: the binary number is in the base 2, that is to say:

a = b_32 * 2^(32)+b_16 * 2^(16) +...+ b_1 * 2 + b_0. Where b_32 to b_0 get values 1 or 0 only.

So get this:

exp(a) = exp(b_32 * 2^(32)+b_16 * 2^(16) +...+ b_1 * 2 + b_0) = exp(b_32 * 2^(32))*exp(b_16 * 2^(16))*...*exp( b_1 * 2) * exp(b_0)

We then get the exponential result as a product of the exp(b_i * 2^ i)= exp(2^i) or 1 (if b_i is nul)

put a_i =exp(2^i) or 1 (if b_i is nul) then the result is exp(a)=a_32*a_16*..*a_2*a_1*a_0

If you want to get the exponential of a number x from 1/128 to 128 with a precision of 1/128 on x then just multiply a by 0x100 (hexa) at start: a=(int)(x * 0x100); that'll do.

In this method, you just need a set of exact values: exp(1/128),exp(1/64), exp(1/32),... exp(0), exp(1), exp(2),...exp(64),exp(128)

And in the End, to get data in between the discret values, say between x1=1/128 and x2=2*1/128, you need to interpolate with a Lagrange polynomial for example with 3 successive values P(x1), P(x2) and P(x3=3*1/128). see the link http://en.wikipedia.org/wiki/Lagrange_interpolation.

 

The Neperien logarithm

 The ln(x) function can be obtained from the Taylor series (very long calculation for good precision), or from a zero finding algorithm of the function f(t)=exp(t)-x. you get ln(x)=t when f(t)=0; The speed of the calculation depends on the zero finding algorithm and the speed of response of the exponential.

The square root

The square root can be calculated with the Taylor series (as long as ln(x) almost),  or by using the job you've made above with sqrt(x)=exp(0.5*ln(x)), ... or with another zero finding algorithm with f(t)= t*t-x, where the solution t=sqrt(x);

The power function

The Power function can be obtained with a basic x^2=x*x ,...x^n=x*..*x (n times) if n is a relative entire Z (divide for n<0), or through x^n=exp( n*ln(x) )... you understand the importance to have a very precise exponential and a good and fast x->ln(x) function.

 

Sine and Cosine functions can be obtained with the same method as the exponential one. That is to say we use a set of cos(1/2048), cos (1/1024), ... cos(1), cos (2), cos(4). We stop at cos (4) because  for values between 0 and 2pi it is enough. We could call that method the 2-Base method. Because it uses the binary system.

Any number x can be decomposed into a n-set of 2^n. It is the binary system. With a(n) coefficients with 0 or 1 binary values. We therefore use the maths identities cos (a+b)= cosa.cosb - sina.sinb, and sin(a+b)=sinacosb + sinb cosa.

With an algorithm we can calculate cos(an.2^n +...+a0) and sin(an.2^n +...+a0). It is this method that is used in Calcreg. It gives fast and exact result.

First cos1=cos( a0 + 2.a1) and sin1=sin (a0 + 2.a1), then cos2=cos1.cos(2^2.a2)-sin1.sin(2^2.a2), sin2=sin1.cos(2^2.a2)+sin(2^2.a2).cos1, same method for cos3, sin3 and... etc... until you've done the whole sum of the x 2-based development. You then get cosn=cos(x) and sinn=sin(x).

 

Créer un site gratuit avec e-monsite - Signaler un contenu illicite sur ce site