Algorithms using real numbers Floating point arithmetic limited - - PDF document

Algorithms using “real” numbers

• floating point arithmetic
• typical computational errors due to limited precision arithmetic
• advice for constructing algorithms using real numbers

1

Algorithms using “real” numbers

• Floating point arithmetic

– limited precision – computational errors – IEEE standard – type double in Java

• Computational errors due to limited precision arithmetic

– typical errors involving summation, convergence, subtraction – how to avoid some of them

• Applications of numerical algorithms

– Ex. π-computation – graphical / geometric computations – scientific computations

2

Floating point arithmetic: limited precision

• A computer can handle only a finite subset of the reals directly

in hardware

• These numbers make up a floating point number system

F(β, s, m, M) characterized by a base β ∈ N\{1} a number of digits s ∈ N a smallest exponent m ∈ Z largest exponent M ∈ Z

• Each floating point number has the form

y = ±d1.d2 . . . ds · βe, with m ≤ e ≤ M and 0 ≤ dk ≤ β − 1.

3

Floating point arithmetic: limited precision

• Each floating point number in F(β, s, m, M) has the form

y = ±d1.d2 . . . ds · βe, with m ≤ e ≤ M and 0 ≤ dk ≤ β − 1.

• The mantissa d1.d2 . . . ds represents d1 +d2 ·β−1 +· · ·+ds ·β−s+1
• the exponent part is βe
• The system includes 0 = 0.0 . . . 0 · βe.
• All other numbers are normalised: 1 ≤ d1 ≤ β − 1

4

Floating point arithmetic: limited precision

• Example: the number system F(2, 3, −2, 1) contains 33 numbers:

−1.11 · 21 . . . −1.00 · 2−2 0.00 1.00 · 2−2 . . . 1.11 · 21

5

Floating point arithmetic: computational errors

• standard demo system is F(10, 4, −99, 99)
• 1.573 og 0.1824 are representable, but the following are not

1.573 + 0.1824 = 1.7554 1.573 − 0.1824 = 1.3906 1.573 × 0.1824 = 0.2869152 1.573 / 0.1824 = 8.6239035 . . .

• Exact arithmetic on a finite subset of reals is not possible.
• strategy for rounding/truncating to a representable number is

needed

6

Floating point arithmetic: computational errors

• rounding/truncation defined by function fl : R → M, where M

is machine numbers

• machine arithmetic operations ⊕, ⊖, ⊗, ⊘ defined by

x ⊕ y = fl(x + y) etc.

• In demo system F(10, 4, −99, 99) define fl by truncation, e.g.

1.573 ⊕ 0.1824 = fl(1.573 + 0.1824) = fl(1.7554) = 1.755 1.573 ⊖ 0.1824 = fl(1.573 − 0.1824) = fl(1.3906) = 1.390 1.573 ⊗ 0.1824 = fl(1.573 × 0.1824) = fl(0.2869152) = .2869 1.573 ⊘ 0.1824 = fl(1.573 / 0.1824) = fl(8.6239035 . . .) = 8.623

7

Floating point arithmetic: computational errors

• Algebraic laws invalid:

(1.418 ⊕ 2937) ⊖ 2936 = 2938 ⊖ 2936 = 2.000 1.418 ⊕ (2937 ⊖ 2936) = 1.418 ⊕ 1.000 = 2.418 1.418 ⊗ (2001 ⊖ 2000) = 1.418 ⊗ 1.000 = 1.418 (1.418 ⊗ 2001) ⊖ (1.418 ⊗ 2000) = 2837 ⊖ 2836 = 1.000

• The usual associative and distributive laws are not valid for

machine arithmetic. Sometimes: (a ⊕ b) ⊖ c = a ⊕ (b ⊖ c) a ⊗ (b ⊖ c) = (a ⊗ b) ⊖ (a ⊗ c)

8

IEEE standard for binary floating point arithmetic

• IEEE standard describes two number systems, approximately

Single precision F(2, 24, −126, 127) Double precision F(2, 53, −1022, 1023)

• The IEEE systems has more numbers:

– closes the representational gap around zero with xtra numbers

• f the form 0.b1b2 . . . b522−1022

– has representations for ±∞ and NaN (= Not a Number)

• The IEEE system has detailed rules for rounding to representable

numbers, ex. – overflow: 1/0 or 21023 · 2 has result ∞ – underflow: 1/∞ or 2−1022−52/2 has result 0 – 0/0 or ∞ − ∞ has result NaN

9

Floating point arithmetic in Java

• Java implementations use IEEE standard for the primitive types

float and double on most computers

• Use modifier strictfp to enforce IEEE standard independent of

the underlying hardware.

• Representable values in type double are

±b1.b2 . . . b53 · 2e, where bi ∈ {0, 1} and −1022 ≤ e ≤ 1023 ±∞ and NaN

• Java syntax for ±∞ and NaN is: Double.NEGATIVE INFINITY,

Double.POSITIVE INFINITY and Double.NaN

• Warning: legal java double literals may not be directly

representable: literal 0.2 is represented as

fl(0.210) = fl(0.00110011 . . .2) = 1.1001100110011001100110011001100110011001100110011010 · 2−3

10

Sources of errors in numerical computations

• Hardware errors

– Pentium bug

• Limited precision arithmetic

– summation: adding numbers of unequal magnitude – subtraction: catastrophic cancellation – convergence: may cycle or diverge – formula: mathematically equivalent formulae have different computational properties

• Logical errors

11

Limited precision arithmetic

• The limited precision of arithmetic in R(10, 4, −99, 99) is

illustrated by computing e = lim

n→∞(1 + 1

n)n k 1 ⊘ 10k 1 ⊕ (1 ⊘ 10k) (1 ⊕ 1 ⊘ 10k)10k 1 1.000E-1 1.100 2.591 2 1.000E-2 1.010 2.591 3 1.000E-3 1.001 2.591 4 1.000E-4 1.000 1.000 5 1.000E-5 1.000 1.000

12

Limited precision arithmetic

• The e-computation in R(2, 24, −126, 127):

k 1 ⊘ 2k 1 ⊕ (1 ⊘ 2k) (1 ⊕ 1 ⊘ 2k)2k 1 0.5 1.5 2.25 2 0.25 1.25 2.4414062 3 0.125 1.125 2.5657845 4 0.0625 1.0625 2.6379285 . . . . . . . . . . . . 22 2.3841858E-7 1.0000002 2.7175977 23 1.1920929E-7 1.0000001 2.7175977 24 5.9604645E-8 1.0 1.0 25 2.9802322E-8 1.0 1.0

13

Limited precision arithmetic

• The e-computation seems to give a reasonable result, if one stops

just before the result degenerates to 1?

• WRONG:

e R(10, 4, −99, 99) 2 .591 R(2, 24, −126, 127) 2.71 75977 R(2, 53, −1022, 1023) 2.7182818 08182473 exact 2.71828182845904523536 ...

• less than half the significant digits are correct!

14

Summation

• The nth harmonic number is defined by

Hn =

n

• i=1

1 i = 1 + 1 2 + 1 3 + · · · + 1 n

• Associative law is not valid.
• Which order of summation is best?

– from the left? – from the right? – pairwise in binary tree?

• experiment in toy system R(10, 4, −99, 99) with truncation

15

Summation

• try forward summation:

Hn = (· · · ((1 ⊕ 1

2) ⊕ 1 3) ⊕ · · · ⊕ 1 n)

n Hn exact ⊕, ⊖, ⊗, ⊘ forward 10 2.928 2.927 100 5.187 5.142 1000 7.485 7.069 10000 9.787 7.069 100000 12.09 7.069 1000000 14.39 7.069

• The computed sum is constant for n ≥ 1000 since

7.069 ⊕ 1 1001 = 7.069 ⊕ 9.990 · 10−4 = 7.069

16

Summation

• try backward summation:

Hn = (1 ⊕ ( 1

2 ⊕ ( 1 3 ⊕ · · · ⊕ 1 n) · · ·))

n Hn exact ⊕, ⊖, ⊗, ⊘ backward forward 10 2.928 2.928 2.927 100 5.187 5.170 5.142 1000 7.485 7.284 7.069 10000 9.787 8.069 7.069 100000 12.09 8.069 7.069 1000000 14.39 8.069 7.069

• By adding terms in increasing order, the error is reduced!

17

Summation

• add terms pairwise in a tree structure

Hn = (· · · ((1 ⊕ 1

2) ⊕ ( 1 3 ⊕ 1 4⊕)) · · · ⊕ 1 n) · · ·))

1

1 2 1 3 1 5 1 7

⊕ ⊕ ⊕ ⊕ ⊕ ⊕

1 4

⊕

1 6 1 8

18

Summation

n Hn exact ⊕, ⊖, ⊗, ⊘ tree-wise backward forward 10 2.928 2.928 2.928 2.927 100 5.187 5.183 5.170 5.142 1000 7.485 7.477 7.284 7.069 10000 9.787 9.778 8.069 7.069 100000 12.09 12.07 8.069 7.069 1000000 14.39 14.36 8.069 7.069

• By adding equal sized terms the error is diminished further!

19

Summation

• similar experiment for java double, R(2, 53, −1022, 1023)

H100000 error +, −, ·, / exact 12.090146129863428 ⊕, ⊖, ⊗, ⊘ “tree-wise” 12.0901461298634 31 −3 · 10−15 ⊕, ⊖, ⊗, ⊘ backward 12.0901461298634 08 20 · 10−15 ⊕, ⊖, ⊗, ⊘ forward 12.090146129863 335 93 · 10−15

• seems better than R(10, 4, −99, 99) ?
• truncation errors add up, but rounding errors may cancel out

(though there is no guarantee!)

20

Summation

• Problem

– compute x1 + x2 + · · · + xn

• Warning:

– ⊕ not associative!

• Advice:

– Add terms in increasing order – Even better: add terms of equal size only

21

Convergence: monotone

• Example: find root of polynomial f(x) using Newton iteration
• compute z = limk→∞ xk for

xk+1 = xk − f(xk) f ′(xk)

f(xk) xk+1 z y = f(x) α xk

22

Convergence: monotone

• Example: squareroot computation, √a
• Find root of f(x) = x2 − a.
• Newton iteration: compute

xk+1 = xk − x2

k − a

2xk using x0 = 1 2(a + 1)

• Does it hold that

lim

k→∞ xk = √a

for limited precision arithmetic?

23

Convergence: monotone

• Computing

√ 2 by Newton iteration R(10, 4, −99, 99) R(2, 53, −1022, 1023) k xk xk 1.500

• 1. 5000000000000000

1 1.416 1.41 66666666666667 2 1.414 1.41421 56862745099 3 1.414 1.41421356237 46899 4 1.414 1.414213562373095 1 5 1.414 1.4142135623730950 6 1.414 1.414213562373095 1 exact √ 2 1.414 1.4142135623730950

• the computation in R(2, 53, −1022, 1023) seems to be cycling?

24

Convergence: monotone

• When should iteration be stopped?
• naive criteria: xk+1 = xk

– WRONG: cycling computation will continue forever

• better criteria: xk+1 is sufficiently close to xk

– what does sufficiently close mean?

25

Convergence: monotone

• Convergence: √a = limk→∞ xk,

for x0 = 1 2(a + 1) xk+1 = xk − x2

k − a

2xk

• Monotone convergence:

x0 > x1 > . . . > xk > xk+1 > . . . > √a

• Approximation (for 0 < ǫ < 1

4):

xk − xk+1 < ǫ · xk ⇒ xk − √a < ǫ · (1 + 3ǫ) · √a ≈ ǫ · √a

• Exact arithmetic: √a is computed with precision ǫ using stop

criteria xk − xk+1 < ǫ · xk

26

Convergence: monotone

• In R(β, s, m, M) the optimal precision is ǫ = β−s, giving the

stopping criteria xk − xk+1 < β−s · xk. This is actually equivalent to ¬(xk > xk+1) for the square root monotone convergence.

• double sqrt(double a) {

xnew = (a + 1)/2; xold = xnew + 1; while (xold > xnew) { xold = xnew; xnew = xold − (x2

• ld − a)/2/xold;}

return xold; }

27

Convergence: monotone

• iteration stops when sequence no longer decreases

R(10, 4, −99, 99) R(2, 53, −1022, 1023) k xk xk 1.500

• 1. 5000000000000000

1 1.416 1.41 66666666666667 2 1.414 1.41421 56862745099 3 1.414 1.41421356237 46899 4 stop 1.414213562373095 1 5 1.4142135623730950 6 1.414213562373095 1 7 stop exact √ 2 1.414 1.4142135623730950

28

Convergence: nonmonotone

• Example: find root of equation 3x2 + 8x − 10 = 0 by iterating

towards fixed point

• Rewrite equation: x = 1

8(10 − 3x2)

• The two roots, z1 and z2 are fixed points for iteration:

xk+1 = 1 8(10 − 3x2

k)

2 4 −4 −2 z2 x0 ∈ x0 = x0 ∈ limk→∞ xk = −∞ limk→∞ xk = z1 limk→∞ xk = z2 −z1 z1

29

Convergence: nonmonotone

• Convergence: z = limk→∞ xk

is root, for x0 = xk+1 = 1 8(10 − 3x2

k)

• Nonmonotone convergence:

x0 < x2 . . . < x2t < . . . < z < . . . < x2t+1 < . . . < x3 < x1

• Approximation (for 0 < ǫ < 1

2):

|xk − xk+1| < ǫ · xk ⇒ |xk − z| < ǫ · (1 + 2ǫ) · z ≈ ǫ · z

• Exact arithmetic: z is computed with precision ǫ using stop

criteria |xk − xk+1| < ǫ · xk

30

Convergence: nonmonotone

• In R(β, s, m, M) the optimal precision is ǫ = β−s, giving the

stopping criteria |xk − xk+1| < β−s · xk. This is equivalent to ¬(|xk − xk+1| > 0) .

• double root(double a) {

xnew = 0; xold = xnew + 1; while (|xold − xnew| > 0) //WRONG! { xold = xnew; xnew = (10 − 3x2

• ld)/8;}

return xk; }

• WRONG: may loop forever!

31 R(10, 4, −99, 99) R(2, 53, −1022, 1023) k xk k xk 0.000 0.0 1 1.250 1 1.25 2 .6642 2 0.6640625 3 1.084 3 1.0846328735351562 4 .8093 4 0.8088393236175762 5 1.004 5 1.0046671057136982 . . . . . . . . . . . . 20 .9271 96 0.9274433277084224 21 .9277 97 0.9274433277084229 22 .9273 98 0.9274433277084226 23 .9276 99 0.9274433277084227 24 .9273 100 0.9274433277084226 25 .9276 101 0.9274433277084227 26 .9273 102 0.9274433277084226 27 .9276 103 0.9274433277084227 . . . . . . . . . . . . exact z .9274 0.9274433277084227 32

Convergence

• Problem

– compute limk→∞ xk using xk+1 = f(xk)

• Warning:

– Using limited precision arithmetic, the series may fail to converge – Naive stopping criteria may lead to looping forever!

• Advice:

– stop iteration when adjacent terms are sufficiently close – the meaning of “sufficiently close” depends on the concrete series – for monotone convergent series it may be ok to loop until strict monotonicity fails

33

Formulae and subtraction

• Example: computing π by convergent series
• π = limk→∞ pk, for pk being half the circumference of regular

k-gon inscribed in unit circle. s6 = 1 s12 ≈ .52 p6 = 3 p12 ≈ 3.1

• Monotonicity: p6 < p12 < p24 < . . . < pk < p2k < . . . < π
• Enough to compute side length sk, since pk = k·sk

2 34

Formulae and subtraction

• computing s2k from sk

sk d s2k

35

Formulae and subtraction

• From Pythagoras:

(sk/2)2 + d2 = s2

2k

(sk/2)2 + (1 − d)2 = 1 leading to s6 = 1 s2k =

• 2 −
• (2 − sk)(2 + sk)

and π = lim

k→∞ pk,

for pk = k · sk 2

36

Formulae and subtraction

• finite precision algorithm:

double BADpi() { snew = 1; sold = ..; k = 6; pnew = 3; pold = pnew − 1; while (pnew > pold) { pold = pnew; sold = snew; snew =

• 2 −
• (2 − sold)(2 + sold);

//BAD! pnew = k · snew/2; k = 2k; } return pold; }

37

k

• (2 − sold)(2 + sold)

2 − √· · · snew pnew 6 − − 1 3 12 1.732 .2680 .5177 3.106 24 1.931 .6900 .2626 3.151 48 1.982 .01800 .1342 3.220 96 1.994 .006000 .07746 3.718 192 1.997 .003000 .05477 5.255 384 1.998 .002000 .04473 8.585 768 1.999 .001000 .03162 12.14 1536 1.999 .001000 .03162 24.28 3072 1.999 .001000 .03162 48.56 6144 1.999 .001000 .03162 97.10 . . . . . . . . . . . . . . .

38

Formulae and subtraction

• What goes wrong?

– it seems pk → ∞ rather than pk → π –

• (2 − sold)(2 + sold) is so close to 2 that the subsequent

subtraction 2 − √· · · leads to a catastrophic cancellation. – The stopping criteria ¬(pnew > pold) doesn’t work.

39

Formulae and subtraction

• Avoid the catastrophic cancellation by rewriting the recursion

formula using x − y = (x + y)(x − y) x + y = x2 − y2 x + y

• i.e.

s2k =

• 2 −
• (2 − sk)(2 + sk) =
• 2 −
• 4 − s2

k

=

• 4−(4−s2

2+√ 4−s2

=

sk

• 2+√

4−s2

• leading to

s2k = sk

• 2 +
• (2 − sk)(2 + sk)

40

Formulae and subtraction

• finite precision algorithm:

double GOODpi() { snew = 1; sold = ..; k = 6; pnew = 3; pold = pnew − 1; while (pnew > pold) { pold = pnew; sold = snew; snew = sold/

• 2 +
• (2 − sold)(2 + sold) ;

pnew = k · snew/2; k = 2k; } return pold; }

41

Formulae and subtraction

k

• (2 − sold)(2 + sold)

2 + √· · · snew pnew 6 − − 1 3 12 1.732 3.732 .5175 3.105 24 1.931 3.931 .2610 3.132 48 1.982 3.982 .1308 3.139 96 1.995 3.995 .06546 3.142 192 1.998 3.998 .03274 3.143 384 1.999 3.999 .01637 3.143 stop

• Correct value π = 3.141592653589793238462643 . . .

42

Formulae and subtraction

• Algorithm BADpi in R(53,2,-1022,1023):

– The stopping criteria ¬(pnew > pold) makes it stop, but not until the value is far past π. (6 out of 17 digits correct) – More iterations does not improve the approximation!

• Algorithm GOODpi in R(53,2,-1022,1023):

– The stopping criteria ¬(pnew > pold) makes it stop, and the value is only slightly past the true value of π.

43 Algorithm BADpi k pk k pk 6 3 12 3.1 058285412302498 786432 3.141592 3038117377 24 3.1 32628613281237 stop - but if continued : 48 3.1 39350203046872 1572864 3.141 6086962248038 96 3.141 03195089053 3145728 3.1415 868396550413 192 3.141 4524722856703 6291456 3.141 6742650217575 384 3.1415 57607912925 12582912 3.141 6742650217575 768 3.1415 83892148936 25165824 3.14 30727401700396 1536 3.14159 04632367617 50331648 3.1 37475099502783 3072 3.141592 1059596714 100663296 3 .092329219213245 6144 3.141592 5165881546 201326592 3 .3541019662496847 12288 3.1415926 186407894 402653184 4.242640687119286 24576 3.1415926 453212157 805306368 6.0 49152 3.1415926 666655567 1610612736 0.0 98304 3.141592 730698579 3221225472 0.0 196608 3.141592 9868306565 6442450944 0.0 393216 3.14159 3669849427 . . . . . . 44 Algorithm GOODpi k pk k pk 6 3 98304 3.141592653 055036 12 3.1 05828541230249 196608 3.141592653 4561036 24 3.1 326286132812378 393216 3.1415926535 56371 48 3.1 393502030468667 786432 3.14159265358 1438 96 3.141 0319508905093 1572864 3.14159265358 77046 192 3.141 4524722854615 3145728 3.141592653589 2713 384 3.1415 576079118575 6291456 3.141592653589 663 768 3.1415 838921483177 12582912 3.1415926535897 61 1536 3.14159 04632280496 25165824 3.1415926535897 856 3072 3.141592 105999271 50331648 3.14159265358979 13 6144 3.141592 5166921563 100663296 3.141592653589793 6 12288 3.1415926 19365383 201326592 3.14159265358979 4 24576 3.1415926 450336897 402653184 3.14159265358979 4 49152 3.14159265 14507664 stop 45

Formulae and subtraction

• Problem

– in some computation two terms of approximately identical size are subtracted: E = E1 − E2, where E1 ≈ E2

• Warning:

– E will be computed with few or no significant digits – results are corrupted and loops may not terminate

• Advice:

– rewrite formula to avoid catastrophic cancellation – sometimes the following simple reformulation is useful: x − y = (x + y)(x − y) x + y = x2 − y2 x + y

46

computing π

• Buffon’s needle (Le Clerc, 1777)
• In/circumscribed polygon (Archimedes, 287-212, B.C.)
• Arctan series (Machin, 1706)
• Modular equations (Ramanujan, 1914; Bayley, Borwein and

Borwein 1989)

47

computing π: Buffon’s needle

(Le Clerc, 1777)

• Throw some needles of unit length on a floor covered with planks
• f unit width
• Count those needles that land across a crack between two planks

Needles: Wooden floor:

48

computing π: Buffon’s needle

• Assuming needles are rotated uniformly randomly, the

probability of a needle landing across a gap is 2

π.

π ≈ 2 · #needles #crossings

• Experiment:

π ≈ 2 · 250000 158847 = 3.14 7682991

• Applet simulator at

http://www.math.csusb.edu/faculty/stanton/probstat/buffon.html

49

computing π: inscribed/circumscribed polygon

Archimedes (287-212 B.C)

• Let an (and bn) be length of circumscribed (and inscribed)

regular 6 · 2n-gons about a circle of diameter 1. bn < π < an

• Mutually recursive formulae:

a0 = 2 √ 3 b0 = 3 an+1 = 2anbn an + bn bn+1 =

• an+1bn
• Archimedes obtained π ≈ 22

7 using n = 4

• Ludolph van Ceulen (1540-1610) calculated 34 digits of π

50

computing π: inscribed/circumscribed polygon

n bn an 3 .0 3 .4641016151377544 1 3.1 058285412302493 3 .215390309173473 2 3.1 32628613281238 3.1 596599420975005 3 3.1 39350203046867 3.14 6086215131435 4 3.141 0319508905093 3.14 27145996453683 5 3.141 452472285462 3.141 873049979824 . . . . . . . . . 23 3.14159265358979 3.14159265358979 62 24 3.14159265358979 13 3.141592653589793 25 3.14159265358979 2 3.14159265358979 22 26 3.14159265358979 2 3.14159265358979 2

51

computing π: Arctan series

Machin, 1706 π 4 = 4 arctan(1 5) − arctan( 1 239) arctan x = x − x3 3 + x5 5 − · · · for |x| ≤ 1

• Machin (1706) calculated 100 digits of π.
• Guilloud and Bouyer (1973) on CDC 7600 computed more than 1

mio digits of π.

52

computing π: Arctan series

n 4n

i=0 (−1)i 2i+1 ( 4 52i+1 − 1 2392i+1 )

3.1 8326359832636 1 3.14 05970293260608 2 3.141 621029325035 3 3.14159 17721821778 4 3.1415926 824044 5 3.14159265 2615309 6 3.141592653 6235554 7 3.14159265358 8603 8 3.141592653589 8366 9 3.14159265358979 27 10 3.14159265358979 44 11 3.14159265358979 44

53

computing π: modular equations

(Ramanujan, 1914; Bayley, Borwein and Borwein 1989) lim

n→∞ αn = 1

π for mutually recursive formulae: α0 = 6 − 4 √ 2 y0 = √ 2 − 1 yn+1 = 1 − (1 − y4

n)1/4

1 + (1 − y4

n)1/4

αn+1 = (1 + yn+1)4αn − 22n+3yn+1(1 + yn+1 + y2

n+1)

• Bayley (1986) computed more than 29 mio digits of π.
• Kanada (1996) computed more than 6400 mio digits of π.

54

computing π: modular equations

n yn

1 αn

0.41421356237309515 2.9142135623730985 1 0.003734885463325165 3.1415926 46213547 2 2.4323154602138167E-11 3.141592653589 8095 3 0.0 3.141592653589 8095 4 0.0 3.141592653589 8095

• Only 13 of 17 digits are correct: The formula leads to a

catastrophic cancellation

55

computing π: modular equations

• rewrite

yn+1 = 1 − (1 − y4

n)1/4

1 + (1 − y4

n)1/4

using x − y = (x + y)(x − y) x + y = x2 − y2 x + y

• The result is

yn+1 = y4

n

(1 + (1 − y4

n)1/4)2(1 + (1 − y4 n)1/2)

that gives a better approximation

56

computing π: modular equations

n yn

1 αn

0.41421356237309515 2.9142135623730985 1 0.0037348854633251377 3.1415926 462135446 2 2.4323113418818725E-11 3.14159265358979 4 3 4.3750867904265306E-44 3.14159265358979 4 4 4.579907470729284E-175 3.14159265358979 4 5 0.0 3.14159265358979 4 6 0.0 3.14159265358979 4

57

computing π

• Buffon’s needle

– Throwing n needles gives approximately 1

2 log10 n digits

– To get 1 additional digit throw 100 times as many needles.

• In/circumscribed polygon, and Arctan series

– n iterations gives O(n) digits – To double the number of digits, double the number of iterations

• Modular equations

– n iterations gives approximately 4n digits – To quadruple the number of digits, make 1 additional iteration http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf

58

Toy system R(4, 10, −99, 99)

• Simulator for R(4, 10, −99, 99) available

eksempler/noter/real/ApproxReal.java

59