Exponential function calculation algorithm

Introduction

In my previous article, I introduced a fast algorithm for calculating natural logarithm. Now let's take a look atAlgorithm. As we know, the exponential function ex can be expanded to the tyleseries:

This series converges to all real numbers x, and converges faster when x approaches zero.

C # Program

The following C # program can be written to add an exp Extension Method for the decimal data type according to the tyleseries expansion of ex described earlier:

 1   Using  System;  2   3   Namespace  Skyiv. Extensions  4   {  5     Static  Class  Decimalextensions  6   {  7       Static   Readonly   Decimal Expmax = 66.542129333754749704054283659 m  ;  8       Static   Readonly   Int [] Mask = { 1 , 2 , 4 , 8 , 16 , 32 , 64  };  9       Static   Readonly   Decimal [] Exps = 10   {  11         2.71828182845904523536028747135 m , // Exp (1)  12         7.38905609893065022723042746058 m , //  Exp (2)  13         54.5981500331442390781102612029 m , //  Exp (4)  14         2980.95798704172827474359209945 m , //  Exp (8)  15         8886110.52050787263676302374078 m , //  Exp (16)  16         78962960182680.6951609780226351 m , //  Exp (32)  17         6235149080811616882909238708.93 m    //  Exp (64)  18   };  19       20       Public   Static  Decimal Exp ( This   Decimal  X)  21   {  22         If (X> expmax) Throw   New Overflowexception ( "  Overflow  "  );  23         If (X <-66 ) Return   0  ;  24         VaR N = ( Int ) Decimal  . Round (X );  25         If (N> 66 ) N -- ;  26         Decimal Z = 1, y = exponential (X-N );  27         For ( Int M = (n < 0 )? -N: N, I = 0 ; I <mask. length; I ++ )  28           If (M & Mask [I])! = 0 ) Z * = Exps [I];  29         Return (N < 0 )? (Y/z): (y *Z );  30   }  31       32       Static   Decimal Exponential ( Decimal  Q)  33 { //  Q (almost) in [-0.5, 0.5]  34         Decimal Y = 1 , T = Q; 35         For ( VaR I = 1 ; T! = 0 ; T * = Q/++ I) Y + = T;  36         Return  Y;  37   }  38   }  39 }

Brief description:

1. The expmax value of Row 1 is an approximate value of the natural logarithm of decimal. maxvalue. It is used to detect whether the exp method overflows (Row 3 ).
2. The exp extension method between 20th and 30 rows is used to calculate the exponential function.
3. This method uses the formula ex + y = exey to divide the X parameter into integer part N and decimal part X-N for calculation.
4. The integer N is decomposed into some sums in the numbers 1, 2, 4, 8, 16, 32, and 64, which are calculated using the previously calculated constants.
5. Row 25th is used to prevent e66.5421 from being broken down into a e67e-0.4579, which overflows when e67 is computed. Instead, it is decomposed into e66e0.5421.
6. The exponential method of rows 32nd to 37 uses the Taylor series to calculate the ex. The closer the parameter q is to zero, the faster the calculation is.
7. This algorithm is still very fast. The number of for loop executions of 35th rows will not exceed 22.

Test procedure

The following is a test program that calls the decimal data type exp extension method:

 1   Using  System;  2   Using  Skyiv. extensions;  3   4  Class  Tester  5   {  6     Static   Void  Main ()  7   {  8       Try  9   {  10         Foreach ( VaR XIn   New   Decimal  [] {  11 - 100 ,- 66 ,- 65 ,- 1 , 0 , 1 , 2.5 m , 16 , 66.5421 m , 67 })  12 Console. writeline ( "  {0,-30}: exp ({1 })  "  , X. exp (), X );  13   }  14       Catch  (Exception ex) {console. writeline (ex. Message );}  15   }  16 }

The running result is as follows:

Work $DMCS tester. CS decimalextensions. CSWork $Mono tester.exe0: exp (-100) 0.0000000000000000000000000000: exp (-66) 0.0000000000000000000000000001: exp (-65) 0.3678794411714423215955237702: exp (-1) 1: exp (0) 2.7182818284590452353602874714: exp (1) 12.182493960703473438070175950: exp (2.5) 8886110.520507872636763023741: exp (16) 79225838488862236701995526357: exp (66.5421) Overflow
 

We can see that the e67 overflow is detected. This is because:

• Decimal. maxvalue = 79,228,162,514,264,337,593,543,950,335
• E67 = 125,236,317,084,221,378,051,352,196,074.4365767534885274...

It can be seen that e67 has exceeded the maximum value of decimal.

Constants calculated in advance

The values of E1, E2, E4, E8, e16, E32, and e64 are stored in the exps static read-only array in rows 9th to 18 of the decimalextensions. CS program. How are these values obtained? This is simple. in Linux, there is a high-precision calculator BC. You can edit the following text file exps_in.txt:

Scale = 30E (1) E (2) E (4) E (8) E (16) E (32) E (64) L (2 ^ 96-1) quit
 

E Indicates exp, L indicates ln, and 296-1 indicates decimal. maxvalue. Run the following command:

 
Work $BC-l exps_in.txt> exps_out.txt
 

The following output exps_out.txt can be obtained:

 
Bytes
 

A little bit of sorting can be used in the above C # program:

• The first seven rows are the power of E.
• The last row is the natural logarithm of decimal. maxvalue.

References

1. Wikipedia: exponential function
2. Wikipedia: Taylor series
3. Linux man pages: BC-an arbitrary precision Calculator Language