# 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: ring=RealField(3)
```

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

method with
the option `truncate=False`

.

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

```
sage: print(ring('1').nextabove().str(truncate=False))
1.2
sage: print(ring('1').nextbelow().str(truncate=False))
0.88
```

Let’s change our precision.

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

You can also specify the rounding mode.

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

```
sage: ring(1/9).str(truncate=False)
'0.11111116'
```

```
sage: ringup(1/9).str(truncate=False)
'0.13'
```

```
sage: ringdown(1/9).str(truncate=False)
'0.10'
```

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

Let’s see exactly what the fraction is.

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

```
sage: n(233017/2097152)
0.111111164093018
```

## Converting to floating point binary (IEEE format)¶

```
sage: x=13.28125
```

Here is `x`

as a binary floating point number.

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

From this binary floating point representation, we can construct `x`

again.

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

```
sage: xbin.n()
13.2812500000000
```

To check, let’s ask Sage what `x`

is as an exact fraction (no rounding
involved).

```
sage: 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: (x*2^(-3)).str(base=2)
'1.1010100100000000000000000000000000000000000000000000'
```

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

```
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'
```

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: c=(3+1023).str(base=2)
sage: c
'10000000010'
```

Note that `c`

has 11 bits, which is exactly what we want.

```
sage: 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: len(mantissa[2:54])
52
```

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

```
sage: sign='0'
```

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

```
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'
```

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

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

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

```
sage: 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: x.sign_mantissa_exponent()
(1, 7476679068876800, -49)
```

```
sage: 7476679068876800*2^(-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: 7476679068876800.str(base=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: ring=RealIntervalField(10)
sage: a=ring(1/9)
sage: 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: 1/a
9.0?
```

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

```
sage: print((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: f(x)=x^2*(sqrt(x^4+16)-x^2)
sage: plot(f,(x,0,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: 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
```