Fast Algorithm for calculating natural logarithm

Introduction

In 1982, tateaki. Sasaki and yasumasa Kanada published a paper: practically Fast Multiple-Precision Evaluation of log (X ). In this four-page paper, they introduced a quick way to calculate the natural logarithm.Algorithm.

C # Program

We know that the system. Math. Log method in the. NET Framework class library can only calculate the natural logarithm of the dobule data type. Now we implement a log extension method for the decimal data type, as shown below:

 1   Using  System;  2   3   Namespace  Skyiv. Extensions  4   {  5    Static   Class  Decimalextensions  6   {  7       Static   Readonly   Decimal PI2 = 6.283185307179586476925286766559 m  ;  8       Static   Readonly   Decimal Ln10 =2.30258509299404568401799145468 m  ;  9       Static   Readonly   Decimal Eps2 = 0.00000000000001 m ; //  1e-14  10       Static   Readonly   Decimal Eps1 = eps2 * eps2; //  1e-28 11       12       Public   Static   Decimal SQRT ( This   Decimal  X)  13   {  14         If (X < 0 ) Throw   New Argumentexception ( " Must be nonnegative  "  );  15         If (X = 0 ) Return   0  ;  16         VaR Y = ( Decimal ) Math. SQRT (( Double  ) X );  17         Return (Y + x/y )/2  ;  18   }  19       20       Public   Static   Decimal Log10 ( This   Decimal  X)  21   {  22         Return Log (X )/ Ln10; 23   }  24       25       Public   Static   Decimal Log ( This   Decimal  X)  26   {  27         If (X <= 0 ) Throw   New Argumentexception ( "  Must be positive  "  );  28         If (X = 1 ) Return   0  ;  29         VaR K = 0  ;  30         For (; X>0.1 m ; K ++) x/= 10  ;  31         For (; X <= 0.01 m ; K --) x * = 10  ;  32         Return K * ln10- Negativelog (X );  33   }  34       35       Static  Decimal Negativelog ( Decimal  Q)  36 { //  Q in (0.01, 0.1]  37         Decimal R = Q, S = Q, n = Q, q2 = Q * q, Q1 = Q2 * Q;  38         For ( VaR P = True ; (N * = Q1)> eps1; S + = N, Q1 * = Q2) 39 R + = (P =! P )? N :- N;  40         Decimal U = 1 - 2 * R, V = 1 + 2 * S, t = u/ V;  41         Decimal A = 1 , B = SQRT ( 1 -T * T ); 42         For (; A-B> eps2; B = SQRT (A * B), a = T) t = (a + B )/ 2  ;  43         Return PI2/(A + B)/V/ V;  44   }  45   }  46 }

In the above program:

• The SQRT extension method is implemented in rows 12th to 18. The Newton iteration method is used to iterate only once (17th rows ). This is because the accuracy of the initial values given by the math. SQRT (double) method in row 16th has reached 15, while that of decimal is 28, and iteration is sufficient.
• Lines 25th to 33 implement the log extension method. This method first transforms the parameter X to the (0.01, 0.1] range (29th to 31 rows), and then calls the negativelog method for calculation.
• The negativelog method of rows 35th to 44 is the core of our algorithm. We use the elliptic θ function (37th to 40 rows) and the arithmetic geometric mean (41st to 43 rows) to quickly calculate the natural logarithm. For the principles of this algorithm, see [1].
• The first row is the 2π value calculated in advance and used in the second row.
• The first row is the value of ln10 calculated in advance, which is used in the second row and the second row.
• The value of ln10 is used in our program, and the value of ln2 is used in reference [1. This is because decimal is better in decimal. The binary type is better for the double type.
• Eps1 = 10-28 and eps2 = 10-14 of rows 10th and 9th use decimal control to calculate accuracy. In reference [1], binary control is used to control the computing accuracy, that is, 2-P and 2-P/2.

Algorithm performance

Add some debugging statements to the above program to report the values of some important variables when the program is running. The results are as follows:

Q: 0.1000000000000000000000000000 K: 2R: 0.0999000009999999000000001000 S: 0.1000000000000000000001000u: 0.8001999980000001999999998000 V: 1.2002000020000002000000002000 loop1: 4b0: 0.895769668060699748813036941a: 0000: 0.9471678269798660936897543331 loop2: 3nl: pushed: 0.0100000000000000000000000001 K: 28r: 0.0099999900000000010000000001 S: 0.01000000000000000000000001u: 0.9800000199999999979999999998 V: 1.0200000200000000020000000002 loop1: 2b0: pushed: 0.65000079528724023734005530455b: 0.6556979528723956496570849678 loop2: 4nl: pushed: 59.867212417845187784467777833
 

The above is the debugging result of calling the LOG method to calculate the natural logarithm for 10 and 100000000000000000000000001 respectively:

• K is the value of the log method when it is executed to 32nd rows. Other variables are the value when the negativelog method is executed in rows 35th to 44.
• Loop1 is the number of for loop executions between 38th and 39 rows.
• Loop2 is the number of for loop executions of 42nd rows.
• Q, R, S, U, V, A, and B are all values at the end of the calculation.
• B0 is the initial value of B (41st rows ). The initial value of A is 1.
• NL is the return value of the negativelog method, that is, the opposite number of the natural logarithm of the parameter q.
• Ln is the return value of the log method, that is, the natural logarithm of parameter X.
• The above debugging results show that this algorithm is very efficient. The number of two for loop executions in the algorithm generally does not exceed 4.

Test procedure

The following is a test program that calls the SQRT, log, and log10 extension methods of the decimal data type:

 1   Using  System;  2   Using Skyiv. extensions;  3   4   Class  Tester  5   {  6     Static   Void  Main ()  7   {  8       Foreach ( VaR X In  New   Decimal [] { 4 / Decimal  . Maxvalue,  9         0.0000001 m , 0.1 m , 1 , 10 , 100000000 , Decimal  . Maxvalue })  10   { 11 Console. writeline ( "  X:  " + X );  12 Console. writeline ( "  SQRT:  " + X. SQRT ());  13 Console. writeline ( "  LN:  " + X. Log ()); 14 Console. writeline ( "  LG:  " + X. log10 ());  15   Console. writeline ();  16   }  17   }  18 }

The running result is as follows:

Work $DMCS tester. CS decimalextensions. CSWork $Mono tester.exeX: 0.0000000000000000000000000001 SQRT: 0.20.0000000001ln:-weight:-28.000000000000000000000000000x: 0.0000001 SQRT: Weight:-weight: 0.1 SQRT: Weight:-weight: 1 SQRT: 1ln: 0lg: 0x: 10 SQRT: Protocol: 1.20.000000000000000000000001x: 100000000 SQRT: 0000ln: bandwidth: 8.20.00000000000000000000000x: 79228162514264337593543950335 SQRT: 281474976710656.00000000000000ln: Protocol: 28.898879583742194740518933893
 

It can be seen that the calculation accuracy of this algorithm can reach 27, and only the last bit may have an error.

References

1. cinii articles: practically Fast Multiple-Precision Evaluation of log (x)
2. Wikipedia: natural logarithm
3. Wikipedia: arithmetic-geometric mean
4. Wikipedia: Newton's method
5. msdn: decimal structure (system)