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 x 0 = 0.671875 x n + 1 = 4x n ( 1 − x n )

Example: Chaos Theory n exact decimal gmpy.mpf float 0 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 1 a 0 + 1 a 1 + 1 a 2 + 1 a 3 + ... a k partial quotients [ a 0 ; a 1 , a 2 , a 3 , . . . , a k ]

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 e 1/n and tan ( 1/n ) – patterns in c.f. e = [ 2 ; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, . . . ] e 1/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 other 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 � � a b c d exy + fx + gy + h = ( x, y ) e f g h

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 � � c + ap d + bp a b change coeffs to g + ep h + fp e f else request q from y and � � b + aq a d + cq c change coeffs to f + eq e h + gq g

Gosper’s algorithm √ π 2 + [ ] [ ] � � 0 1 1 0 0 0 0 1 [ ]

Gosper’s algorithm √ π 2 + [ 3 ] [ ] 3 � � 1 3 0 1 0 1 0 0 [ ]

Gosper’s algorithm √ π 2 + [ 3 ; 7 ] [ ] 7 � � 7 22 1 3 0 7 0 1 [ ]

Gosper’s algorithm √ π 2 + [ 3 ; 7 ] [ 1 ] 1 � � 29 7 4 1 7 0 1 0 [ ]

Gosper’s algorithm √ π 2 + [ 3 ; 7 ] [ 1 ; 2 ] 2 � � 65 29 9 4 14 7 2 1 [ ]

Gosper’s algorithm √ π 2 + [ 3 ; 7 ] [ 1 ; 2 ] � � 14 7 2 1 9 1 1 0 4 [ 4 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7 ] [ 1 ; 2, 2 ] 2 � � 35 14 5 2 19 9 2 1 [ 4 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2 ] 15 � � 530 212 35 14 287 136 19 9 [ 4 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2 ] � � 287 136 19 9 243 76 16 5 1 [ 4 ; 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2 ] � � 243 76 16 5 44 60 3 4 1 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2, 2 ] 2 � � 562 243 37 16 148 44 10 3 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2, 2, 2 ] 2 � � 1367 562 90 37 340 148 23 10 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15 ] [ 1 ; 2, 2, 2, 2, 2 ] 2 � � 3296 13672 217 90 828 340 56 23 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15, 1 ] [ 1 ; 2, 2, 2, 2, 2 ] 1 � � 3513 14572 3296 1367 884 363 828 340 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15, 1 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] 2 � � 8483 35132 7959 3296 2131 884 1996 828 [ 4 ; 1, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15, 1 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] � � 2131 884 1996 828 2090 861 1971 812 3 [ 4 ; 1, 1, 3 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15, 1 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] � � 2090 861 1971 812 41 23 25 16 1 [ 4 ; 1, 1, 3, 1 ]

Gosper’s algorithm √ π 2 + [ 3 ; 7, 15, 1, . . . ] [ 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 2 × [ 1 ] [ ] 1 � � 1 0 1 0 0 1 0 0 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ] [ 1 ] 1 � � 1 1 1 1 1 0 0 0 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ] [ 1 ; 2 ] 2 � � 3 1 3 1 2 1 0 0 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2 ] [ 1 ; 2 ] 2 � � 9 3 3 1 4 2 2 1 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2 ] [ 1 ; 2, 2 ] 2 � � 21 9 7 3 10 4 5 2 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2 ] [ 1 ; 2, 2 ] 2 � � 49 21 21 9 25 10 10 4 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2 ] [ 1 ; 2, 2, 2 ] 2 � � 119 49 51 21 60 25 24 10 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2 ] [ 1 ; 2, 2, 2 ] 2 � � 289 119 119 49 144 60 60 25 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2 ] [ 1 ; 2, 2, 2, 2 ] 2 � � 697 289 287 119 348 144 145 60 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2 ] 2 � � 1681 697 697 289 841 348 348 144 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2, 2 ] 2 � � 4059 1681 1683 697 2030 841 840 348 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2, 2 ] 2 � � 9801 4059 4059 1681 4900 2030 2030 841 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] 2 � � 23661 9801 9799 4059 11830 4900 4901 2030 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] 2 � � 57121 23661 23661 9801 28561 11830 11830 4900 [ ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 2, 2, 2, 2, 2, 2 ] [ 1 ; 2, 2, 2, 2, 2, 2 ] boring, giving up 2 [ 2 ]

Gosper’s algorithm √ √ 2 2 × [ 1 ; 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__