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.

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

To approximate \(\pi\) we could use

\[\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 \(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

\[\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 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

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

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.