Math errors in python

[Alex Martelli]

Heiko Wundram <heikowu at ceosg.de> wrote:

> Am Sonntag, 19. September 2004 09:39 schrieb Gary Herron:
> > That's called rational arithmetic, and I'm sure you can find a package
> > that implements it for you.  However what would you propose for
> > irrational numbers like sqrt(2) and transcendental numbers like PI?
> 
> Just as an example, try gmpy. Unlimited precision integer and rational
> arithmetic. But don't think that they implement anything more than the four
> basic operations on rationals, because algorithms like sqrt and pow become so
> slow, that nobody sensible would use them, but rather just stick to the
> binary arithmetic the computer uses (although this might have some minor
> effects on precision, but these can be bounded).

Guilty as charged, but with a different explanation.  I don't support
raising a rational to a rational exponent, not because it would "become
slow", but because it could not return a rational result in general.
When it CAN return a rational result, I'm happy as a lark to support it:

>>> x = gmpy.mpq(4,9)
>>> x ** gmpy.mpq(1,2)
mpq(2,3)
>>> 

I.e. raising to the power 1/2 (which is the same as saying, taking the
square root) is supported in gmpy only when the base is a rational which
IS the square of some other rational -- and similarly for other
fractional exponents.

Say you're content with finite precision, and you problem is that
getting only a few dozen bits' worth falls far short of your ambition,
as you want _thousands_.  Well, you don't have to "stick to the
arithmetic your computer uses", with its paltry dozens of bits' worth of
precision -- you can have just as many as you wish.  For example:

For example...:

>>> x=gmpy.mpf(2, 2222)
>>> x
mpf('2.e0',2222)
>>> y=gmpy.fsqrt(x)
>>> y
mpf('1.41421356237309504880168872420969807856967187537694807317667973799
073247846210703885038753432764157273501384623091229702492483605585073721
264412149709993583141322266592750559275579995050115278206057147010955997
160597027453459686201472851741864088919860955232923048430871432145083976
260362799525140798968725339654633180882964062061525835239505474575028775
996172983557522033753185701135437460340849884716038689997069900481503054
402779031645424782306849293691862158057846311159666871301301561856898723
723528850926486124949771542183342042856860601468247207714358548741556570
696776537202264854470158588016207584749226572260020855844665214583988939
4437092659180031138824646815708263e0',2222)
>>> 

Of course, this still has bounded accuracy (gmpy doesn't do constructive
reals...):

>>> x-(y*y)
mpf('1.21406321925474744732602075007044436621136403661789690072865954475
776298522118244419272674806546441529118557492550101271984681584381130555
892259118178248950179953390159664508815540959644741794226362686473376767
055696411211498987561487078708187675060063022704148995680107509652317604
479364576039827518913272446772069713871266672454279184421635785339332972
791970690781583948212784883346298572710476658954707852342842150889381157
563045936231138515406709376167997169879900784347146377935422794796191261
624849740964942283842868779082292557869166024095318326003777296248197487
885858223175591943112711481319695526039760318353849240080721341697065981
8471278600062647147473105883272095e-674',2222)

i.e., there IS an error of about 10 to the minus 674 power, i.e. a
precision of barely more than a couple of thousands of bits -- but then,
that IS what you asked for, with that '2222'!-)

Computing square roots (or whatever) directly on rationals would be no
big deal, were there demand -- you'd still have to tell me what kind of
loss of accuracy you're willing to tolerate, though.  I personally find
it handier to compute with mpf's (high-precision floats) and then turn
the result into rationals with a Stern-Brocot algorithm...:

>>> z=gmpy.f2q(y,-2000)
>>> z
mpq(87787840362947056221389837099888119784184900622573984346903816053706
510038371862119498502008227696594958892073744394524220336403937617412073
521953746033135074321986796669379393248887099312745495535792954890191437
233230436746927180393035328284490481153041398619700720943077149557439382
34750528988254439L,62075377226361968890337286609853165704271551096494666
544033975362265504696870569409265091955693911548812764050925469857560059
623789287922026078165497542602022607603900854658753038808290787475128940
694806084715129308978288523742413573494221901565588452667869917019091383
93125847495825105773132566685269L)

If you need the square root of two as a rational number with an error of
less than 1 in 2**-2000, I think this is a reasonable approach.  As for
speed, this is quite decently usable in an interactive session in my
cheap and cheerful Mac iBook 12" portable (not the latest model, which
is quite a bit faster, much less the "professional" Powerbooks -- I'm
talking about an ageing, though good-quality, consumer-class machine!).

gmpy (or to be more precise the underlying GMP library) runs optimally
on AMD Athlon 32-bit processors, which happen to be dirt cheap these
days, so a cleverly-purchased 300-dollars desktop Linux PC using such an
Athlon chip would no doubt let you use way more than these humble couple
thousand bits for such interactive computations while maintaining a
perfectly acceptable interactive response time.


Alex