[Tutor] more precision/predictability in calculations [gmpy module]

[Daniel Yoo]

On Fri, 13 Sep 2002, Prahlad Vaidyanathan wrote:

> def diff2 (fn, var, x0, y0, epsilon=0.0001) :
>     """ Given a 2 variable function, differentiates it partially
>     w.r.t 'var', at the point (x0, y0) """
>     if var == "x" :
>         return (fn(x0+epsilon, y0) - fn(x0, y0))/epsilon
>     elif var == "y" :
>         return (fn(x0, y0+epsilon) - fn(x0, y0))/epsilon
>
> def satisfies_cr_eqns (u, v, z0) :
>     """Checks if the 2 real-valued functions 'u' and 'v' satisfy
>     the Cauchy-Reimann equations at z0"""
>     ux = diff2(u, "x", z0.real, z0.imag) ##; print ux
>     uy = diff2(u, "y", z0.real, z0.imag) ##; print uy
>     vx = diff2(v, "x", z0.real, z0.imag) ##; print vx
>     vy = diff2(v, "y", z0.real, z0.imag) ##; print vy
>     ## return ux == vy and uy == -vx
>     return (ux - vy) < pow(10,-3) and (uy + vx) < pow(10, -3)
>
> """
>
>     Now, I construct 2 functions like so :
>
> """
>
> def f1(a, b) :
>     return pow(math.e, a) * math.cos(b)
>
> def f2(a, b) :
>     return pow(math.e, a) * math.sin(b)
>
> """
>
>     These 2 functions do satisfy the c-r equations for any point in the
> complex plane.

?!?  Waaaa....  *grin*


> But, when I run satisfy_cr_eqns(), it returns true for some points and
> false for others. I presumed that the higher the value of the complex
> number, the lower the precision would be. But, the function turns out to
> be true for complex(99,99), and false for complex(3,3) !!
>
>     I know this may be asking too much of Python, but is there any way I
> can make these calculations more precise ? or, at least more predictable
> ? :-)


Python uses floating point to represent these numbers, so that's where
we're losing precision.  We can probably do a bit better by using the
'gmp' General MultiPrecision library --- it's made to handle hardcore
arbitrary precision arithmetic:

    http://www.swox.com/gmp/


There's are a few Python wrappers around portions of the GMP library, and
we can find a good one in the 'gmpy' module:

    http://gmpy.sourceforge.net/


The gmpy module provides a function called 'mpf()', which constructs
multi-precision floats.  The documentation is a little terse, so here's a
small example:

###
>>> import gmpy
>>> mpf = gmpy.mpf
>>> mpf('3.1415926')
mpf('3.1415926e0')
>>> pi = mpf('3.1415926')
>>> pi / mpf(10**100000)
mpf('3.1415926e-100000')
###

The last example is somewhat neat, considering that it doesn't work (at
least, not easily) in standard Python:

###
>>> 3.1415926 / 10**100000
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
OverflowError: long int too large to convert to float
###


The gmpy module should come with a test file called 'gmpy_test_08_mpf.py'
that showcases the sort of precision we can expect from these numbers.




Best of wishes to you!