summaryrefslogtreecommitdiff
path: root/macros/plain/contrib/misc/pi.tex
diff options
context:
space:
mode:
Diffstat (limited to 'macros/plain/contrib/misc/pi.tex')
-rw-r--r--macros/plain/contrib/misc/pi.tex385
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
+