Lecture 41: Higher Arithmetic with Pari/GP
==========================================

PARI stems from Pascal ARIthmetic (although switched to C,
now a pun on pari de Pascal),
and GP means the Great Programmable calculator,
distributed under the GNU General Public License.

PARI/GP is a widely used computer algebra system designed
for fast computations in number theory (factorizations, algebraic
number theory, elliptic curves...), but also contains a large number
of other useful functions to compute with mathematical entities such
as matrices, polynomials, power series, algebraic numbers, etc.,
and a lot of transcendental functions.

Calculating with GP in Sage
---------------------------

Let us first see how we can call GP in a Sage session.
We first set the precision to 100 decimal places.

::

   gp('default(realprecision, 100)')
   pi100 = gp('Pi')
   print(pi100, 'has type', type(pi100))
   spi100 = pi100.sage()
   print(spi100, 'has type', type(spi100))

The type of ``pi100`` is a ``GpElement``
and after application of the ``sage()`` method
we obtain a real MPFR number.

We can make ASCII plots.
In a Sage worksheet, in Settings the
``Number of word-wrap columns`` had to be changed to 100
for the plot to show up well.

::

   ap = gp.eval("plot(x=0,2*Pi,sin(x))")
   ap

The :index:`ASCII plot` is shown in :numref:`figasciisinplot`.

   .. _figasciisinplot:

   .. figure:: ./figasciisinplot.png
       :width: 75%
       :align: center

       An ASCII plot of sin(x) for x from 0 to 2*Pi.

To approximate :math:`\pi` we could use

.. math::
   \frac{\pi}{4}  =  \sum_{k=0}^\infty \frac{(-1)^k}{2 k + 1}
   =  1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots

but this series converges very slowly.
We check a slow sum to evaluate :math:`pi` below.

::

   gp('default(realprecision, 10)')
   gp('sum4pi = 4*sum(k=0,10000,(-1.0)^k/(2*k+1))')
   sum4pi = gp('sum4pi')
   print(sum4pi)
   print(pi.n(digits=10))

and we see ``3.141692644`` and ``3.141592654.``
We can accelerate a slowly alternating sum.
First we reset the default precision to 100 decimal places.

::

   gp('default(realprecision, 100)')
   print(gp('sum4pi2 = 4*sumalt(k=0,(-1.0)^k/(2*k+1))'))
   print(gp('Pi'))

and we see the same decimal expansion twice.

.. index:: rational approximation

Rational approximations can be obtained with ``bestappr().``

::

   gp('bestappr(Pi)')

As second argument of ``bestappr()`` we can give the size of the denominator.

::

   [gp('bestappr(sqrt(2),10^%d)' % k) for k in range(1,11)]

We obtain the list

::

   [7/5,
   99/70,
   1393/985,
   8119/5741,
   114243/80782,
   665857/470832,
   9369319/6625109,
   131836323/93222358,
   768398401/543339720,
   6333631924/4478554083].

The Cauchy Integral Formula
---------------------------

The number of roots of \ :math:`f(z)` in a 
disk \ :math:`C_{a,r}` in the complex plane
centered at \ :math:`z=a` and with radius \ :math:`r` is

.. math::

   \frac{1}{2 \pi i} \oint_{C_{a,r}} \frac{f'(z)}{f(z)} dx

We first define a polynomial.  Notice the relationship between functions
and polynomials, which is very similar as in Sage.

::

   gp('f(z) = z^3 - z')
   print gp('f')
   gp('df = deriv(f(z),z)')
   print(gp('df'))
   gp('g(x) = substpol(df,z,x)')
   print(gp('g'))

The ``substpol`` is a way to turn an expression into a function,
without retyping the expression.
To be sure that we have functions, we evaluate at a point.

::

   print(gp('f(2)'))
   print(gp('g(2)'))

Now we apply the Cauchy integral formula.

::

   gp('intcirc(z=0, 0.5, g(z)/f(z))')

There is one root in a disk centered at z = 0 of radius 0.5.
Because f is a cubic polynomial, if we enlarge the radius enough,
we know the radius of the disk that contains all roots.

::

   gp('intcirc(z=0, 1.5, g(z)/f(z))')

Expansions and Series
---------------------

The syntax for expansions and :index:`series` is interesting.

::

   print(gp('142 + O(2^10)'))
   print(gp('142 + O(3^10)'))
   print(gp('142 + O(10^10)'))

and we see three different expansions of 142
as 2 + 2^2 + 2^3 + 2^7 + O(2^10), and
1 + 2*3 + 2*3^3 + 3^4 + O(3^10), and the familiar
2 + 4*10 + 10^2 + O(10^10).

With ``truncate()`` we can evaluate the series.

::

   print(gp('c = 142 + O(5^10)'))
   print(gp('truncate(c)'))

The 2 + 3*5 + 5^3 + O(5^10) is turned into 142.

The same syntax applies to :index:`Taylor series`,
about x = 0, of tenth order.

::

   gp('t = sin(x) + O(x^10)')

The ``truncate()`` on t will return a polynomial.

::

   gp('truncate(t)')

Integer Relation Detection
--------------------------

We can detect integer relations between floating-point approximations.
As a test we take the relation ``log(15) = log(3*5) = log(3) + log(5).``

::

   gp('u = [log(15.0), log(3.0), log(5.0)]')
   gp('u')

We have the list ``[2.708050201, 1.098612289, 1.609437912]``
of three floating-point numbers.

::

   gp('r = lindep(u)')

The result ``[-1, 1, 1]~``
indicates that ``-log(15) + log(3) + log(5) = 0.``
The flag ``~`` at the end is there because the output is a column vector.

::

   gp('type(r)')

And we see ``t_COL`` as the type.
We can look for algebraic dependencies between numbers,
which is equivalent to calling lindep on the list ``[1, x, x^2, .. , x^k].``

::

   gp('x = 3^(1/6)')

Let us check if GP can detect from the value
``1.200936955`` that ``x`` is an algebraic number,
that is: that ``x`` is a root of a polynomial.

::

   gp('algdep(x,6)')

On return we have ``x^6 - 3``.

Assignments
-----------

1. Compute the first 10 best rational approximations for 
   the transcendental number :math:`e = \exp(1)`.
   The :math:`k`-th number in the sequence has :math:`k` decimal
   places in the denominator.

2. Is computing the Taylor series of ``sin(x)*cos(x)``
   the same as multiplying the Taylor series of ``sin(x)``
   with the Taylor series of ``cos(x)``?
   In your answer, calculate with Taylor series of order 10.

3. Use the Cauchy integral formula to show that ``sin(3*Pi*x)``
   has 2 roots in a disk of the complex plane centered at the origin
   and with radius 0.5.

4.  Consider ``r = 1.2599210498948732.``  
    Find the algebraic number that is closest to ``r.``