# 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 ASCII plot is shown in Fig. 96.

To approximate \(\pi\) we could use

but this series converges very slowly. We check a slow sum to evaluate \(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.

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 \(f(z)\) in a disk \(C_{a,r}\) in the complex plane centered at \(z=a\) and with radius \(r\) is

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 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 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¶

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

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.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.Consider

`r = 1.2599210498948732.`

Find the algebraic number that is closest to`r.`