"Machin" like formulae to determine log(2) Pascal Sebah (Paris) Let arctanh(x) = 1/2*log((1+x)/(1-x)) = x + x^3/3 + x^5 / 5 +....+ x^(2i+1)/(2i+1)+.... and LL(q) = arctanh(1/q) Note that the function arctanh is similar to the arctan function: arctan(x) = arctanh(ix)/i. It is possible to find with this function many formulae similar to those for arctan. We give two useful relations: log(1+1/p)=2*LL(2p+1) p>=1 and log(1+2/q)=2*LL(q+1) q>=1. By convention we define the efficiency E by the sum of the inverse of the common logarithm of the square of the argument of the LL function.The series LL(q) converges like (1/q^2)^n. As a reference Machin formula for Pi as an efficiency of E = 0.9256. _________________________________________________________________ Some examples: With one term: log(2) = 2*LL(3) , E = 1.048 (here p=1) With two terms: log(2) = 4*LL(6)+2*LL(99), E = 0.8931 (derived from 2=(7/5)^2*(50/49) ) log(2) = 4*LL(7)+2*LL(17), E = 0.9980 (derived from 2=(4/3)^2*(9/8) ) With three terms: log(2) = 14*LL(23)+6*LL(65)-4*LL(485), E = 0.8291 log(2) = 14*LL(31)+10*LL(49)+6*LL(161), E = 0.8577 With four terms: log(2) = 62*LL(127)+20*LL(251)+54*LL(449)+24*LL(4801), E = 0.7704 log(2) = 38*LL(52)-8*LL(161)+10*LL(849)+14*LL(70226), E = 0.7918 log(2) = 28*LL(127)+54*LL(161)+34*LL(244)-10*LL(4801), E = 0.8095 With five terms: log(2) = 144*LL(255)+54*LL(449)+24*LL(4801)+20*LL(31751)+82*LL(32257), E = 0.7541 others .... _________________________________________________________________ Many other formulae are available, this is only a small selection. It's possible to perform separately the computation of each LL, making all formulae parallelisable. _________________________________________________________________ _________________________________________________________________ 29 243 200 digits of log(2) Xavier Gourdon (Paris) Who : Xavier Gourdon (email : Xavier.Gourdon@inria.fr) When : 18 / 12 / 1997 Machine used : SGI R10000 with 256 Mo of memory Timing : 20 hours 58 minutes Formula used : A Salamin-Brent like algorithm based on the Arithmetic-Geometric Mean (AGM) sequence (Reference : algorithm 7.2, page 221 of "Pi and the AGM", Borwein and Borwein, John Wiley and Sons) Statistics on the first 29 000 000 digits : Number of 0 : 2902800 Number of 1 : 2898945 Number of 2 : 2899526 Number of 3 : 2896565 Number of 4 : 2899484 Number of 5 : 2897580 Number of 6 : 2898477 Number of 7 : 2903938 Number of 8 : 2902584 Number of 9 : 2900101 Digits from rank 29 243 181 to rank 29 243 200 (for future verification): ...1621261609 7077552756... A verification has been done by computing log(4) with the same algorithm. The verification used a temporary variable computed for log(2) and took 10 hours 30 minutes on the same machine. The first 10 millions digits are verified with previous record of Patrick Demichel. The program consists in using a Fast Fourier Transform modulo two prime numbers of 59 bits, which were especially choosen to improve the computations modulo these primes. A special thank to Pierre Reiss, who wrote for me some basic routines in assembly code. The saving is about 25 % of the total timings. _________________________________________________________________ _________________________________________________________________ 58 486 400 digits of log(2) Xavier Gourdon (Paris) Who : Xavier Gourdon (email : Xavier.Gourdon@inria.fr) When : between the 24/ 12 / 1997 and 28/12/1997 Machine used : SGI R10000 with 256 Mo of memory Real time : 83 hours 30minutes. CPU time : 53 hours. The difference between the time and the CPU time is due to pagination: the available memory was 256 Mo, the memory needed was 350 Mo. Digits from rank 58 486 381 to rank 58 486 400 (for future verification): ...2153347499 3361339978... The same type of verification as the preceeding record has been done and took 41h30 on the same machine. _________________________________________________________________ _________________________________________________________________ comparison of various ways of computing log(2) Patrick Demichel (Paris) Following are the timings of 3 differents way of computing log(2) method1 : log(2) = log(3/2)+log(4/3) used for my previous record method2 : log(2) = 7LL(31) + 5LL(49) + 3LL(161) ( Pascal Sebah method ) method3 : Xavier Gourdon method described above We can see that method3 is far superior to others methods and this starting at ~100000 digits, under this limit the simple code of method2 should be used. decimals computed method1 timing ratio method2 timing ratio method3 timing ratio 28000 6.10s 5.21s 12.27s 56000 24.36s 3.994 20.76s 3.984 26.51s 2.16 112000 97.60s 4.006 82.51s 3.974 60.19s 2.27 224000 392.8s 4.025 332.7s 4.032 133.52s 2.22 448000 1581s 4.025 1348s 4.051 317.08s 2.37 900000 6324s 5392s 848.79s 2.68 1.8M 26000s 22000s 2251s 2.65 3.6M 28h 25h 5824s 2.58 7.3M 112h 100h 14117s 2.42 14.6M 448h 400h 9h10mn 2.33 29.2M 75j 67j 20h58mn 2.28 58.4M 300j 268j 53h 2.52 _________________________________________________________________ _________________________________________________________________ The first method uses the property that log(2)=log(3/2)+log(4/3) This method was originally used by Leonard Euler to compute manually log(2) The way this is coded in the program below is to limit the number of memory accesses . The number of decimals computed = limit*7 _________________________________________________________________ _________________________________________________________________ #include #include #include #include #define limit 524288 const double dbase = 1000*1000*10; void div(double* a1,double* l1,double* a2,double *to,double di) { register double temp1= 0.0; register double temp2= 0.0; register double temp3= 0.0; register double temp4= 0.0; register double temp5= 0.0; register double val1 = 25.0 ; register double val2 = 49.0 ; register double ival1 = 1.0 / val1 ; register double ival2 = 1.0 / val2 ; register double idi = 1.0 / di; register double t1,t2,t3,t4,t5; for (;;){ temp1 = ( temp1 * dbase ) + *a1; t1 = (double) (long) (temp1 * ival1); *a1 = t1; temp1 = temp1 - ( t1 * val1 ); temp2 = ( temp2 * dbase ) + *a2; t2 = (double) (long) (temp2 * ival2); *a2 = t2; temp2 = temp2 - ( t2 * val2 ); temp3 = ( temp3 * dbase ) + t1; t3 = (double) (long) (temp3 * idi); t5 = *to + t3; temp3 = temp3 - ( t3 * di ); temp4 = ( temp4 * dbase ) + t2; t4 = (double) (long) (temp4 * idi); *to = t5 + t4; temp4 = temp4 - ( t4 * di ); if ( ++a1 > l1 ) return; a2++; to++; } } void adjust(double *a,double *l) { double register carry = 0.0; double register idbase = 1.0 / dbase; for(;;){ double register temp = *a + carry; carry = (double) (long) ( temp * idbase ); if ( temp < 0.0 ) carry -= 1.0; *a = temp - ( carry * dbase ); if ( --a < l ) break; } } #define AVOID_TRASHING 63 #define der limit-1 double mat1[limit+AVOID_TRASHING]; double mat2[limit+AVOID_TRASHING]; double mat3[limit+AVOID_TRASHING]; int main() { char out1[128]; char out2[128]; long lp,lp1,lp2,first; for(lp=0;lp= limit-2 ) break; div(mat1_deb,mat1_fin,mat2_deb,mat3_deb,(double) lp1); } mat3_deb=&mat3[0]; mat3_fin=&mat3[der]; adjust(mat3_fin,mat3_deb); const int digits_per_line = 100; strcpy(out1,""); for(lp2=1;lp2= digits_per_line ){ strncpy(out2,out1,digits_per_line); printf("%s\n",out2); strcpy(out1,&out1[digits_per_line]); } } printf("%s\n",out1); return(0); } _________________________________________________________________ _________________________________________________________________ The second method uses the method proposed by Pascal Sebah This method method is more efficient and we can easily parallelize the computation _________________________________________________________________ _________________________________________________________________ #include #include #include #define limit 262144 const double dbase = 1000*1000*10; double add(double* a,double* c,double* l) { double carry; double temp; carry = 0.0; for(;;){ temp = *a + *c + carry; if ( temp >= dbase ){ carry = 1.0; *a = temp - dbase; }else{ carry = 0.0; *a = temp; } a--; if ( a < l ) break; c--; } return(carry); } double mul(double* a,double* l,double val) { double carry; double temp; carry = 0.0; for (;;){ temp = ( *a * val ) + carry; carry = (double) (long) (temp / dbase); *a = temp - ( carry * dbase) ; a--; if ( a < l ) return(carry); } } void adjust(double *a,double *l) { double register carry = 0.0; double register idbase = 1.0 / dbase; for(;;){ double register temp = *a + carry; carry = (double) (long) ( temp * idbase ); if ( temp < 0.0 ) carry -= 1.0; *a = temp - ( carry * dbase ); if ( --a < l ) break; } } #define AVOID_TRASHING 63 #define der limit-1 double mat1[limit+AVOID_TRASHING]; double mat2[limit+AVOID_TRASHING]; double mat3[limit+AVOID_TRASHING]; void LL(double di) { int lp; for(lp=0;lp al ) break; } di = di*di; idi = 1.0/di; for (;;){ du += 2.0; idu=1.0/du; if ( *dmat1 == 0.0 ){ dmat1++; dmat2++; if ( dmat1 > &mat1[der] ) break; } af=dmat1; bf=dmat2; tp2 = 0.0; tp1 = 0.0; for (;;){ tp1 = ( tp1 * dbase ) + *af; d1 = (double) (long) (tp1 * idi); *af = d1 ; tp1 -= ( d1 * di ); tp2 = ( tp2 * dbase ) + d1; d2 = (double) (long) (tp2 * idu); tp2 -= ( d2 * du ); *bf = *bf + d2;bf++;af++; if ( af > al ) break; } } adjust(&mat2[der],&mat2[0]); mul(&mat2[der],&mat2[0],2.0); } int main() { for(int lp=0;lp= digits_per_line ){ strncpy(out2,out1,digits_per_line); printf("%s\n",out2); strcpy(out1,&out1[digits_per_line]); } } printf("%s\n",out1); return(0); } _________________________________________________________________ _________________________________________________________________