MPFR IN SAGE

• What is MPFR?

  • A high performance C library implementing high precision real numbers, based on GMP
  • Where 'high' means 'arbitrary', limited only by machine memory
  • Developed at LORIA (France)
  • Free, GPL'd, very high performance, widely used, well tested, under active development and so on. (Performance comparison at http://www.medicis.polytechnique.fr/~pphd/mpfr/timings-220.html)
  
  Watch what happens when we do this:

  kyle-schalms-Computer:~/downloads/sage-0.10.12 kyle$ ./sage
  --------------------------------------------------------
  | SAGE Version 0.10.12, Build Date: 2006-01-30-0901    |
  | Distributed under the GNU General Public License V2  |
  | For help type <object>?, <object>??, %magic, or help |
  --------------------------------------------------------
  
  
  sage: 1.111111111111111111111111111111111111111111
   _1 = 1.1111111111111109
  sage: 1.111111111111111111111111111111111111111111.prec()
   _2 = 53
  
  

  Surprised? Stay tuned!

• MPFR in SAGE - History

  • Before (circa Aug 2005), SAGE real numbers were Python floats. Fixed precision. (53 bit?)
  • I needed more!
  • William and I implemented a Pyrex wrapper for MPFR.
  • I was happy.

• Demo

  sage: type(1.1)
   _7 = <type 'mpfr.RealNumber'>
  
  sage: R200=RealField(200)
  sage: R200(1.111111111111111111111111111111111111111111)
  _21 = 1.1111111111111109384097517249756492674350738525390625000000000
  

  Still not what we intended! Try again.

  sage: R200('1.111111111111111111111111111111111111111111')
  _22 = 1.1111111111111111111111111111111111111111109999999999999999990
  

  Much better! Now let's do some real stuff.

  sage: R200(pi)
  _25 = 3.1415926535897932384626433832795028841971693993751058209749445
  sage: R200.pi()
  _27 = 3.1415926535897932384626433832795028841971693993751058209749445
  

  200 digits of pi (both do same thing, but in different ways).

  sage: R200(e)
  _26 = 2.7182818284590452353602874713526624977572470936999595749669679
  

  As much e as you'll ever want.

  sage: RealField(10000)(2.0)
  _27 = (many screenfuls)
  

  What if you use numbers of different precision together?

  sage: One=RealField(60)(1)
  sage: One
  _56 = 1.0000000000000000000
  sage: one=RealField(30)(1)
  _57 = 1.0000000000
  sage: One+one
  _60 = 2.0000000000
  sage: (One+one).prec()
  _61 = 30
  

  The result is in the lower precision.
  
  
  Functions included in MPFR/SAGE:
  • Basic arithmetic

    sage: two=R200(2)
    sage: two+two
    _47 = 4.0000000000000000000000000000000000000000000000000000000000000
    sage: two^(1/2)
    _46 = 1.4142135623730950488016887242096980785696718753769480731766796
    

  • Log, exp

  sage: log(two)
  _65 = 0.69314718055994530941723212145817656807550013436025525412067998
  
  sage: log(R200(e))
  _36 = 1.0000000000000000000000000000000000000000000000000000000000000
  

  Hurrah! Suppose I want to know how many digits are in 200 bits.

  sage: 200*2.0.log10()
  _49 = 60.205999132796236
  

  About 60. What's the log of -1?

  sage: R200(-1).log()
  _37 = NaN
  

  Not a number (as of yet doesn't return complex number). A bit cleaner than

  sage: log(-1)
  
    ***   impossible assignment t_COMPLEX --> t_REAL.
  ------------------------------------------------------------
  Traceback (most recent call last):
    File "", line 1, in ?
    File 
  "/Users/kyle/downloads/sage-0.10.12/local/lib/python2.4/site-packages/sage/misc/functional.py", 
  line 528, in log
      return float(pari(x).log())
  RuntimeError
  

  Some different behaviours:

  sage: log(1)
  _52 = 0.0
  sage: log(1.0)
  _53 = 0.00000000000000000
  sage: log(R200(1.0))
  _54 = 0.00000000000000000000000000000000000000000000000000000000000000
  

  • Trig and inverse trig

  sage: sin(R200(pi/3))
  _64 = 0.86602540378443864676372317075293618347140262690519031402790376
  

  You can also use 'object-oriented' notation:

  sage: R200(pi).sin()
  _40 = 0.00000000000000000000000000000000000000000000000000000000000011419936994248698005287480012038354990556371286030641123282636
  

  printed in fixed point notation, with error

  sage: _.log2()
  _41 = -202.44605909700240181910412473189164935339689258145121229561184
  

  of less than 2^(-200). Good!
  • Hyperbolic (sinh etc)
  • Gamma, zeta

    sage: two.zeta() - R200(pi^2/6)
    _45 = 0.00000000000000000000000000000000000000000000000000000000000000
    

  • And more! dir(<number>) or <number>.<tab> for complete list.

• Pyrex Code

  • Pyrex code for MPFR wrapper found in devel/sage-0.10.12/sage/ext/mpfr.pyx .
  • Sample Pyrex code for square root - a very rough outline:

        ctypedef struct __mpfr_struct:
            pass
        ctypedef __mpfr_struct* mpfr_t
    
        ... more header stuff ...
    
    cdef class RealNumber(element.RingElement):
    
        ....
        cdef mpfr_t value
        ....
    

    Very approximate implementation of constructor:

        def __init__(self, RealField parent, x=0, int base=10, special=None):
            ....
            mpfr_init2(self.value, parent.prec)
            s = str(x)
            mpfr_set_str(self.value, s, base, GMP_RNDZ):
            ....
    

    Actual implementation of square root function:

        def sqrt(self):
            """
            Return the square root of self.
            """
    

    This is the help string. It's what you see when you type help(1.0.sqrt).

            cdef RealNumber x
    

    This declaration is just like a C typedef. It sets no value.

            x = RealNumber(self._parent, None)
    

    Assignment - calls the init method shown above.

            mpfr_sqrt(x.value, self.value, self._parent.rnd)
    

    Where the work happens. self is the input and x is the output.

            return x
    

    All done.

    sage: sqrt(two)
    _76 = 1.4142135623730950488016887242096980785696718753769480731766796
    sage: sqrt(-two)
    _77 = NaN
    

  • To compile:
    • pyrexc mpfr.pyx produces mpfr.c. The part after that is mysterious.
    • better: sage -p recompiles all changed pyrex modules.
  • Pyrex overview at http://www.cosc.canterbury.ac.nz/~greg/python/Pyrex/

• The Future

  • Interval arithmetic with MPFI would be wonderful.
  • Not only is your computation high precision, but you know how high.
  • No wasted bits; don't put 10,000 bits in and hope for 100 good ones out.