Lecture 4: Exact and Floating-Point Numbers
===========================================
The capability to compute exactly and in high precision
is one of the main features of computer algebra.
In this lecture we illustrate how to compute exactly
with integer and rational numbers,
and to compute with multiprecision arithmetic.
We can approximate a transcendental number such as \ :math:`\pi`
as a floating-point number up to any number of decimal places.
Instead of a floating-point approximation,
we may approximate with a continued fraction expansion.
In Sage we can compute continued fractions up to a certain number of
terms or up to a given number of bits in the precision.
The convergents of a continued fraction gives us a sequence of
increasingly more accurate rational approximations for \ :math:`\pi`.
We cover machine precision. Irrational and algebraic numbers are defined.
We introduce the reverse operation to computing an approximation
of an irrational number. Given an approximation for a number
and a tolerance, we may find the smallest polynomial that has this
number as a root. For example, with sympy we compute \ :math:`\sqrt{2}`
from an approximation for the \ :math:`\sqrt{2}`,
as \ :math:`\sqrt{2}` is a solution of \ :math:`x^2 - 2`.
Integer and Rational Numbers
----------------------------
The division of two integer numbers gives a :index:`rational number`.
Selecting its numerator and denominator is straightforward:
::
a = 3/4
print type(a)
print numerator(a), denominator(a)
Any irrational or transcendental can be approximated
with continued fractions.
A :index:`continued fraction` is defined by a list of natural numbers,
which are called :index:`convergents`. For example:
::
x = sqrt(2)
c = continued_fraction(x)
print c
show(c)
which prints
``[1; 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...]``
and the output of ``show(c)`` is displayed
in :numref:`figcontinuedfraction`.
.. _figcontinuedfraction:
.. figure:: ./figcontinuedfraction.png
:align: center
The continued fraction expansion of :math:`\sqrt{2}`.
We can verify that
.. math::
1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2}}}}
approximates \ :math:`\sqrt{2}` up to 3 correct decimal places:
::
f = 1 + 1/(2 + 1/(2 + 1/(2 + 1/2)))
print(f)
print(float(f), float(x))
The continued fraction defines *convergents* which are the
consecutive best rational approximations.
::
c.convergents()
The output is a ``lazy list`` which shows only the first three numbers.
To compute the first five convergents, one takes a slice of length 5
and then converts that slice to a list:
::
v = c.convergents()
v[:5].list()
Each convergent is a :index:`rational approximation`.
The output of the last command is
``[1, 3/2, 7/5, 17/12, 41/29]`` which is a list of
increasingly more accurate approximations converging
to \ :math:`\sqrt{2}.`
The quickest way to get a rational approximation for the square root of 2,
accurate to 2 decimal places, goes as follows
::
QQ(sqrt(2).n(digits=2))
With a list comprehension, we can then compute the sequence
of increasingly more accurate rational approximations for ``sqrt(2)``.
::
[QQ(sqrt(2).n(digits=k)) for k in range(2,10)]
Observe that this sequence is different from the convergents,
because the list above contains rational approximations with
a prescribed accuracy (of ``k`` decimal places), whereas the
convergents give the sequence of the next best rational approximations.
Floating-Point Numbers
----------------------
The default floating-point type is ``float`` as in Python.
The :index:`machine precision` is defined as the *smallest*
positive number we can add
to 1.0 and see a difference from 1.0.
A :index:`double float` has 53 bits in its fraction,
a :index:`single float` has a fraction of 22 bits long.
We can verify the machine precision calculation with the predefined
constants of the numpy package.
::
print('machine precision (double float) :', float(2^(-52)))
print('machine precision (single float) :', float(2^(-23)))
import numpy
from numpy import finfo
print(finfo(float).eps)
print(finfo(numpy.float32).eps)
To compute with higher than double precision, we can make
a :index:`RealField`.
Since the argument is in bits, we multiply with the number of decimal places
we want in our precision with the 2-logarithm of 10.
For 50 decimals, we do
::
nbits = 50*log(10,2)
print('number of bits for 50 decimal places :', int(nbits))
R = RealField(nbits)
print('50-digit approximation of sqrt(2) =', R(x))
We will use this field to verify the accuracy of the convergents
calculating in R with a :index:`list comprehension`
::
c = continued_fraction(sqrt(2))
v = c.convergents()
a = v[:20].list()
print(['%.2e' % R(sqrt(2)-p) for p in a])
The numbers are displayed in scientific format and we see the
exponent gradually decreasing till -15, which indicates that the
last convergent is accurate up to 15 decimal places.
A number as \ :math:`\sqrt{2}` is an :index:`irrational number`.
But it is not a :index:`transcendental number`
because \ :math:`\sqrt{2}` is a root of \ :math:`x^2 - 2 = 0` and
therefore we say that \ :math:`\sqrt{2}` is an `algebraic number`.
With the ``nsimplify`` of ``sympy``, given an approximation for a root,
we can reconstruct the polynomial and thus find the
:index:`closest exact number` of the :index:`algebraic number`.
::
asqrt2 = numerical_approx(sqrt(2),10)
print(asqrt2)
from sympy import nsimplify
print(nsimplify(asqrt2,tolerance=0.001,full=True))
The last statement prints ``sqrt(2)``, obtained from a 10-bit
approximation ``1.4`` for \ :math:`\sqrt{2}`.
Assignments
-----------
1. Explain the outcome of ``3^4^5``. In particular, what is the order
of execution of the two exponentiation operations?
2. Write \ :math:`5^{4^3} - 1` as a product of prime numbers.
3. The greatest common divisor of two integer numbers
\ :math:`a` and \ :math:`b` can be written as a linear combination
(with integer coefficients \ :math:`k` and \ :math:`\ell`)
of \ :math:`a` and \ :math:`b`:
\ :math:`\gcd(a,b) = k a + \ell b`.
In Sage this is achieved with the command ``xgcd``.
Look in the help page of this command to compute
the coefficients of the linear combination of the
greatest common divisor of 12214 and 2012.
Give the command you type in to find these coefficients
and also give the command(s) to verify the result.
4. What is the difference in Sage between
\ :math:`1/3 + 1/3 + 1/3` and \ :math:`1.0/3 + 1.0/3 + 1.0/3`?
Explain.
5. Consecutive rational approximations for \ :math:`\pi`
are 3, 22/7, 333/106, 355/113, ...
In this sequence, what is the next more accurate rational
approximation for \ :math:`\pi`?
How many decimal places are correct in this next rational approximation?
6. Explain the difference between \ :math:`1.0 + 10^{-20}`
and \ :math:`1 + 10^{-20}`.
How can you make Sage return the same `correct` value of these
two sums?
7. The golden ratio is defined as \ :math:`r = \frac{1+\sqrt{5}}{2}`.
Give all Sage commands
(i) to compute a rational approximation for \ :math:`r`
accurate with three decimal places;
(ii) to show that the accuracy of this approximation
is indeed three decimal places;
(iii) to compute a sequence of the first ten consecutive rational
approximations, where the *k*-th number of the sequence is
accurate with *k* decimal places.
8. Consider ``r = 1.2599.``
Find the algebraic number that is closest to ``r``.
9. Consider :math:`\sqrt{3}`.
Compute the first ten terms in the continued fraction expansion.
Write the last element in the corresponding list of convergents.
Compare the floating-point approximation of this rational approximation
for :math:`\sqrt{3}` with the floating-point approximation
for :math:`\sqrt{3}`.
How many decimal places in the rational approximation are correct?
Give the floating-point approximation of the rational approximation
for :math:`\sqrt{3}`. Write only those decimals that are correct.