diff options
Diffstat (limited to 'macros/plain/contrib/misc/pi.tex')
-rw-r--r-- | macros/plain/contrib/misc/pi.tex | 385 |
1 files changed, 385 insertions, 0 deletions
diff --git a/macros/plain/contrib/misc/pi.tex b/macros/plain/contrib/misc/pi.tex new file mode 100644 index 0000000000..ac7f060134 --- /dev/null +++ b/macros/plain/contrib/misc/pi.tex @@ -0,0 +1,385 @@ +% file pi.tex version 0.993 +% +% **** Compute Pi in TeX! **** +% +% +% Author: D. Roegel (roegel@loria.fr) +% +% Version 0.96: 22 July 1996 +% First release +% +% Version 0.97: 22 July 1996 +% Modified by J. Gelinas (jacquesg@clic.net.ca) +% Added one term to get correct last digits +% Added a second optional argument: +% 100 prints decimals 0..100 +% -2 100 prints only decimals 100..101 +% Print exact number of decimals +% +% Version 0.98: 23 July 1996 +% Modified by D. Roegel, so that the input -3 1 gives 141 +% and not 3.141. Also, some unnecessary braces inside \loop...\repeat +% were removed. \ShowResult shortened by introduction of \NextDigit. +% +% Version 0.99: 23 July 1996 +% Modified by D. Roegel. +% Improvements following suggestions by J. Gelinas (jacquesg@clic.net.ca). +% Bug corrected, which made some digits false (\fontdimen\firstpos\xb=0pt +% added in \updatefirstpos). +% \N replaced by \count2 in \ShowResult, this enabling calls +% such as \ShowResult{1/\the\N...}{...}. +% +% Version 0.991: 23 July 1996 +% Modified by D. Roegel. +% Two \UpdateFirstPos were removed. +% +% Version 0.992: 24 July 1996 +% Modified by D. Roegel, following several simplifications suggested +% by J. Gelinas (jacquesg@clic.net.ca). +% \Multiply removed, \Add and \Sub merged. +% \ComputeArcTan shortened. +% +% Version 0.993: 24 July 1996 +% Modified by D. Roegel, following several simplifications suggested +% by J. Gelinas (jacquesg@clic.net.ca). +% \ComputeArcTan1/n computes now completely arctan(1/n). +% +%------------------------------------------------------------------------ +% This programs uses the formula by John Machin: +% +% Pi=16*arctan(1/5)-4*arctan(1/239) +% +% For arctan(x), we use the development +% +% arctan(x)=\sum_{i=0}^\infty [{x^{2i+1}\over 2i+1} - {x^{2i+3}\over 2i+3}] +% +% One array (\xc) is used to store the partial sum up to {x^{2i+1}\over 2i+1}. +% A second array (\xa) us used to store the value of x^{2i+1}. +% A third array (\xb) is used to store {x^{2i+1}\over 2i+1}. +% +% The result is put in the array \xr. +% +% The implementation of arrays uses a trick shown by Tom Rokicki +% in his game of life program (life.tex). +% +% The number of digits you can compute depends on your implementation +% of TeX. I had no trouble computing 5000 digits. +% +% The last digit can be wrong. And if it is a 0 or a 9, this decimal +% and previous ones can be wrong too. However, the absolute error is +% no greater than one unit of the last digit. +% +% If you want to know more about Pi, check the file sci.math.faq. +% See also http://www.primus.com/staff/paulp/useless/pi.html +% and http://www.tu-chemnitz.de/~arndt/joerg.html. +% +%------------------------------------------------------------------------ +% +\newlinechar=`\^^J +\message{^^J***** Computation of Pi with John Machin's formula *****} +\message{^^J** i.e.: pi=16*arctan(1/5)-4*arctan(1/239)} +\message{^^JHow many decimals of pi do you want ? } + +\read16to\nbdigits +\newcount\n +\n\nbdigits +\ifnum\n<0 + \multiply\n-1 + \message{^^JFirst decimal to output ? } + \read16to\firstdigit +\else + \advance\n1 + \def\firstdigit{0} +\fi +\newcount\lastdigit +\lastdigit\firstdigit +\advance\lastdigit\n +\advance\lastdigit-1 + +\newcount\index +\index\lastdigit +\advance\index11 +\divide\index4 +\def\base{10000} +\def\basesp{10000sp} + +\newcount\lastplusone +\lastplusone\index +\advance\lastplusone1 + +% \index is now the index for the last slot in the arrays + +% slot 1 -> integer digits +% slot 2 -> digits 1 to 4 +% slot 3 -> digits 5 to 8 +% ... +% slot \index -> digits (\index-2) * 4 +1 to (\index-1) * 4 +% + +\font\xa=cmr10 at 11truept % array for current values of (1/5)^{2n+1} + % and (1/239)^{2n+1} +\fontdimen\lastplusone\xa=0sp % this creates room +\font\xb=cmr10 at 13truept % array for current values of (1/5)^{2n+1}/(2n+1) + % and (1/239)^{2n+1}/(2n+1) +\fontdimen\lastplusone\xb=0sp % this creates room +\font\xc=cmr10 at 15truept % array for current sums of arctan(1/5) + % and arctan(1/239) +\fontdimen\lastplusone\xc=0sp % this creates room +\font\xr=cmr10 at 17truept % array for the result +\fontdimen\lastplusone\xr=0sp % this creates room + +% (we have each time allocated one more slot than strictly necessary; +% this avoids a test on \lastplusone in \updatefirstpos) + +% \xa, \xb, \xc and \xr are now equal to 0 + +% Some variables (some of them are not strictly necessary, and might be +% replaced by \count's): + +\newcount\dv % will hold dividers +\newcount\firstpos % first non empty slot +\newcount\I % scratch register for loops +\newdimen\carry % for carry (in additions) and borrows (in subtractions) +\newdimen\x % a scratch variable +\newcount\N % counts the terms +\newif\ifcont % flag used to find when an operation on bignums is not done +\newcount\Sdv % value of one digit (used in \ShowResult) +\newcount\dir % toggle for alternating sums + +% Initialization of working arrays + +\def\InitializeArrays{ + { + \I=1 + \loop + \fontdimen\I\xa=0sp + \fontdimen\I\xb=0sp + \fontdimen\I\xc=0sp + \advance\I1 + \ifnum\I<\lastplusone + \repeat + } + } + +% Initialization of the result + +\newcount\I +\I=1 +\loop + \fontdimen\I\xr=0sp + \advance\I1 + \ifnum\I<\lastplusone +\repeat + + +% divide array #1 by #2 beginning at slot \firstpos and up to \index; +% result is in array #3 +% Maximum carry is 9999, so we need to be able to store 99999999. +% \dimen's can hold up to +/- 2,147,483,647 sp and +% \count's up to +/- 2,147,483,647, so it fits. + +\def\Divide#1#2into#3{% + \carry0sp + \I\firstpos + { + \loop + \x=\fontdimen\I#1 + \multiply\carry\base + \advance\x\carry + \carry\x + \divide\x#2 + \fontdimen\I#3=\x + \multiply\x#2 + \advance\carry-\x + \advance\I1 + \ifnum\I<\lastplusone + \repeat + } + } + +% Add or Subtract #2 to #3, depending on #1. array #2 is not modified. + +\def\Add#1#2to#3{% + \carry0sp + \I\index + { + \loop + \x=\fontdimen\I#3 + \advance\x by #1\fontdimen\I#2 + \advance\x by #1\carry + \fontdimen\I#3=\x + \carry\x + \divide\carry\base + \multiply\carry\base + \advance\x-\carry + \divide\carry\base + \ifdim\x<0sp + \advance\x\basesp + \advance\carry1sp + \fi + \fontdimen\I#3=\x + \advance\I-1 + \ifnum\I<\firstpos \ifnum\carry=0 \contfalse + \else \conttrue \fi + \else \conttrue + \fi + \ifcont + \repeat + } + } + + +% Multiply array #1 by #2. Result is in #1. +% This macro is not used in this program, and only remains +% here for didactic and historical reasons. + +\def\Multiply#1#2{% + \carry0sp + \I\index + \loop + \x=\fontdimen\I#1 + \multiply\x by #2 + \advance\x by \carry + \fontdimen\I#1=\x + \carry\x + \divide\carry\base + \multiply\carry\base + \advance\x-\carry + \fontdimen\I#1=\x + \divide\carry\base + \advance\I-1 + \ifnum\I<\firstpos \ifnum\carry=0 \contfalse + \else \conttrue \fi + \else \conttrue + \fi + \ifcont + \repeat + } + +% update value of \firstpos; in the worst case, \firstpos +% gets increased by 2. + +\def\updatefirstpos{ + \ifdim\fontdimen\firstpos\xa=0sp + \ifnum\firstpos<\lastplusone + \fontdimen\firstpos\xb=0sp + \advance\firstpos1 + \fi + \fi + } + +\def\UpdateFirstPos{\updatefirstpos\updatefirstpos} + +% Compute arctan(1/#1). Arrays \xa and \xb are used. Result is in \xc. + +\def\ComputeArcTan1/#1{% + \firstpos1 + \InitializeArrays + % initialize \xa with 1: + \fontdimen1\xa=1sp + \Divide\xa{#1}into\xa + % now, \xa contains 1/#1 + \Add1\xa to\xb + % \xb contains 1/#1 + \firstpos=2 + \dv=#1 + \multiply\dv\dv + \N1 + \dir-1 + \message{^^JI am now computing the following sum: arctan(1/#1)=1/#1} + \Add1\xb to\xc % first term + \loop + \Divide\xa\dv into\xa + \UpdateFirstPos + \ifnum\firstpos<\lastplusone + \advance\N2 + \Divide\xa\N into\xb + \ifnum\dir>0 + \message{+1/\the\N*1/#1^\the\N} + \else + \message{-1/\the\N*1/#1^\the\N} + \fi + \Add\dir\xb to\xc + \multiply\dir-1 + \repeat + } + +% Extract the next digit from \count1 + +\def\NextDigit#1{ + \ifnum\count2<\lastdigit + \Sdv\count1 + \divide\Sdv#1 + \advance\count2 1 + \ifnum\count2>\firstdigit + \edef\res{\res\number\Sdv} + \fi + \multiply\Sdv#1 + \advance\count1-\Sdv + \fi + } + +% Display the result; the digits are scanned one at a time, +% and those between \firstdigit and \lastdigit are extracted +% and saved in \res. + +\def\ShowResult#1#2{{ + \count0=2 + \count2=1 + \advance\lastdigit1 + \def\res{} + \loop + \count1=\fontdimen\count0#2 + \NextDigit{1000} + \NextDigit{100} + \NextDigit{10} + \NextDigit{1} + \advance\count0by 1 + \ifnum\count2<\lastdigit + \repeat + \advance\lastdigit-1 % for correct display in next line + \message{^^J#1\res} + } +} + +\ComputeArcTan1/5 + +%\ShowResult{arctan(1/5)=0.}\xc + +\message{^^JI multiply it by 4...} +\firstpos2 +\Add1\xc to\xc % alternative method is to write \Multiply\xc4 + % or to change \ComputeArcTan so that the initial value is 4 +\Add1\xc to\xc +\message{done} + +% add 4*arctan(1/5) to \xr +\Add1\xc to\xr + +%\ShowResult{4*arctan(1/5)=0.}\xr + +\ComputeArcTan1/{239} + +%\ShowResult{arctan(1/239)=0.}\xc + + +%\ShowResult{4*arctan(1/239)=0.}\xc + +\message{^^JI subtract arctan(1/239) to 4*arctan(1/5)...} + +\firstpos2 +\Add{-1}\xc to\xr +\message{done} + +\message{^^JAnd finally, I multiply it by 4 giving:} +\Add1\xr to\xr +\Add1\xr to\xr + +\ifnum\firstdigit=0 + \ShowResult{pi[0 .. \number\lastdigit]=3.}\xr +\else + \ShowResult{pi[\firstdigit.. \number\lastdigit]=}\xr +\fi + +\bye + |