| T.R | Title | User | Personal
Name | Date | Lines |
|---|
| 92.1 | | HARE::STAN | | Thu Jul 05 1984 21:40 | 3 |
| +0unqn"e20un'buhbuvhkqn<jbuqqn*10/3uiqi<\+2*10+(qq*qi)uablkqi*2-1ujqa/qjuq
qa-(qq*qj)-2\10i-1%i>qq/10utqh+qt+48uwqw-58"e48uw%v'qv"nqv^t'qwuvqq-(qt*10)uh>
|
| 92.2 | | HARE::STAN | | Thu Jul 05 1984 21:44 | 27 |
| IMPLICIT INTEGER (A-Z)
C
C PROGRAM TO CALCULATE PI
C ORIGINAL ALGORITHM BY STANLEY RABINOWITZ
C Please do not publish this without my permission. - STAN -
C
DIMENSION VECT(2000)
DATA VECT/2000*2/
DATA HOLD/0/
DATA ZERO/'0'/
NUM=500
LEN=10*NUM/3
DO 2 N=1,NUM
Q=0
DO 3 L=LEN,1,-1
X=10*VECT(L)+Q*L
Q=X/(2*L-1)
VECT(L)=X-Q*(2*L-1)
3 CONTINUE
X=Q/10
R=Q-X*10
TYPE 100, HOLD+X+ZERO
100 FORMAT('+',A1$)
HOLD=R
2 CONTINUE
STOP
END
|
| 92.3 | | HARE::STAN | | Fri Jul 06 1984 14:12 | 3 |
| A good book to look at is A History of Pi by Petr Beckman,
St. Martin's Press: 1971. I think it's available in Dalton's
in the Mall of New Hampshire.
|
| 92.4 | | CASTOR::[RABAHY] | | Fri Jul 06 1984 15:10 | 6 |
| RE: .1
What is that, the object form of reply 2? I don't think it's APL although
it has been a while.
David.
|
| 92.5 | | HARE::STAN | | Fri Jul 06 1984 20:53 | 1 |
| It's a source file. Not APL though. It would be much shorter in APL.
|
| 92.6 | | AURORA::HALLYB | | Sat Jul 07 1984 17:03 | 8 |
| In APL, of course, it would only be one line long.
Back in the early 60s didn't the RAND corporation compute PI to one million
places, and use that output for their hit book "1,000,000 random digits"?
As I recall, Scientific American quoted Dr. Matrix as saying that, when
properly interpreted, PI contains the entire history of the human race.
Interesting that the quote is true not just for pi but any irrational...
|
| 92.7 | | RANI::LEICHTERJ | | Sat Jul 07 1984 21:49 | 17 |
| The RAND book of random numbers was NOT based on PI. Those digits were produced
using some sort of electronic random signal source - most likely a noisy
diode, I'd guess.
No one knows how "random" the digits of PI are. There are all sorts of
results that will tell you that all but "very few" (in some sense, usually
either a countable set or a set of measure 0) numbers must have random
expansions, for any reasonable definition of random - but that proves
nothing about any PARTICULAR number. (In fact, if you use the Kolmogorov-
Chaitin definition of randomness, in which a number is "random" if there
is no program much shorter than the number itself that prints the number
out, then PI is certainly NOT random, since there are very short programs
(e.g. reply .1) that print as many digits of PI as you like.)
The same argument, BTW applies to whether numbers "contain all the past
history of the human race".
-- Jerry
|
| 92.8 | | XENON::GAUDREAU | | Wed Jul 11 1984 13:03 | 145 |
| For Your Amusement or Edification - PI to 10,000 places using Stan's algorithm.
When I get Petr Beckmann's book from the company library, I'll check each
and every digit but I don't think I'll find anything wrong here...
Joe
-=-
3.141592653589793238462643383279502884197169399375105820974944592307816406
286208998628034825342117067982148086513282306647093844609550582231725359408
128481117450284102701938521105559644622948954930381964428810975665933446128
475648233786783165271201909145648566923460348610454326648213393607260249141
273724587006606315588174881520920962829254091715364367892590360011330530548
820466521384146951941511609433057270365759591953092186117381932611793105118
548074462379962749567351885752724891227938183011949129833673362440656643086
021394946395224737190702179860943702770539217176293176752384674818467669405
132000568127145263560827785771342757789609173637178721468440901224953430146
549585371050792279689258923542019956112129021960864034418159813629774771309
960518707211349999998372978049951059731732816096318595024459455346908302642
522308253344685035261931188171010003137838752886587533208381420617177669147
303598253490428755468731159562863882353787593751957781857780532171226806613
001927876611195909216420198938095257201065485863278865936153381827968230301
952035301852968995773622599413891249721775283479131515574857242454150695950
829533116861727855889075098381754637464939319255060400927701671139009848824
012858361603563707660104710181942955596198946767837449448255379774726847104
047534646208046684259069491293313677028989152104752162056966024058038150193
511253382430035587640247496473263914199272604269922796782354781636009341721
641219924586315030286182974555706749838505494588586926995690927210797509302
955321165344987202755960236480665499119881834797753566369807426542527862551
818417574672890977772793800081647060016145249192173217214772350141441973568
548161361157352552133475741849468438523323907394143334547762416862518983569
485562099219222184272550254256887671790494601653466804988627232791786085784
383827967976681454100953883786360950680064225125205117392984896084128488626
945604241965285022210661186306744278622039194945047123713786960956364371917
287467764657573962413890865832645995813390478027590099465764078951269468398
352595709825822620522489407726719478268482601476990902640136394437455305068
203496252451749399651431429809190659250937221696461515709858387410597885959
772975498930161753928468138268683868942774155991855925245953959431049972524
680845987273644695848653836736222626099124608051243884390451244136549762780
797715691435997700129616089441694868555848406353422072225828488648158456028
506016842739452267467678895252138522549954666727823986456596116354886230577
456498035593634568174324112515076069479451096596094025228879710893145669136
867228748940560101503308617928680920874760917824938589009714909675985261365
549781893129784821682998948722658804857564014270477555132379641451523746234
364542858444795265867821051141354735739523113427166102135969536231442952484
937187110145765403590279934403742007310578539062198387447808478489683321445
713868751943506430218453191048481005370614680674919278191197939952061419663
428754440643745123718192179998391015919561814675142691239748940907186494231
961567945208095146550225231603881930142093762137855956638937787083039069792
077346722182562599661501421503068038447734549202605414665925201497442850732
518666002132434088190710486331734649651453905796268561005508106658796998163
574736384052571459102897064140110971206280439039759515677157700420337869936
007230558763176359421873125147120532928191826186125867321579198414848829164
470609575270695722091756711672291098169091528017350671274858322287183520935
396572512108357915136988209144421006751033467110314126711136990865851639831
501970165151168517143765761835155650884909989859982387345528331635507647918
535893226185489632132933089857064204675259070915481416549859461637180270981
994309924488957571282890592323326097299712084433573265489382391193259746366
730583604142813883032038249037589852437441702913276561809377344403070746921
120191302033038019762110110044929321516084244485963766983895228684783123552
658213144957685726243344189303968642624341077322697802807318915441101044682
325271620105265227211166039666557309254711055785376346682065310989652691862
056476931257058635662018558100729360659876486117910453348850346113657686753
249441668039626579787718556084552965412665408530614344431858676975145661406
800700237877659134401712749470420562230538994561314071127000407854733269939
081454664645880797270826683063432858785698305235808933065757406795457163775
254202114955761581400250126228594130216471550979259230990796547376125517656
751357517829666454779174501129961489030463994713296210734043751895735961458
901938971311179042978285647503203198691514028708085990480109412147221317947
647772622414254854540332157185306142288137585043063321751829798662237172159
160771669254748738986654949450114654062843366393790039769265672146385306736
096571209180763832716641627488880078692560290228472104031721186082041900042
296617119637792133757511495950156604963186294726547364252308177036751590673
502350728354056704038674351362222477158915049530984448933309634087807693259
939780541934144737744184263129860809988868741326047215695162396586457302163
159819319516735381297416772947867242292465436680098067692823828068996400482
435403701416314965897940924323789690706977942236250822168895738379862300159
377647165122893578601588161755782973523344604281512627203734314653197777416
031990665541876397929334419521541341899485444734567383162499341913181480927
777103863877343177207545654532207770921201905166096280490926360197598828161
332316663652861932668633606273567630354477628035045077723554710585954870279
081435624014517180624643626794561275318134078330336254232783944975382437205
835311477119926063813346776879695970309833913077109870408591337464144282277
263465947047458784778720192771528073176790770715721344473060570073349243693
113835049316312840425121925651798069411352801314701304781643788518529092854
520116583934196562134914341595625865865570552690496520985803385072242648293
972858478316305777756068887644624824685792603953527734803048029005876075825
104747091643961362676044925627420420832085661190625454337213153595845068772
460290161876679524061634252257719542916299193064553779914037340432875262888
963995879475729174642635745525407909145135711136941091193932519107602082520
261879853188770584297259167781314969900901921169717372784768472686084900337
702424291651300500516832336435038951702989392233451722013812806965011784408
745196012122859937162313017114448464090389064495444006198690754851602632750
529834918740786680881833851022833450850486082503930213321971551843063545500
766828294930413776552793975175461395398468339363830474611996653858153842056
853386218672523340283087112328278921250771262946322956398989893582116745627
010218356462201349671518819097303811980049734072396103685406643193950979019
069963955245300545058068550195673022921913933918568034490398205955100226353
536192041994745538593810234395544959778377902374216172711172364343543947822
181852862408514006660443325888569867054315470696574745855033232334210730154
594051655379068662733379958511562578432298827372319898757141595781119635833
005940873068121602876496286744604774649159950549737425626901049037781986835
938146574126804925648798556145372347867330390468838343634655379498641927056
387293174872332083760112302991136793862708943879936201629515413371424892830
722012690147546684765357616477379467520049075715552781965362132392640616013
635815590742202020318727760527721900556148425551879253034351398442532234157
623361064250639049750086562710953591946589751413103482276930624743536325691
607815478181152843667957061108615331504452127473924544945423682886061340841
486377670096120715124914043027253860764823634143346235189757664521641376796
903149501910857598442391986291642193994907236234646844117394032659184044378
051333894525742399508296591228508555821572503107125701266830240292952522011
872676756220415420516184163484756516999811614101002996078386909291603028840
026910414079288621507842451670908700069928212066041837180653556725253256753
286129104248776182582976515795984703562226293486003415872298053498965022629
174878820273420922224533985626476691490556284250391275771028402799806636582
548892648802545661017296702664076559042909945681506526530537182941270336931
378517860904070866711496558343434769338578171138645587367812301458768712660
348913909562009939361031029161615288138437909904231747336394804575931493140
529763475748119356709110137751721008031559024853090669203767192203322909433
467685142214477379393751703443661991040337511173547191855046449026365512816
228824462575916333039107225383742182140883508657391771509682887478265699599
574490661758344137522397096834080053559849175417381883999446974867626551658
276584835884531427756879002909517028352971634456212964043523117600665101241
200659755851276178583829204197484423608007193045761893234922927965019875187
212726750798125547095890455635792122103334669749923563025494780249011419521
238281530911407907386025152274299581807247162591668545133312394804947079119
153267343028244186041426363954800044800267049624820179289647669758318327131
425170296923488962766844032326092752496035799646925650493681836090032380929
345958897069536534940603402166544375589004563288225054525564056448246515187
547119621844396582533754388569094113031509526179378002974120766514793942590
298969594699556576121865619673378623625612521632086286922210327488921865436
480229678070576561514463204692790682120738837781423356282360896320806822246
801224826117718589638140918390367367222088832151375560037279839400415297002
878307667094447456013455641725437090697939612257142989467154357846878861444
581231459357198492252847160504922124247014121478057345510500801908699603302
763478708108175450119307141223390866393833952942578690507643100638351983438
934159613185434754649556978103829309716465143840700707360411237359984345225
161050702705623526601276484830840761183013052793205427462865403603674532865
105706587488225698157936789766974220575059683440869735020141020672358502007
245225632651341055924019027421624843914035998953539459094407046912091409387
001264560016237428802109276457931065792295524988727584610126483699989225695
968815920560010165525637567
-=*=- EOF -=*=-
|
| 92.9 | | HARE::STAN | | Wed Jul 11 1984 17:12 | 3 |
| This is correct to at least 9900 places. Unfortunately,
my edition of Petr Beckmann's book has the last line of
pi (places 9901 to 10000) cut off at the bottom.
|
| 92.10 | | HARE::STAN | | Wed Jul 11 1984 17:12 | 1 |
| Now what's the longest repeated string?
|
| 92.11 | | HARE::STAN | | Wed Jul 11 1984 19:26 | 13 |
| Here are some references to pi algorithms:
R. P. Brent, Fast Multiple-Precision Evaluation of Elementary Functions,
J. ACM, 23(1976)242-251.
E. Salamin, Computation of Pi using Arithmetic-Geometric Mean,
Math. Comput. 135(1976)565-570.
D. Shanks and J. W. Wrench, Jr. Calculation of pi to 100,000 decimals.
Math Comput. 16(1962)76-79.
J. W. Wrench, Jr. The Evolution of Extended Decimal Approximations to Pi.
The mathematics teacher, 53(1960)644-650.
|
| 92.12 | | HARE::STAN | | Wed Jul 11 1984 19:42 | 30 |
| The current (July 1984) issue of Siam Review has a good article
on calculating Pi:
J. M. Borwein and P. B. Borwein, The Arithmetic-Geometric Mean and Fast
Computation of Elementary Functions, Siam Review, 26(1984)351-366.
Here is their algorithm: (like most algorithms for computing pi,
it needs a multiprecision package):
Define a recursion as follows:
a = sqrt(2) b = 0 p = 2 + sqrt(2).
0 0 0
a = ( sqrt(a + 1/sqrt(a ) ) / 2
n+1 n n
b = sqrt(a ) ( (b + 1)/(b + a ) )
n+1 n n n n
p = p b ( (1 + a )/(1 + b ) ) .
n+1 n n+1 n+1 n+1
Then p converges exponentially to pi and | p - pi | < 1/10**(2**n).
n n
Twenty iterations of this algorithm will give pi to 2 million places.
Each iteration requires about ten multiprecise operations. The algorithm
is very stable with all the operations being performed on numbers
between 1/2 and 7.
|
| 92.13 | | HARE::STAN | | Wed Jul 11 1984 19:44 | 5 |
| Traditional pi algorithms require a multiprecision package
and usually calculate pi from Taylor series for arctangents.
A typical formula is
pi = 24 arctan(1/8) + 8 arctan(1/57) + 4 arctan(1/239) .
|
| 92.14 | | HARE::STAN | | Wed Jul 11 1984 21:32 | 22 |
| C Traditional calculation of pi using the formula
C
C PI/4=4*ARCTAN(1/5)-ARCTAN(1/239)
C
C and using the Brent Multiprecision Package (MP).
C
COMMON B, T, M, LUN, MXR, R
INTEGER B, T, R(258)
INTEGER PI(64), C(110), TEMP(64)
CALL MPSET (6, 110, 64, 258)
CALL MPART1(5,TEMP)
CALL MPMULI(TEMP,4,TEMP)
CALL MPART1(239,PI)
CALL MPSUB(TEMP,PI,PI)
CALL MPMULI(PI,4,PI)
C
C CONVERT TO PRINTABLE FORMAT (F102.100) AND WRITE OUT
C
CALL MPOUT (PI,C,102,100)
WRITE (LUN,10) C
10 FORMAT(1X,60A1)
END
|
| 92.15 | | TURTLE::GILBERT | | Wed Jul 11 1984 21:52 | 6 |
| The longest repeated sequence of digits in the 10,000-place value of pi
(problem and value from previous responses) is "7111369".
- Gilbert
Yes, I have a (poor) computer algorithm for this.
|
| 92.16 | | HARE::STAN | | Thu Jul 12 1984 02:59 | 3 |
| I calculated pi to 10,000 places using the Borwein algorithm
(program is in next response) and got the same answer that my
program gave.
|
| 92.17 | | HARE::STAN | | Thu Jul 12 1984 03:00 | 61 |
| C Calculation of pi using the Arithemetic-Geometric Mean
C algorithm of Borwein and Borwein.
C
C This uses the Brent Multiprecision Package (MP).
C
COMMON BB, T, M, LUN, MXR, R
PARAMETER (PREC=10000)
INTEGER BB, T, R(PREC*4+20)
INTEGER PI(PREC), A(PREC), B(PREC)
INTEGER NEXTA(PREC), NEXTB(PREC)
INTEGER TEMP(PREC)
INTEGER TEMP2(PREC)
INTEGER SUM(PREC)
INTEGER SQRTA(PREC)
INTEGER C(10010)
CALL MPSET (6, 10010, PREC, PREC*4+20)
C
C Set A[0]=SQRT(2), B[0]=0, PI[0]=2+SQRT(2)
C
CALL MPCIM(2,A)
CALL MPSQRT(A,A)
CALL MPCIM(0,B)
CALL MPCIM(2,PI)
CALL MPADD(A,PI,PI)
DO 50 N=1,14
TYPE *, 'N=', N
C
C Compute A[n+1], B[n+1], Pi[n+1] from A[n], B[n], Pi[n].
C
C
C A[n+1]=( sqrt(A[n])+1/sqrt(A[n]) ) / 2
C
CALL MPSQRT(A,SQRTA)
CALL MPREC(SQRTA,TEMP)
CALL MPADD(SQRTA,TEMP,NEXTA)
CALL MPDIVI(NEXTA,2,NEXTA)
C
C B[n+1]=sqrt(A[n]) * (B[n]+1)/(B[n]+A[n])
C
CALL MPADDI(B,1,TEMP)
CALL MPADD(A,B,SUM)
CALL MPDIV(TEMP,SUM,TEMP)
CALL MPMUL(TEMP,SQRTA,NEXTB)
C
C Pi[n+1]=B[n+1] * Pi[n] * (A[n+1]+1)/(B[n+1]+1)
C
CALL MPADDI(NEXTA,1,TEMP)
CALL MPADDI(NEXTB,1,TEMP2)
CALL MPDIV(TEMP,TEMP2,TEMP)
CALL MPMUL(TEMP,NEXTB,TEMP)
CALL MPMUL(TEMP,PI,PI)
CALL MPADDI(NEXTA,0,A)
CALL MPADDI(NEXTB,0,B)
50 CONTINUE
C
C CONVERT TO PRINTABLE FORMAT (F10002.10000) AND WRITE OUT
C
CALL MPOUT (PI,C,10002,10000)
WRITE (LUN,10) C
10 FORMAT(1X,80A1)
END
|
| 92.18 | | XENON::GAUDREAU | | Thu Jul 12 1984 17:04 | 246 |
| Here is the Pascal version of Stan's Algorithm -
Joe
-=-
Program Pi_Maker (Input, Output);
{==============================================================================
PI Maker - Joseph Gaudreau on 10-July-1984
This program uses an algorithm created by Stanley Rabinowitz of the
Spitbrook DEC (ZK02-3, Hare::Stan) to compute the irrational number PI
to any desired accuracy. The technique is not yet commented on and Stan
has asked not to have his algorithm published as he has not yet done so
himself.
This is the original Fortran program that outputs PI -
IMPLICIT INTEGER (A-Z)
C
C PROGRAM TO CALCULATE PI
C ORIGINAL ALGORITHM BY STANLEY RABINOWITZ
C Please do not publish this without my permission. - STAN -
C
DIMENSION VECT(2000)
DATA VECT/2000*2/
DATA HOLD/0/
DATA ZERO/'0'/
NUM=500
LEN=10*NUM/3
DO 2 N=1,NUM
Q=0
DO 3 L=LEN,1,-1
X=10*VECT(L)+Q*L
Q=X/(2*L-1)
VECT(L)=X-Q*(2*L-1)
3 CONTINUE
X=Q/10
R=Q-X*10
TYPE 100, HOLD+X+ZERO
100 FORMAT('+',A1$)
HOLD=R
2 CONTINUE
STOP
END
Sample output of the Pascal program -
3.141592653589793238462643383279502884197169399375105820974944592307816406
286208998628034825342117067982148086513282306647093844609550582231725359408
128481117450284102701938521105559644622948954930381964428810975665933446128
475648233786783165271201909145648566923460348610454326648213393607260249141
273724587006606315588174881520920962829254091715364367892590360011330530548
820466521384146951941511609433057270365759591953092186117381932611793105118
5480744623799627495673518857527248912279381830119491
==============================================================================}
Const
Line_Width = 75; { Output line length }
Num_Digits = 500; { Number of digits to compute }
Max_Calc = 2000; { Four times Num_Digits }
Type
Vector = Array[1 .. Max_Calc] of Integer; { Work space }
Store = Array[1 .. Num_Digits] of Integer; { Holds number }
Var
Take : Store; { Holds the number }
{=============================================================================}
{ Init workspace with two's }
{=============================================================================}
Procedure Zero ( Var Vect : Vector; { Works space of 4*Num_Digits }
Maxl : Integer ); { Should equal Max_Calc }
Var
Loop : Integer; { Loop variable }
Begin
For Loop := 1 to Maxl do { Loop through all elements and }
Vect[Loop] := 2 { init it }
{End FOR}
End; {Zero}
{=============================================================================}
{ Compute PI }
{=============================================================================}
Procedure PI_Calc ( Var Take : Store; { Base 10 of PI }
Maxl : Integer ); { Compute # }
Var
L,
N,
Q,
R,
X,
Len,
Pos, { Position within Take array }
Hold : Integer;
Vect : Vector;
Begin
Zero (Vect, Max_Calc); { Clear out array }
Pos := 1; { 1st element of storage }
Hold := 0;
Len := 10 * Maxl div 3;
For N := 1 to Maxl do
Begin
Q := 0;
For L := Len downto 1 do
Begin
X := 10 * Vect[L] + Q * L;
Q := X div (2 * L - 1);
Vect[L] := X - Q * (2 * L - 1)
End;
{End FOR}
X := Q div 10;
R := Q mod 10;
Take[Pos] := Hold + X; { Store a single digit 0-10 }
Pos := Pos + 1; { Next storage location }
Hold := R
End
{End FOR}
End; {Pi_Calc}
{=============================================================================}
{ Convert digit 10 }
{=============================================================================}
Procedure Fix_It (Var Take : Store; { Holds "PI" }
Pos : Integer ); { Current position }
Begin
Take[Pos] := 0; { Current - 10 = 0 }
Take[Pos - 1] := Take[Pos - 1] + 1; { Carry the one }
If Take[Pos - 1] = 10 then { Check for possible carry }
Fix_It (Take, Pos - 1) { if yes, recurse & fix }
{End IF}
End; {Fix_IT}
{=============================================================================}
{ Fix the number PI to true base 10 }
{=============================================================================}
Procedure Fix ( Var Take : Store; { Holds "PI" }
Maxl : Integer ); { Maximum position }
Var
Loop : Integer; { Loop variable }
Begin
For Loop := 1 to Maxl do { Loop through all digits }
If Take[Loop] = 10 then { If ten, must convert }
Fix_It (Take, Loop) { do it }
{End IF}
{End FOR}
End; {Fix}
{=============================================================================}
{ Print out PI }
{=============================================================================}
Procedure Print_IT ( Take : Store; { Holds "PI" }
Maxl, { Maximum array size }
Width : Integer ); { Width of output }
Var
Len, { Current line length }
Loop : Integer; { Loop variable }
Begin
Len := 3; { Length of a space, three and decimal pnt. }
Write (Output, ' 3.'); { Used only for formatting convenience,
does not represent "cheating" }
For Loop := 2 to Maxl do { Loop through remaining #'s }
Begin
Write (Output, Take[Loop]:1 ); { Write a single number }
Len := Len + 1; { Increment length }
If Len >= Width then { Test for newline }
Begin
Len := 0; { Zero line length }
Writeln (Output) { New line }
End
{End IF}
End;
{End FOR}
Writeln (Output); { Blank lines at end }
Writeln (Output)
End; {Print_It}
{=============================================================================}
{ Main program }
{=============================================================================}
Begin
Pi_Calc (Take, Num_Digits); { Compute PI }
Fix (Take, Num_Digits); { Convert PI to normal base 10 }
Print_It (Take, Num_Digits, Line_Width) { Print PI }
End. {Main}
|
| 92.19 | | HARE::STAN | | Thu Jul 12 1984 19:55 | 4 |
| Note that the algorithm I gave in 92.1 and 92.2 is intended to be
a real small piece of code for computing pi. I have faster methods
for computing pi using a similar algorithm, but it takes more code.
I may include it here when I get back from vacation.
|
| 92.20 | | HARE::STAN | | Thu Jul 12 1984 19:55 | 39 |
| Here are some more references to calculating pi:
Ballentine, J. P. The best (?) formula for computing pi to a thousand
places. American Mathematical Monthly. 46(1939)499-501.
Camp, C. C. A new calculation of pi. American Mathematical Monthly.
33(1926)472-473.
Frame, J. S. A series useful in the computation of pi. American Mathematical
Monthly. 42(1935)499-501.
Guinard, A. P. An asymptotic series for computing pi. Math. Gazette.
29(1945)214-218.
Reitwiesner, G. W. An ENIAC determination of pi and e to more than 2000
decimal places. Mathematical Tables and Other Aids to Computation.
4(1950)11-15.
Schepler, Herman. The chronology of pi. Mathematics Magazine. 23(1950)
165-170,216-228,279-283.
Eves, Howard. The latest about pi. Math Teacher. 55(1962)129-130.
Genuys, F. Dix milles decimals de pi. Chiffres. 1(1958)17-22.
Nicholson, S. C. and Jeneel, J. Some comments on a NORC computation of pi.
Mathematical Tabels and other Aids to Computation. 9(1955)162-164.
Wrench, J. W. Jr. and Eves, Howard, eds. The evolution of the extended
decimal approximation to pi. Math Teacher. 53(1960)644-650.
Knuth, Donald E. Art of Computer Programming, Vol 2 (1971) page 248.
Maier, Bruce. Comparison between various methods for determining pi with
respect to running time in a computer. School Science and
Mathematics. 72(1972)777-781.
Shanks, Daniel and J. W. Wrench Jr. Calculation of pi to 100,000 decimals.
Mathematics of Computation. 16(1962)76-99.
|
| 92.21 | | HARE::STAN | | Fri Jul 13 1984 04:57 | 8 |
| In the next response, I give a revised program that runs 4 times faster
than the original program. I used it to calculate pi to 10,000 places
(by changing the NUM=504 to NUM=10004) and timed it against the Borwein
algorithm. The Borwein algorithm took 2:26.30, my algorithm took 0:58:47.
It's much faster than traditional algorithms for such small number of
digits (10,000). Unfortunately, it won't beat any records (as in millions
of digits) because it gets much slower for that many digits and also
it page faults to death.
|
| 92.22 | | HARE::STAN | | Fri Jul 13 1984 04:57 | 23 |
| IMPLICIT INTEGER (A-Z)
C
C PROGRAM TO CALCULATE PI
C ORIGINAL ALGORITHM BY Stanley Rabinowitz
C
PARAMETER (NUM=504,LEN=10*NUM/3)
DIMENSION VECT(LEN)
DATA VECT/LEN*2/, CARRY,WIDTH/0,0/
DO 2 N=1,NUM/4
Q=0
DO 3 L=LEN,1,-1
X=10000*VECT(L)+Q*L
Q=X/(2*L-1)
3 VECT(L)=X-Q*(2*L-1)
X=Q/10000
R=Q-X*10000
TYPE 100, X+CARRY
100 FORMAT('+',I4.4,$)
WIDTH=WIDTH+4
IF (WIDTH.GE.80) TYPE *
IF (WIDTH.GE.80) WIDTH=0
2 CARRY=R
END
|
| 92.23 | | LAMBDA::VOSBURY | | Mon Jul 16 1984 11:39 | 156 |
| { DO NOT PUBLISH WITHOUT STAN'S PERMISSION
What follows is my VAX PASCAL version of Stan's PI algorithm. My intent was
to better display the idea behind the algorithm, so there are some simple
changes that would speed it up some. I one-plussed it so that it will output
PI in any desired radix. (The line to change should be obvious.) The output
format routine is another example of the same idea at work.
What is that idea? Suppose I have a number in [0,1) represented in one number
system (binary, floating point, some 'funny' system) and I want to convert it
into, say, a decimal representation. First I multiply my number by 10 its
number system; i.e. using its notions of '10' and 'multiply'. This gives me a
number in [0,10) in that number system. Now I somehow separate it into its
integer and fractional parts. The integer part, converted into decimal, is
the most significant digit of my decimal representation and the fractional
part, of course, is a another number in [0,1) still expressed in the original
number system. Now I give the handle another crank; 'multiply' the fraction
by '10', separate into integer and fractional parts, save the fraction,
convert the integer into decimal. P.R.N., as the physicians say.
The output routine illustrates this idea, where a binary to 'Radix' conversion
is taking place. Now suppose I had an infinite precision floating point
representation of PI and an infinite precision floating point multiply. Using
this idea, I could output as many decimal places of PI as I liked. The
problem is, how did I get that infinite precision floating point value of PI
to begin with? Suppose instead, I had a number system in which it was easy to
express the value of PI to any degree of precision I wished and easy to
multiply it by a small constant.
In effect, this is what Stan has done. He started with the following series:
PI 1 1 * 2 1 * 2 * 3 1 * 2 * 3 * 4
-- = 1 + - + ----- + --------- + ------------- + ...
2 3 3 * 5 3 * 5 * 7 3 * 5 * 7 * 9
W0 W1 W2 W3 W4 WN
He now creates a mixed radix positional number system, using the terms of the
series as weights; i.e.
a.bcde... = W0 * a + W1 * b + W2 * c + W3 * d + W4 * e + ...
In this system PI has a simple representation:
PI = 2.222222...
A couple of problems remain: how to multiply by a constant and how to make
sure the 'fractional' part is in [0,1). The multiplication part is easy:
N*(a.bcd...) = W0 * (N*a) + W1 * (N*b) + W2 * (N*c) + W3 * (N*d) ...
The second part is handled by putting constraints on the size of the
coefficients in each of the positions. In the decimal system, each coefficient
is constrained to the range 0 to 9. In this system, the coefficient in the
n-th position (n>=1) is constrained to lie in the range 0 to 2*n.
This puts the 'fractional' part in the range [0,2), which is the best that
can be done in this system, given the desire for integer coefficients. This
accounts for the occasional excess in the carry out which showed up as a ':'
in Stan's original program and necessitates keeping the generated digits
around in case they have to be 'fixed-up'.
Multiplication is done as usual from 'right' to 'left' with any overflow in
one position carried over into the next position. Since a '1' in one position
has a different weight than a '1' in the next position, the carry has to be
adjusted by the difference in the two weights.
I've only tested this program for a couple of decimal radix based cases. I
have no way to check, say, a hexadecimal representation. If anyone who does
and cares to check it out, finds that I've made an error, please let me know.
The injunction not to publish this without Stan's permission still applies.
This is his algorithm. I had to explain it to myself and thought I'd pass
along my thoughts. I re-reading what I've written, I notice in a couple of
places I've said things like "Stan did <thus-or-so>." This is, of course
presumptuous on my part. Stan did whatever Stan did. What I've written is my
after the fact reconstruction.
I'll leave you with the following questions:
How did I determine that the 'fractional' part is in the range [0,2)?
(Hint: I didn't formally sum the series.)
What would I have to change in the program in order to generate e, the
base of the natural logarithms.
Does all this shed any light on Note 67?
Mike Vosbury.
}
program StansPI( input, output );
const
Places= 500;
Radix= 10;
Power= 4;
Line_length= 72;
Multiplier= Radix ** Power;
Array_max= ( Places + Power ) div Power;
Array_min= - ( 4 * ( Places + Power ) );
type
Index= Array_min..Array_max;
var
A: array[ Index ] of unsigned;
n, k: Index;
Carry, Modulus, Factor, temp: unsigned;
Column: 0..Line_length-1;
begin
for k:= Array_min to 0 do A[ k ]:= 2;
for k:= 1 to Array_max do A[ k ]:= 0;
for n:= 1 to Array_max do
begin
Carry:= 0;
for k:= -Array_min downto 1 do
begin
Modulus:= 2 * k + 1;
Factor:= k;
temp:= A[ -k ] * Multiplier + Carry;
Carry:= temp div Modulus * Factor;
A[ -k ]:= temp rem Modulus;
end;
for k:= 0 to n do
begin
temp:= A[ k ] * Multiplier + Carry;
Carry:= temp div Multiplier;
A[ k ]:= temp rem Multiplier;
end;
end;
Column:= 0;
for k:= Array_max downto 1 do
for n:= 1 to Power do
begin
temp:= A[ k ] * Radix;
Carry:= temp div Multiplier;
A[ k ]:= temp rem Multiplier;
if (Carry<10)
then write( chr( ord('0') + Carry ) )
else write( chr( ord('A') + Carry - 10 ) );
Column:= ( Column + 1 ) mod Line_length;
if (Column=0) then writeln;
end;
end.
|
| 92.24 | | HARE::STAN | | Sun Jul 22 1984 02:09 | 27 |
| Amazing! Mike has indeed figured out how I arrived at this algorithm.
I didn't think anyone would be able to figure it out from the code.
His analysis is essentially correct.
Now that the jig is out, here are some technical references:
The algorithm was motivated by the simpler algorithm for e that Mike alludes
to. This algorithm circulated around the RT-11 group as folklore
about 8 years ago. I managed to track down the original source as:
A. H. J. Sale, The calculation of e to many significant digits.
Comput. J. 11(Aug 1968)229-230.
I then modified the algorithm to use the series for pi/2 that
was given in the previous note (which I got from Jolley, Summation of
Series) adding in the stuff to handle the fact that the digits produced
were no longer in the range 0 to b-1 (in the base b). I still hope to
publish this sometime.
Algorithms such as this are known as Spigot functions in the literature,
because they produce a digit at a time from the left at (roughly) constant
speed (up to a point).
A more general algorithm for producing values for lots of transcendental
functions by a similar method can be found in Algorithm 393 of the
collected algorithms from CACM: Special Series Summation with Arbitrary
Precision, by S. Kamal Abdali.
|
| 92.25 | | HARE::STAN | | Mon Jul 23 1984 03:23 | 6 |
| According to the Mathematical Intelligencer (vol 6 issue 3, 1984),
Dr. Yasumasa Kanada at the University of Tokyo and his colleagues
have computed pi to 2^24 digits in May 1983. (approximately
16 million digits). They are planning to compute pi to 2^25 digits.
The longest repetition they found in pi was 314159 of length 6 which
repeats 6 times.
|
| 92.26 | | HARE::STAN | | Mon Jul 23 1984 03:30 | 15 |
| It was not clear from the article whether this "longest sequence"
was the longest in length or the longest in most number of times
repeated. It read as if they meant longest in length; but this
couldn't be (see previous response by Peter Gilbert). It is
possible that they meant the longest repetition beginning with
31415926535...
Anyhow, 314159 appears at the following places:
1259351
1761051
6467324
6518294
9753731
9973760
|
| 92.27 | | TURTLE::GILBERT | | Mon Jul 23 1984 17:26 | 2 |
| Certainly, by the pigeon-hole principle, an arbitrarily long repeated sequence
can always be found (given enough digits in the expansion of the number).
|
| 92.28 | | HARE::STAN | | Sun Jul 29 1984 18:12 | 8 |
| Lynn has previously published the first 100 hexadecimal digits of pi in:
Lynn Yarbrough, Precision Calculations of e and pi Constants
Communications of the ACM 10(1967)537.
--------------
Knuth has published 45 octal digits of pi and various other constants
as an appendix to his Art of Computer Programming.
|
| 92.29 | | GOLLY::STAN | | Wed Dec 12 1984 02:30 | 59 |
| From: ROLL::USENET "USENET Newsgroup Distributor" 11-DEC-1984 22:18
To: HARE::STAN
Subj: USENET net.math newsgroup articles
Newsgroups: net.math
Path: decwrl!decvax!tektronix!hplabs!hao!seismo!mcvax!vu44!tstorm
Subject: Re: multiple-precision arithmetic
Posted: Sun Dec 9 08:26:36 1984
Here is another (faster) algorithm to calculate pi. The algorithm that
was posted by mcvax!steven seems to be working with unnecessary
big numbers, although it is in some sense "incremental".
This algorithm is not "incremental".
[The algorithm by mcvax!steven is not posted in MATH.NOT. - STAN -]
:::::::::::::::::::::::::::
First a B implementation of the algorithm:
:::::::::::::::::::::::::::
HOW'TO PI n:
\ Print the first n decimals of pi. mcvax!vu44!tstorm
WRITE 'Approximation of pi in `n` digits' /
WRITE (arcam(48,18,n))+(arcam(32,57,n))+(arcam(-20,239,n)) /
YIELD arcam(a, m, n):
\ Calculate (10**n)*a*arctan(1/m), mcvax!vu44!tstorm
PUT floor(((10**n)*a)/m), -m*m, 0, 1, 0 IN t, mc, s, kk, i
WHILE t<>0:
PUT s+floor(t/kk) IN s
PUT floor(t/mc), kk+2, i+1 IN t, kk, i
IF m=18:
WRITE 'Error < `i+i` units of the last decimal' /
RETURN s+ceiling(i/2)
:::::::::::::::::::::::::::
And now a BC implementation
Old BC versions use the =OP notation the version I'm working with
uses OP= notation. Change line 5 if necessary.
:::::::::::::::::::::::::::
define t(a, m, n){
/* calculate (10^n)*a*arctan(1/m), mcvax!vu44!tstorm */
t = ((10^n)*a)/m; c = -m*m; s = 0; k = 1
while (t != 0){
s += t/k; t /= c; k += 2
}
return(s)
}
define p(n){
"Approximation of pi. Number of digits: "; n
scale = 0; return(t(48, 18, n)+t(32, 57, n)+t(-20, 239, n))
}
:::::::::::::::::::::::::::
--
Theo van der Storm, 52 20'N / 4 52'E, {seismo|decvax}!mcvax!vu44!tstorm
|
| 92.30 | 33 million digits of PI | CLT::GILBERT | it's a word you hear everyday | Sat Jul 19 1986 12:58 | 17 |
| Newsgroups: net.math
Path: decwrl!pyramid!hplabs!hao!ames!eugene
Subject: New record for digits of PI (Japan)
Posted: 17 Jul 86 22:19:23 GMT
Organization: NASA-Ames Research Center, Mtn. View, CA
We have just received a letter from Japan that a newer record for
computation of digits of Pi was accomplished. Previously David Bailey
here at Ames did a 30 million digit computation on the Cray-2.
The new computation was done on an older Hitachi 810 supercomputer
using extended storage. The new record is 33 million digits.
Dave replied, "This means war!"
--eugene miya
NASA Ames Research Center
{hplabs,ihnp4,dual,hao,decwrl,tektronix,allegra}!ames!aurora!eugene
eugene@ames-aurora.ARPA
|
| 92.31 | arctan formulas for PI | CTCADM::ROTH | If you plant ice you'll harvest wind | Tue Jan 03 1989 11:17 | 56 |
| Here's a nice twist on the use of arctan formulas for the computation
of PI - a simple way of getting arbitrarily fast convergence.
It might be interesting to mechanize this in a LISP or MAPLE program to
get high order rational approximations (such as continued fraction
expansions) to PI.
The fastest series (non-AGM) approach to PI I know of is a series
due to Ramunjun - it gives about 10 digits per term, and is not too
complicated.
- Jim
Newsgroups: sci.math
Path: decwrl!sun!pitstop!sundc!seismo!uunet!mcvax!lambert
Subject: Re: Machin's approximation to Pi
Posted: 31 Dec 88 20:54:31 GMT
Organization: CWI, Amsterdam
In article <3696@s.cc.purdue.edu> ags@s.cc.purdue.edu (Dave Seaman) writes:
) >
) > pi 1 1
) > --- = 4 atan(-) - atan(---)
) > 4 5 239
) >
) When you are evaluating speed of convergence, it is the dominant term (1/5,
) in this case), that matters. That term converges much more slowly than the
) term with 1/239.
If the purpose is to produce a decimal representation of pi, the arithmetic
must be done in base 10 (or some power), otherwise the output conversion is
the bottleneck. Then 1/5 really helps in doing the arithmetic. Further
speed-up can be obtained by using
1 1 1
arctan (-) = 2 arctan (---) - arctan (-----------) ,
q 2 q 4 q^3 + 3 q
so that
1 1 1
arctan (-) = 2 arctan (--) - arctan (---) .
5 10 515
This halving trick can be continued to get arbitrarily fast convergence.
I have not given thought to what this means in terms of the actual
computation cost (anyway, `superconvergent' methods to compute pi are
known), but you can obtain, for example,
pi/4 = - arctan(1/239) + 256 arctan(1/320) - 4 arctan(1/515)
- 8 arctan(1/4030) - 16 arctan(1/32060) - 32 arctan(1/256120)
- 64 arctan(1/2048240) - 128 arctan(1/16384480)
--Lambert Meertens, CWI, Amsterdam; lambert@cwi.nl
|
| 92.32 | Well, since you asked... | AKQJ10::YARBROUGH | I prefer Pi | Fri Feb 09 1990 17:49 | 36 |
| OK, Pi lovers; here's the latest in Pi calculation for those of you who
want to burn vax cycles long into the night. The rewards can be Billions of
Pi digits. All you need is a multiprecision arithmetic package and a place to
store the results. A neat part about the method is automatic checkpointing:
you can stop anywhere and resume later without having to repeat anything.
The method was developed by the Chudnovsky brothers, David and Gregory,
according to a report in the January "Discover" Magazine, and is notable
because it works solely with rationals. Only when you want to print
something do you need to perform a long square root, a long division, and a
long data conversion to look at the 3.14... fraction. If you want, and have
the storage to burn, you can do the square root (of a constant) once for
all time.
The algorithm looks like this:
c1 = 13591409/(2*7*9*11*19*127*163)
c2 = ((640320)^3/2)/(8*27*7*11*19*127*163)
For n=0 to as big as you like:
num=(-1)^n*(c1+n)*(6n)! {long integer}
den=(3n)!*(n!^3)*(640320)^3n {Long integer}
sum = sum + num/den {Long rational}
Pi = c2/sum
Each new term in the sum adds about 14 digits of precision to the value of
Pi, and each term is independent of the rest and is the ratio of two
integers. So save what you have of the sum at any stage and restart by
picking up with the last n and the last partial sum.
I've worked it out in MAPLE and it works like a charm. The Chudnovsky's
have used the method to generate BILLIONS of digits of Pi. The only
practical limit to the number of digits is the space to store them.
Lynn Yarbrough
|
| 92.33 | How did they come up with this? | VMSDEV::HALLYB | The Smart Money was on Goliath | Fri Feb 09 1990 19:24 | 6 |
| So c2 starts out with .0078326, right?
Could you show results of the first couple of iterations?
How do you determine how many digits of precision you have? (n*14?)
John
|
| 92.34 | More details on Chudnovsky's method | AKQJ10::YARBROUGH | I prefer Pi | Mon Feb 12 1990 14:41 | 103 |
| Regretably, the article gave neither the origin of the method or the
addresses of the Chudnovskys. I don't know any more about either.
On the following pages are the MAPLE program I used, and its results for
n=0..4.
It's interesting to contemplate how fast this might go on a 3-parallel
processor, if you let 1 calculate num, 2 calculate den, and 3 calculate
sum for the *previous* values of sum, num, and den.
Lynn Yarbrough
# Chudnovsky's approximations [P] to Pi.
c1:= 13591409/(2*7*9*11*19*127*163):
`c1=`, evalf(c1,30);
c2:= ((640320)^(3/2))/(8*27*7*11*19*127*163):
`c2=`, evalf(c2,30);
sum:=0;
for n from 0 to 4 do;
num:=(-1)^n*(c1+n)*(6*n)!:
den:=(3*n)!*(n!^3)*(640320)^(3*n):
sum:= sum + num/den;
P:= evalf(c2/sum,14*(n+1));
`Pi=`, evalf(Pi,14*(n+1));
od;
> # Chudnovsky's approximations [P] to Pi.
c1=, .02493195446879352309804436449729
c2=, .0783260449987954763530818832734
sum := 0
[n=0]
13591409
num := -----------
545140134
den := 1
13591409
sum := -----------
545140134
P := 3.1415926535897
Pi=, 3.1415926535898
[n=1]
22349261720
num := - -------------
30285563
den := 1575224475844608000
386174605303287977741
sum := -------------------------
15489142890368134348800
P := 3.141592653589793238462643383
Pi=, 3.141592653589793238462643383
[n=2]
381498051571200
num := -----------------
393319
den := 397013143887987201995994027786240000000
11265031300432781139756997469231103677
sum := -----------------------------------------
451831055384488536753074001779097600000
P := 3.14159265358979323846264338327950288419717
Pi=, 3.14159265358979323846264338327950288419717
[n=3]
7617308700708016128000
num := - ------------------------
393319
den := 1418372775126561539510626880564903301220881987010560000000000
275771130733485243994916313389977808668129333888013
sum := -------------------------------------------------------
11060951161236225691202257974594411667453378560000000
P := 3.1415926535897932384626433832795028841971693993751058210
Pi=, 3.1415926535897932384626433832795028841971693993751058210
[n=4]
51695833046067011734732800000
num := -------------------------------
20701
den :=
31458317598414079330987587340593945086396822885242319073473492051558400000000000000
213213457884315375924611324013015017490452292719365714507508085126983
sum := ------------------------------------------------------------------------
8551814826679207764567978244525141156892574605621809983258624000000000
P := 3.141592653589793238462643383279502884197169399375105820974944592307817
Pi=, 3.141592653589793238462643383279502884197169399375105820974944592307817
|
| 92.35 | | ALLVAX::ROTH | It's a bush recording... | Mon Feb 12 1990 16:08 | 5 |
| Isn't that a series due to Ramunjun? He published quite a few
astounding series in a paper titled something like "Modular
Equations and Approximations to PI" that is widely quoted.
- Jim
|
| 92.36 | Notes from The Constant Society | ALLVAX::JROTH | It's a bush recording... | Fri Mar 09 1990 13:50 | 140 |
| The following may be of interest - it came from a newsgroup.
The series converging to an integer is surprising; probably representing
it as a continued fraction would explain why...
- Jim
Article 9247
Path: ryn.esg.dec.com!shlump.nac.dec.com!decuac!haven!aplcen!uakari.primate.wisc.edu!zaphod.mps.ohio-state.edu!van-bc!ubc-cs!uw-beaver!milton!aesop
From: aesop@milton.acs.washington.edu (Jeff Boscole)
Newsgroups: sci.physics,sci.math
Subject: 50,000 digits of 'e' (revisited) .. . . . . .
Summary: What is found depends upon what was sought....
Keywords: relative irrationalities, transcendental wiff'n-poof
Message-ID: <2285@milton.acs.washington.edu>
Date: 7 Mar 90 18:21:55 GMT
Reply-To: aesop@milton.acs.washington.edu (Jeff Boscole)
Organization: University of Washington, Seattle
Lines: 154
Xref: ryn.esg.dec.com sci.physics:10900 sci.math:9247
----------------------------------------
In regards the question of relative "irrationalities," for pi & e:
1 2 3
A A A
/ 2 \ / 2 * 4 \ / 4 * 6 * 6 * 8 \
pi = 2 * ( --- ) * ( ------- ) * ( --------------- ) *
\ 1 / \ 3 * 3 / \ 5 * 5 * 7 * 7 /
4
A
/ 8 * 10 * 10 * 12 * 12 * 14 \
( ---------------------------- ) . . . . .
\ 9 * 9 * 11 * 11 * 13 * 13 /
Where A=1 . (J. Wallis, 1616-1703)
-----------------------------------------------------------------------
1 2 3
B B B
/ 2 \ / 2 * 4 \ / 4 * 6 * 6 * 8 \
e = 2 * ( --- ) * ( ------- ) * ( --------------- ) *
\ 1 / \ 3 * 3 / \ 5 * 5 * 7 * 7 /
4
B
/ 8 * 10 * 10 * 12 * 12 * 14 \
( ---------------------------- ) . . . . .
\ 9 * 9 * 11 * 11 * 13 * 13 /
Where B = 1/2 . (N. Pippenger _A.M.M._ May 1980 87:5, p. 391)
-----------------------------------------------------------------------
So both pi & e have a similar representational form in this manner.
An answer for "relative irrationality" between pi & e depends upon how
one approaches the inquiry.
-------------------------------------------------------------------------
Two intriguing results (among many possible) show how pi & e work
harmoniously with each other as :
-2 +3 -4
/ 2 \ / 2 \ / 2 \
pi * e = 6 * ( 1 + --- ) * ( 1 + --- ) * ( 1 + --- ) *
\ 2 / \ 3 / \ 4 /
+5
/ 3 \
( 1 + --- ) . . . . .
\ 5 /
(Z.A. Melzak _A.M.M._ v.68 1961, p. 39 )
----------------------------------------------------------------------
and:
1 2 3 4
pi / 3 \ / 2 \ / 5 \ / 4 \
---- = 2 * ( --- ) * ( --- ) * ( --- ) * ( --- ) . . . . .
e \ 1 / \ 4 / \ 3 / \ 6 /
(Z.A. Melzak _A.M.M._ v.68 1961, p. 39 )
------------------------------------------------------------------------
The preceeding information was brought to you by:
the Constant Society
PO Box 45513
Seattle, WA 98145
or contact: niet@milton.u.washington.edu
-----------------------------------------------------------------------
If 'n' is an odd-square, i.e. n = ( 2k + 1 ) * ( 2k + 1 ), k integer,
then: ( n-1 terms )
n-1
----
\
\ / - ( k * k ) * pi \ is very nearly an integer
S = ) exp ( --------------------- ) , in fact, for n = 9 :
/ \ n / S = ~1.000000000001050
/
----
k=1
-----------------------------------------------------------------------
If 'n' is an even-square, i.e. n = 4 * k * k, with k an integer, then:
n-1
----
\
\ / - ( k * k ) * pi \ is very near to an
S = ) exp ( --------------------- ) , integer plus 1/2 .
/ \ n / For n = 16, S differs
/ from 3.5 in the
---- 22nd decimal place.
k=1
--------------------------------------------------------------------------
-----------------------------------------------------------------------fin
|
| 92.37 | | BEING::EDP | Always mount a scratch monkey. | Mon Jun 29 1992 16:52 | 17 |
| <<< RUSURE::NOTES1:[NOTES$LIBRARY]MATH.NOTE;7 >>>
-< Mathematics at DEC >-
================================================================================
Note 1634.0 Request for C/C++ PI Computation Program: No replies
TALLIS::PERKINS 11 lines 29-JUN-1992 08:46
--------------------------------------------------------------------------------
Does anyone have a C or C++ version of a PI computation program?
I don't care which version of C/C++ or where the program actually runs on,I
can convert it to the system that I want to run it one.
What I like is a version which can do an infinite number of digits.
Either post a reply here or you can send me mail at:
TALLIS::PERKINS
-- Thanks, Eric
|
| 92.38 | PI = 6 ARCSIN (1/2), C version | GVAADG::DUBE | E = MC2 , mainly on Sunday | Thu Jul 23 1992 00:22 | 308 |
|
Hi, Eric,
I have written, on Macintosh, a "Think C" version of a PI computation. As
it uses dynamic memory allocation, the only limitation ( in theory ) is the
longint overflow, which can occur beyond 16 millions digits. Anyway, in practice
I think that you cannot go beyond 100000 on a VAX, as it takes 5 minutes for
10000 digits, about 30 minutes for 100000 digits, and I guess more than 3000 minutes
for 1 million digits. On Macintosh IISi, with Cache-memory set to 1024K, you can
get 10000 digits in 41 minutes.
Please find herebelow the VAX-VMS version of this C program. It just differs from
the Mac version by its long integers printout format : "%-d" on VMS, versus "%li" on Mac.
I use the expansion of PI = 6 Arc Sin (1/2). My very first version of this program was
in PASCAL, which is my second mother language. And as my first mother language is
French, this C version suffers from two iniquities : some French wordings as you will
see, and a kind of direct translation from PASCAL techniques.
Please fill free to use this program, or any idea you could find with it, without any
restriction.
-- Friendly , ##### Remy #####
/* C program to compute Pi to high precision
( in theory, up to 16 millions digits, in you have time enough ... )
( in practice, you can get 10000 digits in 6 minutes on VAX 8500
or 10000 digits in 43 minutes on Macintosh IISi )
Feel free to use part or whole of this program, without any restriction.
##### Remy DUBE #####
simply skip this page if you are just interested with the program.
or have a look at it, patient reader, if you want an overview on the
approach.
Still there ! OK, ready, set, go .....
The method starts from the relation SIN (PI/6) = 1/2
hence PI = 6 ARC_SIN ( 1/2 )
infinity
Which gives the expansion PI = SUM ( U )
i=0 i
with the recurrency :
U = 3
0
U = U * (2i-1)^2 / ( 4 * (2i) * (2i+1) )
i i-1
To demonstrate the above statement, on must compute the
derivative of Arc Sin (x)
d/dx ( ARC_SIN (x) ) = 1 / SQRT ( 1 - x^2 )
infinity
-> We get d/dx ( ARC_SIN x ) = SUM ( V )
i=0 i
with V = (2i)! x^2i / ( (2^2i) * (i!)^2 )
i
the above serie is convergent for ABS (x) < 1
-> within this interval, so for -1 < x < 1 , we may integrate
each term of the expansion, which gives :
infinity
ARC_SIN x = SUM ( W ) + k
i=0 i
with W = (2i)! x^(2i+1) / ( (2^2i) * (i!)^2 * (2i+1) )
i
and k = integration constant = ARC_SIN 0 = 0
=> With x = 0.5, we obtain :
PI = 3 SUM ( (2i)! / ( 2^(4i) * (i!)^2 * (2i+1) )
or PI = SUM (U ), with U = 3
0
and U = U * (2i-1)^2 / ( 4 (2i) (2i+1) )
i i-1
let's designate k the factor (2i-1)^2 / ( 4 (2i) (2i+1) )
i
so : U = U * k
i i-1 i
Now, the following program uses this expansion, with some tricks :
trick No 1 :
----------
as we need only the n first terms of the serie, the SUM is computed as
follows, in descending order :
PI = (...(((k + 3) * k + 3) * .... *) k + 3) * k + 3) * k + 3
n n-1 3 2 1
trick No 2 :
----------
the k factor is decomposed in two terms : k' = (2i-1) / (4i)
k'' = (2i-1) / (4i+2)
trick No 3 :
----------
each element p->value contains a cluster of several digits,
( in the program above, the number of digits in a cluster in called "ensemble" ).
This term "ensemble" is determined dynamically during initialization.
In order to avoid an overflow with long integer calculations, one must have :
log (ensemble ) <= (2^28) * log (4) / precision
10 10
so, with the lowest possible ensemble ( 1 is the lowest ), we get the max
precision that can be given by this program : max precision = 16 millions digits
End_of_tricks
End_of_bla_bla, ouf !
*/
/* !!!!! This is the VAX/VMS version of the program
!!!!! To run this program on a Macintosh, simply replace the printout
!!!!! format "%-d" by the standard "%li" for all long integers printings
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define ZERO 0
#define CLK_TCK 60
#define minimum(a,b) (a<b ? a : b)
typedef struct group
{ long value;
struct group *next;
struct group *previous;
};
extern clock_t clock();
long ensemble;
struct group *p, *first, *last;
long precision; /* number of digits to be computed */
long imax;
long istop;
long weight;
double log_4;
clock_t start, finish, duration;
FILE *file_var;
void ask_number_of_digits()
{ char s[30];
printf("How many digits do you want to compute ? : ");
gets(s);
precision = atol(s);
printf("Computing PI with %-d digits \n\n", precision);
} /* ask_number_of_digits */
void initialise()
{ long i;
log_4 = log10(4);
imax = precision / log_4;
for ((weight=1, ensemble=0); weight <= pow(2,28) / (10*imax); (weight *=10, ensemble++));
for (i=0; i <= (precision+ensemble-1)/ensemble; i++)
{ p = (struct group *) malloc(sizeof(struct group));
p->value = 0;
p->next = ZERO;
if (first == ZERO)
{ first = p;
p->previous = ZERO;
}
else
{ last->next = p;
p->previous = last;
}
last = p;
}
file_var = fopen("pi.dat","w");
fprintf(file_var," PI WITH %-d digits \n \n",precision);
}
void multiply_and_divide ( num, den )
long num;
long den;
{ long temp = 0;
long i ;
for ( (p = first, i = istop+1); i; ( p = p->next, i-- ) )
{ temp *= weight;
temp += num * p->value;
p->value = temp/den;
temp -= den * p->value;
}
}
void compute()
{
long i;
for (i=imax; i; i--)
{ istop = minimum ( (precision - i*log_4 + 2 + ensemble)/ensemble,
(precision+ensemble-1)/ensemble
);
multiply_and_divide(2*i-1, 4*i);
multiply_and_divide(2*i-1, 4*i+2);
first->value += 3;
/* printf("%-d \r",i); !!!!! keep it if you want to see the progession */
}
}
void print_result()
{ long i = 1;
long j = 0;
long count = 0;
long carry = 0;
long temp;
long heures, minutes;
float secondes;
for (p = last; p != ZERO; p = p->previous)
{ p->value = p->value + carry;
carry = 0;
while (p->value >= weight)
{ p->value = p->value - weight;
carry++;
}
}
secondes = (double) (finish - start) / CLK_TCK;
minutes = secondes / 60; secondes = secondes - minutes*60;
heures = minutes / 60; minutes = minutes - heures*60;
fprintf(file_var," Calculation time : %-d hours %-d minutes %f secondes \n \n",
heures, minutes, secondes
);
fprintf(file_var,"%-d.",first->value);
for (p = first->next; p != ZERO; p = p->next)
{ for (j=1; j<=ensemble ; j++)
{ p->value = p->value * 10;
temp = p->value / weight;
p->value = p->value - ( weight * temp);
if (count <= precision) fprintf(file_var,"%-d",temp);
count++;
if (i >= 100)
{ i = 0;
fprintf(file_var,"\n");
fprintf(file_var," ");
}
i++;
}
}
}
main ()
{
ask_number_of_digits();
initialise();
start = clock();
compute();
finish = clock();
print_result();
}
|
| 92.39 | | AUSSIE::GARSON | | Thu Jul 23 1992 03:22 | 10 |
| re .-1
...a lot faster than my first PI computation which was overnight for
10000 digits (in the days when on *a VAX* meant only one thing)...and
as it was I discovered a bug in the output routine which meant that the
result saved in the file was useless. )-:
Seriously, you may want to search in this conference (possibly even
this topic) for the Borwein(?) method for computing PI, which has very
fast convergence.
|
| 92.40 | how to compile .38 | CLARID::DEVAL | | Fri Jul 24 1992 07:45 | 76 |
| should I link pi.c with any library ? is this math.h missing on my system ???
any help appreciated.JcD
$ cc pi
$ link pi
%LINK-W-NUDFSYMS, 10 undefined symbols:
%LINK-I-UDFSYM, ATOL
%LINK-I-UDFSYM, C$MAIN
%LINK-I-UDFSYM, CLOCK
%LINK-I-UDFSYM, FOPEN
%LINK-I-UDFSYM, FPRINTF
%LINK-I-UDFSYM, GETS
%LINK-I-UDFSYM, LOG10
%LINK-I-UDFSYM, MALLOC
%LINK-I-UDFSYM, POW
%LINK-I-UDFSYM, PRINTF
%LINK-W-USEUNDEF, undefined symbol PRINTF referenced
in psect $CODE offset %X00000010
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol GETS referenced
in psect $CODE offset %X0000001A
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol ATOL referenced
in psect $CODE offset %X00000024
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol PRINTF referenced
in psect $CODE offset %X0000003E
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol LOG10 referenced
in psect $CODE offset %X00000078
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol POW referenced
in psect $CODE offset %X000000A5
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol POW referenced
in psect $CODE offset %X000000CC
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol MALLOC referenced
in psect $CODE offset %X000000FC
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol FOPEN referenced
in psect $CODE offset %X00000150
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol FPRINTF referenced
in psect $CODE offset %X00000170
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol FPRINTF referenced
in psect $CODE offset %X00000376
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol FPRINTF referenced
in psect $CODE offset %X0000038C
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol FPRINTF referenced
in psect $CODE offset %X0000039F
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol C$MAIN referenced
in psect $CODE offset %X00000412
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol PRINTF referenced
in psect $CODE offset %X00000422
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol GETS referenced
in psect $CODE offset %X0000042C
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol ATOL referenced
in psect $CODE offset %X00000436
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol PRINTF referenced
in psect $CODE offset %X0000044D
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol CLOCK referenced
in psect $CODE offset %X0000045B
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
%LINK-W-USEUNDEF, undefined symbol CLOCK referenced
in psect $CODE offset %X00000470
in module PI file QUALITY$DISK:[DEVAL]PI.OBJ;1
|
| 92.41 | | BEING::EDP | Always mount a scratch monkey. | Fri Jul 24 1992 12:18 | 18 |
| Re .40:
Those are all standard C functions; you need to include the VAX C
run-time library in your link. I think you can do this with either
LINK PI,SYS$LIBRARY:VAXCRTL.OLB/LIB
or
LINK PI,SYS$LIBRARY:VAXCRTL.EXE/SHARE
The former command includes a copy of the library in your image, so it
can be run on system where the library is not present. The latter
results in a smaller image but requires the system to have the
shareable library available at run time.
-- edp
|
| 92.42 | | AUSSIE::GARSON | | Fri Jul 24 1992 23:23 | 13 |
| re .41
> LINK PI,SYS$LIBRARY:VAXCRTL.EXE/SHARE
Oh, if only this would work...
Actually you need to create a file, say, C.OPT with contents
SYS$LIBRARY:VAXCRTL/SHARE
and then
LINK PI,C/OPT
|
| 92.43 | | WLW::KIER | My grandchildren are the NRA! | Sat Jul 25 1992 23:56 | 7 |
| This is easier:
$LINK PI,SYS$INPUT/OPT
SYS$LIBRARY:VAXCRTL/SHARE
^Z
Mike
|
| 92.44 | | AMCFAC::RABAHY | dtn 456-5689 | Mon Jul 27 1992 20:08 | 13 |
| ... just for fun, I got the following on an Alpha (your mileage may
vary)
PI WITH 10000 digits
Calculation time : 0 hours 0 minutes 48.383335 secondes
3.141592653589793...
Question #1: Is "secondes" simply a typo?
Question #2: For such things, does it make sense to modify the
algorithm for the target platform?
|
| 92.45 | yes, .38 can be adapted for target platform | GVAADG::DUBE | E = MC2 , mainly on Sunday | Tue Jul 28 1992 23:17 | 86 |
|
RE .44
Answer #1 : "secondes" is yes a typo , as the right english spelling is
"seconds". Sorry for this fault. I own you one beer.
However, be aware that these 48 seconds are meant
as elapsed calculation time. So "48 secondes" for 10000
digits on Alpha platform looks a reliable measurement, as
it takes 6 minutes on VAX 8550. One more think to be aware
of : this measurement uses the C function clock(), assuming
that each clock tick corresponds to 1/60 second ( see the
instruction #define CLK_TCK 60 ). Such setting of CLK_TCK
gave a true measurement of the elapsed time, both on
VAX ( 6 minutes ) and on Macintosh ( 41 minutes ).
Now, this clock tick interval can take different values
on other C implementation. For example, in "C , a reference
manual" of Samuel P.Harbison & Guy L.Steel Jr., we have an
example with CLK_TCK defined to 1000000.
Answer #2 : adapting the algorithm in .38 for a specific platform
looks a very good idea, especially for Alpha. I had not yet
the opportunity to work on Alpha, and I know few about it,
except that it's a fast 64 bits chip. Suppose that C
implementation on Alpha can deal with 64 bits integers.
Let's call "super int" such hypothetic integer type. OK,
if "super int" won't exist on Alpha-C, forget about the
following. Now if such "super int" can exist, we just have
to use such type for the variables :
"value" in struct group,
and "temp" in the function "multiply_and_divide".
The overflow calculation then becomes :
log (ensemble ) <= (2^60) * log (4) / precision
10 10
=> in function initialise, one must replace :
"weight <= pow(2,28)/ (10*imax)"
by "weight <= pow(2,60)/ (10*imax)"
This way for a computation of 10000 digits, the cluster of
digits "ensemble" will be increased to 12 digits per
cluster (as compared to 4 digits per cluster with long int).
For 1 million digits, the cluster size will be increased to
10 digits (as compared to 2 digits with long int).
So now, let's evaluate the computation time for 1 million
digits with such configuration, taking "48 secondes with
long int for 10000 digits" for granted.
-> 48 secondes with cluster = 4 indicates that we should get
less than 20 secondes with cluster = 10, always for
10000 digits.
-> for a same cluster size, the computation time varies as
the square of computed digits. So for 1 million digits,
always with cluster = 10, the computation time should be
56 hours.
p.s. : by the way, forget the estimation of timing I gave
in .38 for 1 million digits with "classical VAX". I simply
misused the square law, so I gave there a too optimistic
estimation. Now I own you a second beer, ... hips !
Fun : to compute more than the actual world record ( 33 millions
digits, if I remember correctly ), we have to wait
for "hyper int" integers, based on 128 bits machines.
RE .39 Sure, the algorithm of Borwein, as presented in .12, looks
very faster. So I would like to find somewhere a
demonstration of their formula. Could somebody give an
entry point to their relations ? Thanks in advance.
Challenge ( for myself ) : I will try, with the algorithm in .38, a
computation of PI with 100000 digits on Macintosh. With
cluster group = 3 digits, I forecast 4 days of computation.
I just need one miracle : find some 4 consecutive days
without any power cut in french electrical supply.
-- Friendly
##### Remy #####
|
| 92.46 | | AUSSIE::GARSON | | Wed Jul 29 1992 10:24 | 10 |
| > to compute more than the actual world record ( 33 millions
> digits, if I remember correctly )
I would have thought that the record is higher than this.
> RE .39 Sure, the algorithm of Borwein, as presented in .12, looks
> very faster. So I would like to find somewhere a
> demonstration of their formula.
There appears to be a program in .19 but it uses the Brent MP package.
|
| 92.47 | | BEING::EDP | Always mount a scratch monkey. | Wed Jul 29 1992 12:56 | 7 |
| Re .46:
That note was written 6 years ago; I think the record is now at least a
billion.
-- edp
|
| 92.48 | | AMCFAC::RABAHY | dtn 456-5689 | Wed Jul 29 1992 15:35 | 4 |
| RE: .45
On Alpha, CLOCKS_PER_SEC is 100. So, the correct time is actually ~29
secondes. I haven't tried to use 64 bit ints yet.
|
| 92.49 | 64-bit ints | CX3PT2::KOWTOW::J_MARSH | | Wed Jul 29 1992 17:52 | 6 |
| 92.50 | printing PI computed with 64 bits integers | GVAADG::DUBE | E = MC2 , mainly on Sunday | Thu Jul 30 1992 12:18 | 38 |
| Re .49
>>>You can do 64-bit calculations but I have not yet found a way to print
>>>the results! (This is using DEC C T1.0-045C-2369)
Maybe you can just use the long int formats for printing the result.
In function "print_the_result" ( referring here to algorithm .38 ),
you have to define the variable "temp" as a "long int", not as a 64 bits
int. In fact, such variable "temp" only contains one single digit of PI.
This way, the actual instruction fprintf(file_var,"%-d",temp) could still
be used without any adaptation.
So, for 64 bit calculation, I should suggest the following :
use 64 bits integers for the huge variables :
. value ( in "typedef struct group" )
. weight ( global variable )
. temp ( in function "multiply_and_divide" )
use 32 bits integers for the output variables :
. hours ( in function "print_result" )
. minutes ( in function "print_result" )
. temp ( in function "print_result" )
So, be aware that variable "temp" is :
"64 bits integer" in function "multiply_and_divide" (huge variable)
"32 bits integer" in function "print_result" (printout)
For other integer variables, I think you can use either 32 bits
or 64 bits. However, just in order to avoid a lot of conversions during
the calculation, you have better use 64 bits integers.
Hopefully, this could work properly
##### Remy #####
|
| 92.51 | PI Calculations on a 33-Mhz 80486: | TALLIS::PERKINS | Eric R. Perkins, DTN: 226-6085, (508) 486-6085, TALLIS::PERKINS | Wed Sep 02 1992 15:29 | 20 |
| I took the program given in 92.38 and built a Microsoft C QuickWin
applications. This means that it run as a windows application without
the overheard of dealing with Windows. This also means that the
application get access to all of Windows memory (on my system it is 32MB
physical with anywhere from 4MB to 128MB virtual).
The time I got to do 10,000 was 19 min 3.1 seconds.
Over the weekend, if I have time to get it started I will launch a
1 million digit calculation.
When I get a chance I will also re-build the program for Windows-NT and
see what performance increase there is in going to full 32-bits.
Question:
How can I change the program to go from long to unsigned long. This
should allow me to calculate more digits? or not?
-- Eric
|
| 92.52 | PI Results: | TALLIS::PERKINS | Eric R. Perkins, DTN: 226-6085, (508) 486-6085, TALLIS::PERKINS | Tue Sep 08 1992 13:42 | 23 |
| I ran the program from 92.38 on the 33Mhx 486 running
Windows-NT (32-bit OS) and the results were much better than
under Windows 3.1.
Under NT the program took about 3 min to compute 10K digits.
After looking over the program I notices that for each group that
is allocated, there is 32-bits for the value, and then two pointers.
On a 32-bit OS this means that for each structure there is 32-bit
allocated for the value and then 64-bit for the two pointers.
Wouldn't it be easier just to allocate a huge array to store all
of the 32-bit values?
Also, could someone point out the place where long integer will
overflow causing the 16 million digit limit? It was not
obvious by looking at the code.
I am going to see if I can come up with a C program of the method given
in note 92.32.
-- Eric
|
| 92.53 | | HANNAH::OSMAN | see HANNAH::IGLOO$:[OSMAN]ERIC.VT240 | Tue Sep 08 1992 18:48 | 2 |
| Perhaps more importantly, how did you check your digits ?
|
| 92.54 | :-) | VMSDEV::HALLYB | Fish have no concept of fire. | Tue Sep 08 1992 18:52 | 5 |
| > Perhaps more importantly, how did you check your digits ?
Taking a cue from Nasser, he probably spell-checked them.
John
|
| 92.55 | limit to 16 millions digits on 32 bits machines | GVAADG::DUBE | | Wed Sep 16 1992 23:37 | 113 |
|
re .52
>>Under NT the program took about 3 min to compute 10K digits.
really fast, twice faster than Vax 8550. Anyway for
1 Million digits, it will take 42 days. I will explain
this estimation some paragraphs below.
>>Wouldn't it be easier just to allocate a huge array to store all
>>of the 32-bit values?
yes, indeed. My PASCAL version ran that way. For the C version,
I chose this dynamic allocation using pointeurs for academic
reason, providing a method with no physical limitation, excepted
for the overflow limitation beyond 16 million digits.
>>Also, could someone point out the place where long integer will
>>overflow causing the 16 million digit limit? It was not
>>obvious by looking at the code.
There it is :
the structure allocates 32 bits for the value. Depending on the
precision required, each element gathers several digits, let's say
a "cluster" of digits. Let's call "ensemble" the number of digits in
each cluster.
Now the overflow is to be watched in the function multiply_and_divide.
If you are requiring n digits, the program has to compute a number imax
of terms of the expansion. We get imax = n / log (4) , where "log" means
"decimal logarithm". So imax = n / 0.602
The global variable "weight" is equal to 10^ensemble. For example, if
ensemble = 4, then weight = 10000
We easily can verify, from the calling sequence, that in the function
multiply_and_divide ( num, den ), we have :
num < 2 imax
den < 4 imax
num < den / 2
Now suppose that, before each loop of this function, the 2 following
conditions are satisfied :
(clause 1) : temp < den
(clause 2) : p->value < 2 weight
Let's run the four instructions of the loop :
temp *= weight; {temp < den * weight, now}
temp += num * p-> value; {temp < den * weight + 2 num * weight}
{so : temp < 2 den * weight}
p->value = temp/den; {p->value < 2 weight, so clause 2 verified}
temp -= den * p->value; {temp < den, so clause 1 verified}
So the overflow is to be watched at the second instruction, carrying the
greatest integer value of the computation :
temp < 2 den * weight
with den < 4 imax, or den < 4 n / log(4), we get :
temp < 8 weight * n / log(4)
The best we can do ( for high precision ) is to minimize the cluster
size "ensemble" to ensemble = 1 digit per cluster. This gives
weight = 10. As temp has to remain within signed integers limits
( temp < 2^31 ), this gives, for n, a maximum value which is a little
bit greater than 16 millions digits.
Your remark about unsigned integers is excellent. Using ulong type for
the variable "temp" allows you to compute 32 Millions digits.
You wanted to compute 1 Million digits. You can do it, anyway be aware it
will last for 42 days. Let's evaluate :
3 minutes for 10K digits was your timing. The cluster size depends on
the required precision. For 10K digits, the cluster size was 4 ( you
can verify it using the overflow calculation above). Now for 1 million
digits, you get a cluster size = 2 ( each time you multiply your
precision by 10, then you decrease the cluster size by 1 )
-> so, from a cluster size point of view, your program will be
2 times slower.
-> you will have to compute 100 times more terms of the expansion
-> each term of the expansion will be 100 times greater ( more digits )
So your 3 minutes are finally multiplied by 20000. Which gives 42 days.
I did a little bit less on Macintosh IIsi : 100000 digits, reached
within 96 hours and 35 minutes.
So for several million digits, Borwein formule is better. However,
I cannot understand the origin of such formules, so I feel a little
bit frustrated with them. Could somebody give their mathematical proof ?
re .53 :
>>Perhaps more importantly, how did you check your digits ?
everybody knows them by heart :)
Friendly
##### Remy #####
p.s. : 100K digits to be emailed on request.
|