This page presents tiny C programs written by the author of the present page. These programs compute some familiar mathematical constants to determine their decimal expression in a few thousands digit accuracy. The source codes given in this page are all written for 32 bit ANSI C compilers. In particular, int variables are assumed to have 32 bit width.
A simple way of changing accuracy is shown for each program, so that you can easily modify the code to increase its accuracy. It should be noted however that then the time required for the computation also increases rapidly (approximately in propotional to the square of the accuracy measured in the number of digits to be computed) due to the algorithm adopted in these programs.
Further programs and related materials can be found at web sites listed in Related Links.
int a[9985],b=10000,d=1,i,n,x;main(){for(n=b;i?x=x/d*i,d=--i?
i-~i:b:(i=n-=4*printf("%.4d",x/d));a[i]=x%d)x+=*a?a[i]*b:2;}
This program computes pi to the 2488th decimal place. The source code is
very short, only of 121 characters long.
The program uses spigot algorithm applied to a series expansion formula
for arctan x found by Euler.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||||
|---|---|---|---|---|---|
| 2488 (original) | a[9985] | b=10000 | n=b | :b: | %.4d |
| 6800 | a[16*1703-15] | n=16*1703 | |||
| 30000 | a[12*10004-11] | b=1000 | n=12*10004 | :b*b: | %.3d |
| 120000 | a[12*40004-11] | n=12*40004 | |||
int*a,*q,w[73*33+4*303],b=33,n=303,d,e,h=1,j,k,m=10000,p,r,t,x,
y,f[]={23,21,15,14,251,449,-4801,8749,241,97,67,107};main(){for
(q=w;p=f[r]*b,k=r++<4;q+=p)for(;k<p;k+=3)q[k]=2;for(;9/h++;)for
(x=y=j=0;++j<n;y=y%m*m){if(h<3)for(r=0,a=w;p=f[r]*b,r<4;q[4*j+r
++]=d,a+=p)for(d=0,k=p;k-=t=1;a[k]=d%(e=k%3?f[r+4]:(t=k/3*2-1)+
2),d=d/e*t)d+=a[k]*m;for(r=4;r--;x+=(72499379*(h==8|!(h&3)+!(h&
1))+91945470*(!(h-9)+!(h%3))+35752517*!(h%5)+165880261*!(h-7))%
f[r+8]*q[4*j+r]);y+=x/m;printf(j<2?"log %d":j<3?" = %d.":j+1-n?
"%.4d":"%.4d\n",j<2?h:y/m);x=x%m*m;}}
This program of 541 characters long computes log 2, log 3, ... , log 10
simultaneously to the 1200th decimal place.
The program uses spigot algorithm applied to the Taylor series of arctanh x
with x = 1/251, 1/449, 1/4801, 1/8749.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||
|---|---|---|---|
| 1200 (original) | w[73*33+4*303] | b=33 | n=303 |
| 12000 | w[73*327+4*3003] | b=327 | n=3003 |
| 120000 | w[73*3263+4*30003] | b=3263 | n=30003 |
int*p,w[2048],n=825,x;main(){for(p=w;++p-w>>11?p=
w,x=n--:1;)n?x=x%n*10+*p,*p=x/n:printf("%d",*p);}
This program computes e to the 2048th decimal place (the integer part of e is not shown).
The source code is very short, only of 98 characters long.
The algorithm used is a traditional one, closely related to the Horner representation
of a quite well-known series expansion (spigot algorithm is not used).
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||
|---|---|---|---|
| 2048 (original) | w[2048] | n=825 | >>11? |
| 4096 | w[4096] | n=1495 | >>12? |
| 8192 | w[8192] | n=2729 | >>13? |
| 16384 | w[16384] | n=5017 | >>14? |
| 131072 | w[1<<17] | n=32179 | >>17? |
| 1048576 | w[1<<20] | n=214152 | >>20? |
The first 3 cases in the table, i.e., the cases in which n is initialized below 215/10, have compatibility also with 16 bit C compilers.
int*d,*e,w[2256],b=10000,i,j,k,p,r=2,x;main(){while(j?k
&1?x=x%i*b+*--e,*--d=x/i:(*d=(x+=k-2?b-*d-1+*e++:(~r?r?
2916:3456:3072)**d)%b,d++,x/=b),j--:k?d=w+k/2*3*p,j=k-k
/4,x=k--&1,e=w+(x?!k:6)*p,j*=p:++i<33*p?k=4:~r?p=376,i=
0,w[5*p]=r--?r?23:58:36:printf("%.4d",w[--p])&&p>60);}
This program of 274 characters long computes Euler's constant to the 1264th
decimal place. The program uses a series expansion of the exponential
integral Ei(-x) evaluated at three points (x = 2916, 3456 and 3072 in the
original code). The evaluation at three point is to eliminate logarithm term
from the series expansion in a certain linear combination. The series is
numerically summed with a conventinal method (spigot algorithm is not used).
The author would like to mention that the task of executing this program
is somewhat heavy; it takes for 10 to 30 seconds with a 1GHz class processor,
and no messages are shown unitl the computation is finished.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||||
|---|---|---|---|---|---|
| 1264 (original) | w[2256] | 2916:3456:3072 | p=376 | 23:58:36 | p>60 |
| 3800 | w[6756] | 8748:10368:9216 | p=1126 | 26:66:41 | p>176 |
| 7116 | w[12012] | 16384:18432:17496 | p=2002 | 71:98:28 | p>223 |
| 10672 | w[18012] | 24576:27648:26244 | p=3002 | 74:102:29 | p>334 |
| 21348 | w[36024] | 49152:55296:52488 | p=6004 | 79:109:31 | p>667 |
int a[4*8819],b=8819,n=2*667+1,d,e,g,j,k,m=1000,r,t,u=7000;main
(){for(;k<3*b;k+=2)t=k+2&6?k&8?-u:u:0,a[k+b]=k&2?-6*t:6*t,k<b?a
[k]=t:0;for(u*=40,d=a[b]-*a;--n;d=r?printf("%.3d",(d-=e)/u),d%u
:d*m+e)for(e=r=n&1,k=r?b:3*b;t=k--/2;e=e/g*t)g=r^!(k-r&3)?2:1,g
*=t+1,e+=a[j=k+!r*b]*m,a[j]=e%g;}
This program of 285 characters long computes Catalan's constant
to the 2000th decimal place.
The program uses spigot algorithm applied to Broadhurst's formula.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||
|---|---|---|---|
| 2000 (original) | a[4*8819] | b=8819 | n=2*667+1 |
| 5087 | a[4*22500] | b=22500 | n=2*1696+1 |
| 10166 | a[4*45000] | b=45000 | n=2*3389+1 |
| 33875 | a[4*150000] | b=150000 | n=2*11292+1 |
int a[9922],b=9922,n=667,d,e,k,m=1000,t,u=7000;main(){for(;k<b;
k+=3)a[k]=k%6?-u:u;for(u*=8;n--;d+=*a*m,printf("%.3d",d/u),*a=d
%u)for(d=0,k=b;t=(--k+2)/3;d=d/e*t)d+=a[k]*m,a[k]=d%(e=k%3?t+1:
t*4+2);}
This program of 197 characters long computes zeta(3)
to the 2000th decimal place.
The program uses spigot algorithm applied to Apery's formula.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||
|---|---|---|---|
| 2000 (original) | a[9922] | b=9922 | n=667 |
| 10043 | a[50000] | b=50000 | n=3348 |
| 50183 | a[250000] | b=250000 | n=16728 |
| 150527 | a[750000] | b=750000 | n=50176 |
int*p,a[23*4392],b=15*4392,n=667,d,e,f,g,h,k,m=1000,q=313255,r,
t,c[]={31,173,1,-1614,284,0,-31,-173,-1,-6212,-457,-1,-31,-173,
-1,-1614,284,0,31,173,1,74552,-111,1};main(){for(;r+3||(d+=f/q*
64,f=f%q*m,e=m*m,h++<3||printf("%.3d",d/e)+n--,d=d%e*m,r=0,n);r
--)for(p=r?p+=~r?b/3:b:a,f+=(r?~r?-369:14:72)**(h?p:c-r),e=0,k=
b/(1-2*r);t=(--k+4)/5;p[k]=e%g,*p=e/g*t)g=k%5<1-r-t%2?2*t+2:t+1
,e=*p+(h?p[k]:k%5?0:c[k/5%8*3-r])*m;}
This program of 415 characters long computes zeta(5)
to the 2000th decimal place.
The program uses spigot algorithm applied to Broadhurst's formula.
The accuracy of computation can be modified as in the following table.
| accuracy [digits] | replacements in the code | ||
|---|---|---|---|
| 2000 (original) | a[23*4392] | b=15*4392 | n=667 |
| 5210 | a[23*11500] | b=15*11500 | n=1737 |
| 15824 | a[23*35000] | b=15*35000 | n=5275 |
The web site given by X. Gourdon and P. Sebah. This site presents a lot of contents about mathematical constants. In the Tiny programs for constants computation page you can find C programs of the same kind as those in the present page.
Jasonp's web site. In the Pi Program Pile page you can find a lot of programs for computing pi.