Marcin Ciura Google Poland Krakw, May 12, 2011 Reals in Python - - PowerPoint PPT Presentation

Exact Real Arithmetic in Python

Marcin Ciura

Google Poland Kraków, May 12, 2011

Reals in Python

float, math

◮ IEEE 754 double arithmetic ◮ built-in type, standard library

Reals in Python

decimal.Decimal

◮ decimal floating-point arithmetic ◮ since Python 2.4

Reals in Python

gmpy.mpf

◮ multiple-precision arithmetic ◮ third-party library

Example: Chaos Theory

logistic map x0 = 0.671875 xn+1 = 4xn(1 − xn)

Example: Chaos Theory

n exact decimal gmpy.mpf float 0.671875 0.671875 0.671875 0.671875 1 0.881835 0.881835 0.881835 0.881835 2 0.416805 0.416805 0.416805 0.416805 3 0.972314 0.972314 0.972314 0.972314 4 0.107675 0.107675 0.107675 0.107675

Example: Chaos Theory

n exact decimal gmpy.mpf float 53 0.717980 0.717980 0.717980 0.749943 54 0.809938 0.809938 0.809937 0.750112 55 0.615751 0.615751 0.615755 0.749775 56 0.946406 0.946406 0.946403 0.750449 57 0.202886 0.202886 0.202897 0.749100

Example: Chaos Theory

n exact decimal gmpy.mpf float 68 0.910298 0.910298 0.893725 0.000935 69 0.326622 0.326622 0.379919 0.003740 70 0.879760 0.879760 0.942322 0.014904 71 0.423127 0.423127 0.217403 0.058727 72 0.976362 0.976362 0.680555 0.221116

Example: Chaos Theory

n exact decimal gmpy.mpf float 91 0.035427 0.023866 0.241684 0.487988 92 0.136690 0.093186 0.733092 0.999422 93 0.472024 0.338009 0.782670 0.002306 94 0.996869 0.895036 0.680388 0.009206 95 0.012483 0.375783 0.869839 0.036486

Example: Chaos Theory corollary

for some applications, multiple-precision arithmetic is not good enough

we must go deeper

Reals in Python

◮ float ◮ decimal.Decimal ◮ gmpy.mpf

round-off errors (begot numerical analysis) fast

Exact Reals

◮ continued fractions ◮ scaled integers ◮ Möbius transformations ◮ other methods

no round-off errors slowish

Outline

So far

◮ Introduction

Now

◮ Continued fractions ◮ Gosper’s algorithm ◮ Python implementation ◮ Live demo

Continued Fractions

a0 + 1 a1 + 1 a2 + 1 a3 + 1 ... ak partial quotients [a0; a1, a2, a3, . . . , ak]

Continued Fractions

byproduct of Euclid’s gcd algorithm 96 = 1 × 65 + 31 65 = 2 × 31 + 3 31 = 10 × 3 + 1 3 = 3 × 1 96 65 = [1; 2, 10, 3]

Continued Fractions

◮ 1 : 1 correspondence with reals ◮ finite c.f. – rationals ◮ infinite c.f. – irrationals ◮ truncated c.f. – best approximations ◮ usually ∼ 0.97027 partial quotients per digit

(e.g. 971 p.q. for initial 1000 digits of π)

Continued Fractions

quadratic irrationals – periodic c.f. √ 2 = [1; 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, . . .] − √ 2 = [−2; 1, 1, 2, 2, 2, 2, 2, 2, 2, . . .] √ 3 = [1; 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, . . .]

Continued Fractions

e1/n and tan(1/n) – patterns in c.f. e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, . . .] e1/2 = [1; 1, 1, 1, 5, 1, 1, 9, 1, 1, 13, . . .] tan(1/2) = [0; 1, 1, 4, 1, 8, 1, 12, 1, 16, . . .]

Continued Fractions

• ther irrationals – irregular c.f.

π = [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, . . .]

3

√ 2 = [1; 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, . . .]

Continued Fractions

◮ Gosper, 1972: +, −, ×, / ◮ Newton, 1671:

√

◮ Brouncker, 1656: π ◮ Ciura, 2007: tan, arctan, exp, log

Gosper’s algorithm

bihomographic function z(x, y) = axy + bx + cy + d exy + fx + gy + h =

• a b c d

e f g h

• (x, y)

Gosper’s algorithm

special cases of bihomographic function x + y = 0xy + 1x + 1y + 0 0xy + 0x + 0y + 1 x − y = 0xy + 1x − 1y + 0 0xy + 0x + 0y + 1 xy = 1xy + 0x + 0y + 0 0xy + 0x + 0y + 1 x/y = 0xy + 1x + 0y + 0 0xy + 0x + 1y + 0

Gosper’s algorithm

lazy evaluation of bihomographic function if [redacted] then emit r and change coeffs to

• e

f g h a − er b − fr c − gr d − hr

• else if [redacted] then request p from x and

change coeffs to

• c + ap d + bp a b

g + ep h + fp e f

• else request q from y and

change coeffs to

• b + aq a d + cq c

f + eq e h + gq g

Gosper’s algorithm

+ π [ ] √ 2 [ ]

• 0 1 1 0

0 0 0 1

• [ ]

Gosper’s algorithm

+ π [3] √ 2 [ ]

• 1 3 0 1

0 1 0 0

• [ ]

3

Gosper’s algorithm

+ π [3; 7] √ 2 [ ]

• 7 22 1 3

7 0 1

• [ ]

7

Gosper’s algorithm

+ π [3; 7] √ 2 [1]

• 29 7 4 1

7 0 1 0

• [ ]

1

Gosper’s algorithm

+ π [3; 7] √ 2 [1; 2]

• 65 29 9 4

14 7 2 1

• [ ]

2

Gosper’s algorithm

+ π [3; 7] √ 2 [1; 2]

• 14 7 2 1

9 1 1 0

• [4]

4

Gosper’s algorithm

+ π [3; 7] √ 2 [1; 2, 2]

• 35 14 5 2

19 9 2 1

• [4]

2

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2]

• 530 212 35 14

287 136 19 9

• [4]

15

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2]

• 287 136 19 9

243 76 16 5

• [4; 1]

1

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2]

• 243 76 16 5

44 60 3 4

• [4; 1, 1]

1

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2, 2]

• 562 243 37 16

148 44 10 3

• [4; 1, 1]

2

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2, 2, 2]

• 1367 562 90 37

340 148 23 10

• [4; 1, 1]

2

Gosper’s algorithm

+ π [3; 7, 15] √ 2 [1; 2, 2, 2, 2, 2]

• 3296 13672 217 90

828 340 56 23

• [4; 1, 1]

2

Gosper’s algorithm

+ π [3; 7, 15, 1] √ 2 [1; 2, 2, 2, 2, 2]

• 3513 14572 3296 1367

884 363 828 340

• [4; 1, 1]

1

Gosper’s algorithm

+ π [3; 7, 15, 1] √ 2 [1; 2, 2, 2, 2, 2, 2]

• 8483 35132 7959 3296

2131 884 1996 828

• [4; 1, 1]

2

Gosper’s algorithm

+ π [3; 7, 15, 1] √ 2 [1; 2, 2, 2, 2, 2, 2]

• 2131 884 1996 828

2090 861 1971 812

• [4; 1, 1, 3]

3

Gosper’s algorithm

+ π [3; 7, 15, 1] √ 2 [1; 2, 2, 2, 2, 2, 2]

• 2090 861 1971 812

41 23 25 16

• [4; 1, 1, 3, 1]

1

Gosper’s algorithm

+ π [3; 7, 15, 1, . . .] √ 2 [1; 2, 2, 2, 2, 2, 2, . . .]

ad nauseam

[4; 1, 1, 3, 1, . . .]

Gosper’s algorithm

Gosper’s algorithm

× √ 2 [ ] √ 2 [ ]

• 1 0 0 0

0 0 0 1

• [ ]

Gosper’s algorithm

× √ 2 [1] √ 2 [ ]

• 1 0 1 0

0 1 0 0

• [ ]

1

Gosper’s algorithm

× √ 2 [1] √ 2 [1]

• 1 1 1 1

1 0 0 0

• [ ]

1

Gosper’s algorithm

× √ 2 [1] √ 2 [1; 2]

• 3 1 3 1

2 1 0 0

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2] √ 2 [1; 2]

• 9 3 3 1

4 2 2 1

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2] √ 2 [1; 2, 2]

• 21 9 7 3

10 4 5 2

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2] √ 2 [1; 2, 2]

• 49 21 21 9

25 10 10 4

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2] √ 2 [1; 2, 2, 2]

• 119 49 51 21

60 25 24 10

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2] √ 2 [1; 2, 2, 2]

• 289 119 119 49

144 60 60 25

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2] √ 2 [1; 2, 2, 2, 2]

• 697 289 287 119

348 144 145 60

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2]

• 1681 697 697 289

841 348 348 144

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2]

• 4059 1681 1683 697

2030 841 840 348

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2]

• 9801 4059 4059 1681

4900 2030 2030 841

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2, 2]

• 23661 9801 9799 4059

11830 4900 4901 2030

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2, 2]

• 57121 23661 23661 9801

28561 11830 11830 4900

• [ ]

2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2, 2]

boring, giving up

[2] 2

Gosper’s algorithm

× √ 2 [1; 2, 2, 2, 2, 2, 2] √ 2 [1; 2, 2, 2, 2, 2, 2]

boring, giving up

[2; ∞] ∞

Gosper’s algorithm

◮ irrational arguments, rational result –

undecidable

◮ give up after 100 resultless iterations ◮ risk of incorrect result at most

1.5 × 10−31 × number of partial quotients

What’s needed

What’s needed

◮ integer bignums ◮ lazy evaluation ◮ garbage collection

Python has got it all

◮ long type ◮ generators ◮ garbage collection

Implementation

◮ source generators: c.f. for rational numbers,

canned sequences of p.q., and π

◮ intermediate generators: bihomographic

function, homographic function, exp, log, tan, atan, and sqrt

◮ sinks: type conversions, __cmp__, __str__,

and __repr__

Implementation

◮ functions and methods return (chained)

generators

• 1 0

0 2

• tan
• 0 2 0 0

1 0 0 1

• sin

◮ calling a sink makes it request p.q. from its

inputs, they request p.q. from theirs, etc.

◮ generated partial quotients are cached

Implementation

◮ cf.cf ◮ http://tnij.org/continued-fractions/ ◮ 2000 lines of Python ◮ function parity with math

Like a Computer Algebra System

>>> (sqrt(2) - 1) * (sqrt(2) + 1) 1.

Like a Computer Algebra System

>>> sqrt(5 + 2*sqrt(6)) == sqrt(2) + sqrt(3) True

Like a Computer Algebra System

>>> [sin(x)**2 + cos(x)**2 for x in ... [pi * random.random(), ... pi * random.random(), ... pi * random.random()]] [1., 1., 1.]

Like a Computer Algebra System

>>> [cos(3*x)/cos(x) == ... cos(x)**2 - 3*sin(x)**2 for x in ... [pi * random.random(), ... pi * random.random(), ... pi * random.random()]] [True, True, True]

Like a Computer Algebra System

>>> [x == x + cf(10)**-1000 for x in ... [cf(random.random()), ... cf(random.random()), ... cf(random.random())]] [False, False, False]

. . . But Not Quite

>>> sqrt(2) == sqrt(2) + cf(10)**-1000 True

Possible Applications

◮ chaos theory ◮ computational geometry ◮ computational number theory ◮ educational applications

Live Demo