Sage Quickstart for Numerical Analysis

This Sage quickstart tutorial was developed for the MAA PREP Workshop “Sage: Using Open-Source Mathematics Software with Undergraduates” (funding provided by NSF DUE 0817071). It is licensed under the Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA).

Sage includes many tools for numerical analysis investigations.

The place to begin is the default decimal number types in Sage.

Basic Analysis

• The RealField class using arbitrary precision (implemented with MPFR).

• The default real numbers (RR) is RealField(53) (i.e., 53 bits of precision).

• But you can make them live at whatever precision you wish.

Sage

sage: ring=RealField(3)

Python

>>> from sage.all import *
>>> ring=RealField(Integer(3))

To print the actual number (without rounding off the last few imprecise digits to only display correct digits), call the .str() method:

Sage

sage: print(ring('1').nextabove())
1.2

Python

>>> from sage.all import *
>>> print(ring('1').nextabove())
1.2

Sage

sage: print(ring('1').nextabove().str())
1.2
sage: print(ring('1').nextbelow().str())
0.88

Python

>>> from sage.all import *
>>> print(ring('1').nextabove().str())
1.2
>>> print(ring('1').nextbelow().str())
0.88

Let’s change our precision.

Sage

sage: ring=RealField(20)
sage: print(ring('1').nextabove().str())
1.0000019
sage: print(ring('1').nextbelow().str())
0.99999905

Python

>>> from sage.all import *
>>> ring=RealField(Integer(20))
>>> print(ring('1').nextabove().str())
1.0000019
>>> print(ring('1').nextbelow().str())
0.99999905

You can also specify the rounding mode.

Sage

sage: ringup=RealField(3,rnd='RNDU')
sage: ringdown=RealField(3,rnd='RNDD')

Python

>>> from sage.all import *
>>> ringup=RealField(Integer(3),rnd='RNDU')
>>> ringdown=RealField(Integer(3),rnd='RNDD')

Sage

sage: ring(1/9).str()
'0.11111116'

Python

>>> from sage.all import *
>>> ring(Integer(1)/Integer(9)).str()
'0.11111116'

Sage

sage: ringup(1/9).str()
'0.13'

Python

>>> from sage.all import *
>>> ringup(Integer(1)/Integer(9)).str()
'0.13'

Sage

sage: ringdown(1/9).str()
'0.10'

Python

>>> from sage.all import *
>>> ringdown(Integer(1)/Integer(9)).str()
'0.10'

Sage

sage: ring(1/9).str(base=2)
'0.00011100011100011100100'

Python

>>> from sage.all import *
>>> ring(Integer(1)/Integer(9)).str(base=Integer(2))
'0.00011100011100011100100'

Let’s see exactly what the fraction is.

Sage

sage: ring(1/9).exact_rational()
233017/2097152

Python

>>> from sage.all import *
>>> ring(Integer(1)/Integer(9)).exact_rational()
233017/2097152

Sage

sage: n(233017/2097152)
0.111111164093018

Python

>>> from sage.all import *
>>> n(Integer(233017)/Integer(2097152))
0.111111164093018

Converting to floating point binary (IEEE format)

Sage

sage: x=13.28125

Python

>>> from sage.all import *
>>> x=RealNumber('13.28125')

Here is x as a binary floating point number.

Sage

sage: x.str(base=2)
'1101.0100100000000000000000000000000000000000000000000'

Python

>>> from sage.all import *
>>> x.str(base=Integer(2))
'1101.0100100000000000000000000000000000000000000000000'

From this binary floating point representation, we can construct x again.

Sage

sage: xbin=2^3 +2^2 + 2^0+2^-2+2^-5; xbin
425/32

Python

>>> from sage.all import *
>>> xbin=Integer(2)**Integer(3) +Integer(2)**Integer(2) + Integer(2)**Integer(0)+Integer(2)**-Integer(2)+Integer(2)**-Integer(5); xbin
425/32

Sage

sage: xbin.n()
13.2812500000000

Python

>>> from sage.all import *
>>> xbin.n()
13.2812500000000

To check, let’s ask Sage what x is as an exact fraction (no rounding involved).

Sage

sage: x.exact_rational()
425/32

Python

>>> from sage.all import *
>>> x.exact_rational()
425/32

Now let’s convert x to the IEEE 754 binary format that is commonly used in computers. For IEEE 754, the first step in getting the binary format is to normalize the number, or express the number as a number between 1 and 2 multiplied by a power of 2. For our x above, we multiply by \(2^{-3}\) to get a number between 1 and 2.

Sage

sage: (x*2^(-3)).str(base=2)
'1.1010100100000000000000000000000000000000000000000000'

Python

>>> from sage.all import *
>>> (x*Integer(2)**(-Integer(3))).str(base=Integer(2))
'1.1010100100000000000000000000000000000000000000000000'

The fraction part of the normalized number is called the mantissa (or significand ).

Sage

sage: f=(x*2^(-3)).frac() # .frac() gets the fraction part, the part after the decimal
sage: mantissa=f.str(base=2)
sage: mantissa
'0.10101001000000000000000000000000000000000000000000000'

Python

>>> from sage.all import *
>>> f=(x*Integer(2)**(-Integer(3))).frac() # .frac() gets the fraction part, the part after the decimal
>>> mantissa=f.str(base=Integer(2))
>>> mantissa
'0.10101001000000000000000000000000000000000000000000000'

Since we multiplied by \(2^{-3}\) to normalize the number, we need to store this information. We add 1023 in order to not have to worry about storing negative numbers. In the below, c will be our exponent.

Sage

sage: c=(3+1023).str(base=2)
sage: c
'10000000010'

Python

>>> from sage.all import *
>>> c=(Integer(3)+Integer(1023)).str(base=Integer(2))
>>> c
'10000000010'

Note that c has 11 bits, which is exactly what we want.

Sage

sage: len(c)
11

Python

>>> from sage.all import *
>>> len(c)
11

Evaluating mantissa[2:54] will give the first 52 binary digits after the decimal point of the mantissa. Note that we don’t need to store the leading 1 before the decimal point because it will always be there from the way we normalized things. This lets us get 53-bit precision using only 52 bits of storage.

Sage

sage: len(mantissa[2:54])
52

Python

>>> from sage.all import *
>>> len(mantissa[Integer(2):Integer(54)])
52

Since the original number was positive, our sign bit is zero.

Sage

sage: sign='0'

Python

>>> from sage.all import *
>>> sign='0'

So here is our 64-bit double-precision floating point number.

Sage

sage: sign+' '+c+' '+mantissa[2:54] # the [2:] just chops off the '0.', since we just need to store the digits after the decimal point
'0 10000000010 1010100100000000000000000000000000000000000000000000'

Python

>>> from sage.all import *
>>> sign+' '+c+' '+mantissa[Integer(2):Integer(54)] # the [2:] just chops off the '0.', since we just need to store the digits after the decimal point
'0 10000000010 1010100100000000000000000000000000000000000000000000'

Sage

sage: len(sign+c+mantissa[2:54]) # it's 64 bits!
64

Python

>>> from sage.all import *
>>> len(sign+c+mantissa[Integer(2):Integer(54)]) # it's 64 bits!
64

Here we convert back to our original number from the floating point representation that we constructed.

Sage

sage: ((-1)^(int(sign)) * 2^(int(c,base=2)-1023)*(1+RR(mantissa[:54], base=2)))
13.2812500000000

Python

>>> from sage.all import *
>>> ((-Integer(1))**(int(sign)) * Integer(2)**(int(c,base=Integer(2))-Integer(1023))*(Integer(1)+RR(mantissa[:Integer(54)], base=Integer(2))))
13.2812500000000

Sage

sage: x
13.2812500000000

Python

>>> from sage.all import *
>>> x
13.2812500000000

So they agree!

Sage uses a cutting-edge numerical library, MPFR, to carry out precise floating point arithmetic using any precision a user specifies. MPFR has a slightly different convention for normalization. In MPFR, we normalize by multiplying by an appropriate power of 2 to make the mantissa an integer, instead of a binary fraction. This allows us to use big integer libraries and sophisticated techniques to carry out calculations at an arbitrary precision.

Sage

sage: x.sign_mantissa_exponent()
(1, 7476679068876800, -49)

Python

>>> from sage.all import *
>>> x.sign_mantissa_exponent()
(1, 7476679068876800, -49)

Sage

sage: 7476679068876800*2^(-49)
425/32

Python

>>> from sage.all import *
>>> Integer(7476679068876800)*Integer(2)**(-Integer(49))
425/32

Note that the mantissa here has the same zero/nonzero bits as the mantissa above (before we chopped off the leading 1 above).

Sage

sage: 7476679068876800.str(base=2)
'11010100100000000000000000000000000000000000000000000'

Python

>>> from sage.all import *
>>> Integer(7476679068876800).str(base=Integer(2))
'11010100100000000000000000000000000000000000000000000'

Interval Arithmetic

Sage also lets you compute using intervals to keep track of error bounds. These basically use the round up and round down features shown above.

Sage

sage: ring=RealIntervalField(10)
sage: a=ring(1/9)
sage: a
0.112?

Python

>>> from sage.all import *
>>> ring=RealIntervalField(Integer(10))
>>> a=ring(Integer(1)/Integer(9))
>>> a
0.112?

The question mark notation means that the number is contained in the interval found by incrementing and decrementing the last digit of the number. See the documentation for real interval fields for details. In the above case, Sage is saying that 1/9 is somewhere between 0.111 and 0.113. Below, we see that 1/a is somewhere between 8.9 and 9.1.

Sage

sage: 1/a
9.0?

Python

>>> from sage.all import *
>>> Integer(1)/a
9.0?

We can get a more precise estimate of the interval if we explicitly print out the interval.

Sage

sage: print((1/a).str(style='brackets'))
[8.9843 .. 9.0157]

Python

>>> from sage.all import *
>>> print((Integer(1)/a).str(style='brackets'))
[8.9843 .. 9.0157]

Included Software

Scipy (included in Sage) has a lot of numerical algorithms. See the Scipy docs.

Mpmath is also included in Sage, and contains a huge amount of numerical stuff. See the mpmath codebase.

The Decimal python module has also been useful for textbook exercises which involved rounding in base 10.

Plotting with precision

Sometimes plotting involves some rather bad rounding errors because plotting calculations are done with machine-precision floating point numbers.

Sage

sage: f(x)=x^2*(sqrt(x^4+16)-x^2)
sage: plot(f,(x,0,2e4))
Graphics object consisting of 1 graphics primitive

Python

>>> from sage.all import *
>>> __tmp__=var("x"); f = symbolic_expression(x**Integer(2)*(sqrt(x**Integer(4)+Integer(16))-x**Integer(2))).function(x)
>>> plot(f,(x,Integer(0),RealNumber('2e4')))
Graphics object consisting of 1 graphics primitive

We can instead make a function that specifically evaluates all intermediate steps to 100 bits of precision using the fast_callable system.

Sage

sage: R=RealField(100) # 100 bits
sage: g=fast_callable(f, vars=[x], domain=R)
sage: plot(g,(x,0,2e4))
Graphics object consisting of 1 graphics primitive

Python

>>> from sage.all import *
>>> R=RealField(Integer(100)) # 100 bits
>>> g=fast_callable(f, vars=[x], domain=R)
>>> plot(g,(x,Integer(0),RealNumber('2e4')))
Graphics object consisting of 1 graphics primitive