Saturday, March 12, 2005

Computing Pi using BigDecimal

Introduction

Example program for how to use class java.math.BigDecimal.
See also: Fibonacci Numbers for a program that use the BigDecimal class.

Warning: Source code only compiles with jdk-1_5_0 or later.

Discussion

The program computes π using the Gauss-Legendre Algorithm (Smith 2005).

import java.math.BigDecimal;
import static java.math.BigDecimal.*;

class Pi {
    private static final BigDecimal TWO = new BigDecimal(2);
    private static final BigDecimal FOUR = new BigDecimal(4);
    
    public static void main(String[] args) {
        System.out.println(pi(Integer.parseInt(args[0])));
    }
    
    // Gauss-Legendre Algorithm
    public static BigDecimal pi(final int SCALE) {
        BigDecimal a = ONE;
        BigDecimal b = ONE.divide(sqrt(TWO, SCALE), SCALE, ROUND_HALF_UP);
        BigDecimal t = new BigDecimal(0.25);
        BigDecimal x = ONE;
        BigDecimal y;
        
        while (!a.equals(b)) {
            y = a;
            a = a.add(b).divide(TWO, SCALE, ROUND_HALF_UP);
            b = sqrt(b.multiply(y), SCALE);
            t = t.subtract(x.multiply(y.subtract(a).multiply(y.subtract(a))));
            x = x.multiply(TWO);
        }
        
        return a.add(b).multiply(a.add(b)).divide(t.multiply(FOUR), SCALE, ROUND_HALF_UP);
    }
    
    // the Babylonian square root method (Newton's method)
    public static BigDecimal sqrt(BigDecimal A, final int SCALE) {
        BigDecimal x0 = new BigDecimal("0");
        BigDecimal x1 = new BigDecimal(Math.sqrt(A.doubleValue()));
        
        while (!x0.equals(x1)) {
            x0 = x1;
            x1 = A.divide(x0, SCALE, ROUND_HALF_UP);
            x1 = x1.add(x0);
            x1 = x1.divide(TWO, SCALE, ROUND_HALF_UP);
        }
        
        return x1;
    }
}
[user]$ javac Pi.java
[user]$ java Pi 1000
3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805
321712268066130019278766111959092164202002

Of course, the last few digits are not correct as you may verify by computing more digits. The last 4 digits should be 1989 and not 2002.

The approximation π with 808 places was first obtained by Wrench Ferguson in Sept 1947 using a desk calculator, introducing the calculation of π into the computer age (O'Connor & Robertson 2000).

Square root of 2

Square roots can also be calculated easily. The square root of two can be calculated this way: System.out.println(sqrt(TWO, 1000)); with approximately 1000 places.

1.414213562373095048801688724209698078569671875376948073176679737990732478462107
03885038753432764157273501384623091229702492483605585073721264412149709993583141
32226659275055927557999505011527820605714701095599716059702745345968620147285174
18640889198609552329230484308714321450839762603627995251407989687253396546331808
82964062061525835239505474575028775996172983557522033753185701135437460340849884
71603868999706990048150305440277903164542478230684929369186215805784631115966687
13013015618568987237235288509264861249497715421833420428568606014682472077143585
48741556570696776537202264854470158588016207584749226572260020855844665214583988
93944370926591800311388246468157082630100594858704003186480342194897278290641045
07263688131373985525611732204024509122770022694112757362728049573810896750401836
98683684507257993647290607629969413804756548237289971803268024744206292691248590
52181004459842150591120249441341728531478105803603371077309182869314710171111683
916581726889419758716582152128229518488472

Programming with BigDecimal

You may acces static members by writing import static java.math.BigDecimal.*; This is called static import, and you use it in order to avoid writing BigDecimal.ONE and BigDecimal.ROUND_HALF_UP all over the place. Now, you may just write code such as BigDecimal a = ONE;

You need some constants for 2 and 4 in BigDecimal, so we define them. ONE is already defined (since 1.5):

private static final BigDecimal TWO = new BigDecimal(2);
private static final BigDecimal FOUR = new BigDecimal(4);

You can't use operators and write stuff like x = x * 2; Such code is called syntactic sugar, and would just complicate the language if allowed. Instead, you have to use methods and write x = x.multiply(TWO);

You need to know the number of decimal places, SCALE, when performing division such as BigDecimal b = ONE.divide(sqrt(TWO, SCALE), SCALE, BigDecimal.ROUND_HALF_UP);

The program does not use MathContext objects because it slows the program down as much as 10 times (java version 1.5.0)! You should probably write a.divide(TWO, SCALE, BigDecimal.ROUND_HALF_UP); instead of:

import java.math.MathContext;
...
private static MathContext mc = new MathContext(1000);
...
a.divide(TWO, mc);

Pi to 10000 places

Just to test the code, here is Pi to 10000 places:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805
32171226806613001927876611195909216420198938095257201065485863278865936153381827
96823030195203530185296899577362259941389124972177528347913151557485724245415069
59508295331168617278558890750983817546374649393192550604009277016711390098488240
12858361603563707660104710181942955596198946767837449448255379774726847104047534
64620804668425906949129331367702898915210475216205696602405803815019351125338243
00355876402474964732639141992726042699227967823547816360093417216412199245863150
30286182974555706749838505494588586926995690927210797509302955321165344987202755
96023648066549911988183479775356636980742654252786255181841757467289097777279380
00816470600161452491921732172147723501414419735685481613611573525521334757418494
68438523323907394143334547762416862518983569485562099219222184272550254256887671
79049460165346680498862723279178608578438382796797668145410095388378636095068006
42251252051173929848960841284886269456042419652850222106611863067442786220391949
45047123713786960956364371917287467764657573962413890865832645995813390478027590
09946576407895126946839835259570982582262052248940772671947826848260147699090264
01363944374553050682034962524517493996514314298091906592509372216964615157098583
87410597885959772975498930161753928468138268683868942774155991855925245953959431
04997252468084598727364469584865383673622262609912460805124388439045124413654976
27807977156914359977001296160894416948685558484063534220722258284886481584560285
06016842739452267467678895252138522549954666727823986456596116354886230577456498
03559363456817432411251507606947945109659609402522887971089314566913686722874894
05601015033086179286809208747609178249385890097149096759852613655497818931297848
21682998948722658804857564014270477555132379641451523746234364542858444795265867
82105114135473573952311342716610213596953623144295248493718711014576540359027993
44037420073105785390621983874478084784896833214457138687519435064302184531910484
81005370614680674919278191197939952061419663428754440643745123718192179998391015
91956181467514269123974894090718649423196156794520809514655022523160388193014209
37621378559566389377870830390697920773467221825625996615014215030680384477345492
02605414665925201497442850732518666002132434088190710486331734649651453905796268
56100550810665879699816357473638405257145910289706414011097120628043903975951567
71577004203378699360072305587631763594218731251471205329281918261861258673215791
98414848829164470609575270695722091756711672291098169091528017350671274858322287
18352093539657251210835791513698820914442100675103346711031412671113699086585163
98315019701651511685171437657618351556508849099898599823873455283316355076479185
35893226185489632132933089857064204675259070915481416549859461637180270981994309
92448895757128289059232332609729971208443357326548938239119325974636673058360414
28138830320382490375898524374417029132765618093773444030707469211201913020330380
19762110110044929321516084244485963766983895228684783123552658213144957685726243
34418930396864262434107732269780280731891544110104468232527162010526522721116603
96665573092547110557853763466820653109896526918620564769312570586356620185581007
29360659876486117910453348850346113657686753249441668039626579787718556084552965
41266540853061434443185867697514566140680070023787765913440171274947042056223053
89945613140711270004078547332699390814546646458807972708266830634328587856983052
35808933065757406795457163775254202114955761581400250126228594130216471550979259
23099079654737612551765675135751782966645477917450112996148903046399471329621073
40437518957359614589019389713111790429782856475032031986915140287080859904801094
12147221317947647772622414254854540332157185306142288137585043063321751829798662
23717215916077166925474873898665494945011465406284336639379003976926567214638530
67360965712091807638327166416274888800786925602902284721040317211860820419000422
96617119637792133757511495950156604963186294726547364252308177036751590673502350
72835405670403867435136222247715891504953098444893330963408780769325993978054193
41447377441842631298608099888687413260472156951623965864573021631598193195167353
81297416772947867242292465436680098067692823828068996400482435403701416314965897
94092432378969070697794223625082216889573837986230015937764716512289357860158816
17557829735233446042815126272037343146531977774160319906655418763979293344195215
41341899485444734567383162499341913181480927777103863877343177207545654532207770
92120190516609628049092636019759882816133231666365286193266863360627356763035447
76280350450777235547105859548702790814356240145171806246436267945612753181340783
30336254232783944975382437205835311477119926063813346776879695970309833913077109
87040859133746414428227726346594704745878477872019277152807317679077071572134447
30605700733492436931138350493163128404251219256517980694113528013147013047816437
88518529092854520116583934196562134914341595625865865570552690496520985803385072
24264829397285847831630577775606888764462482468579260395352773480304802900587607
58251047470916439613626760449256274204208320856611906254543372131535958450687724
60290161876679524061634252257719542916299193064553779914037340432875262888963995
87947572917464263574552540790914513571113694109119393251910760208252026187985318
87705842972591677813149699009019211697173727847684726860849003377024242916513005
00516832336435038951702989392233451722013812806965011784408745196012122859937162
31301711444846409038906449544400619869075485160263275052983491874078668088183385
10228334508504860825039302133219715518430635455007668282949304137765527939751754
61395398468339363830474611996653858153842056853386218672523340283087112328278921
25077126294632295639898989358211674562701021835646220134967151881909730381198004
97340723961036854066431939509790190699639552453005450580685501956730229219139339
18568034490398205955100226353536192041994745538593810234395544959778377902374216
17271117236434354394782218185286240851400666044332588856986705431547069657474585
50332323342107301545940516553790686627333799585115625784322988273723198987571415
95781119635833005940873068121602876496286744604774649159950549737425626901049037
78198683593814657412680492564879855614537234786733039046883834363465537949864192
70563872931748723320837601123029911367938627089438799362016295154133714248928307
22012690147546684765357616477379467520049075715552781965362132392640616013635815
59074220202031872776052772190055614842555187925303435139844253223415762336106425
06390497500865627109535919465897514131034822769306247435363256916078154781811528
43667957061108615331504452127473924544945423682886061340841486377670096120715124
91404302725386076482363414334623518975766452164137679690314950191085759844239198
62916421939949072362346468441173940326591840443780513338945257423995082965912285
08555821572503107125701266830240292952522011872676756220415420516184163484756516
99981161410100299607838690929160302884002691041407928862150784245167090870006992
82120660418371806535567252532567532861291042487761825829765157959847035622262934
86003415872298053498965022629174878820273420922224533985626476691490556284250391
27577102840279980663658254889264880254566101729670266407655904290994568150652653
05371829412703369313785178609040708667114965583434347693385781711386455873678123
01458768712660348913909562009939361031029161615288138437909904231747336394804575
93149314052976347574811935670911013775172100803155902485309066920376719220332290
94334676851422144773793937517034436619910403375111735471918550464490263655128162
28824462575916333039107225383742182140883508657391771509682887478265699599574490
66175834413752239709683408005355984917541738188399944697486762655165827658483588
45314277568790029095170283529716344562129640435231176006651012412006597558512761
78583829204197484423608007193045761893234922927965019875187212726750798125547095
89045563579212210333466974992356302549478024901141952123828153091140790738602515
22742995818072471625916685451333123948049470791191532673430282441860414263639548
00044800267049624820179289647669758318327131425170296923488962766844032326092752
49603579964692565049368183609003238092934595889706953653494060340216654437558900
45632882250545255640564482465151875471196218443965825337543885690941130315095261
79378002974120766514793942590298969594699556576121865619673378623625612521632086
28692221032748892186543648022967807057656151446320469279068212073883778142335628
23608963208068222468012248261177185896381409183903673672220888321513755600372798
39400415297002878307667094447456013455641725437090697939612257142989467154357846
87886144458123145935719849225284716050492212424701412147805734551050080190869960
33027634787081081754501193071412233908663938339529425786905076431006383519834389
34159613185434754649556978103829309716465143840700707360411237359984345225161050
70270562352660127648483084076118301305279320542746286540360367453286510570658748
82256981579367897669742205750596834408697350201410206723585020072452256326513410
55924019027421624843914035998953539459094407046912091409387001264560016237428802
10927645793106579229552498872758461012648369998922569596881592056001016552563756
78

References

No comments: