An Introduction to GP/PARI

Bryden Cais

October 5, 2002

1 What is PARI ?

In brief, PARI is a sophisticated calculator. What makes it different from other computation packages (like MAPLE, MATHEMATICA, MATLAB etc. is that PARI was designed by Number Theorists and professional mathematicians, not by entrepreneurs and engineers. Why should this matter? For starters, PARI has a lot of built in functions that are really useful to number theorists that some of the other programs either don’t have at all or have implemented in a cumbersome manner (I am thinking, for example, of PARI’s ability to work with p-adic numbers). More importantly,though, is that one gets the impression that PARI “thinks” like a mathematician. For example, when you work with fractions, PARI keeps them as quotients of integers, rather than expanding in decimal notation (unless you tell it to). But rather than give you more examples, let’s dive right in and see PARI at work.

2 Starting, Stopping and Entering Input

How you start PARI depends, of course, on what machine you are running it on. For a UNIX box, like we are using, the command (at the usual command prompt) is gp or gp-stat. Once PARI has started up, it will give you its own prompt, something like gp> or ?. To quit PARI type \q. All the usual arithmetic operations work just fine: for example, try

    567+129

and

    121*1729

Remenber to hit the enter key after you have finished typing your input and want PARI to start processing. Let’s try division. Input

    1728/156

Cool, PARI automatically reduces fractions to lowest common terms. You might try adding two fractions, like

    567/129+144/13

Again, we get the fraction in reduced form. What if you wanted it as a decimal expansion? The way to tell PARI that you are thinking in decimals is to put a decimal point after the numbers you enter. For example,

    567./129

or

    567/129.

both produce

    4.395348837209302325581395348

If you want more precision, thats easy; just input \p 50 for 50 places of accuracy, or \p 100 for 100 and so on. Try

    exp(Pi*sqrt(163))

Notice that since PARI is case sensitive, this is not the same thing as Exp(Pi*sqrt(163)) or exp(pi*sqrt(163)) etc., so be careful to get the cases right. What is going on? It looks like e^{π\sqrt{163}} is nearly an integer. This surely cannot be. To get more accuracy, type \p 100 and try again. Sure enough, its clear that e^{π\sqrt{163}} is not an integer. But why so close? Stick with number theory long enough and you’ll find out; the connections are deep and startling. OK, enough with the basics, lets move on.

3 Number Theory

Although PARI has a lot of mathematical functionality, we will be primarily interested in its number theory functions. Input ? and observe that PARI gives you a list of help topics. Now input ?4. PARI now lists all of the number theory related functions; to find out about a specific one, input a ? followed by the name. Try it. Let’s take a closer look at some of the most useful functions for our purposes.

3.1 factor

Clearly a useful command. Try

        factor(144169225289)

PARI gives the output

        [40163 1]
        [3589603 1].

This is a two by two matrix. In general, the output of factor(n) will be a r × 2 matrix, with the r factors of n being listed in the first column, and their corresponding multiplicities in the second column. The above matrix tells us that

144169225289 = 40163 · 3589603

and that both factors are prime. Try some more examples to get a feel for factoring integes. Just so you know, factor can factor other things, like polynomials over ℚ. Try it.

3.2 factormod, Mod, lift

This little function will factor a polynomial f(x) modulo a prime p. For example, if we want to factor x^{3} + 3x + 7 modulo 11, we input

        factormod(x^3+3*x+7,11)

and PARI returns

        [Mod(1, 11)*x + Mod(10, 11) 1]
        [Mod(1, 11)*x^2 + Mod(1, 11)*x + Mod(4, 11) 1].

Yikes! What does this mean? Those Mod’s are just saying that the first component is being considered modulo the second. For example, x=Mod(7,23) considers x as an integer modulo 23. Observe that Mod(7,23)=\Mod(30,23)=Mod(37,23) etc. To get rid of those annoying Mod’s, input lift(%) (here, the % symbol refers to the last output. Similarly, %k refers to k outputs ago). Then PARI returns the more sensible result

        [x + 10 1]
        [x^2 + x + 4 1],

telling us that

x^{3} + 3x + 7 ≡ (x + 10)(x^{2} + x + 4) mod 11.

Mod’s aren’t all bad, though, as they can be very useful for doing arithmetic mod n. For example, if you want to determine the inverse of 23 modulu 58, just input Mod(1/23,58) to find that the inverse is 53. Mod also works with polynomials. We’ll explore this later.

3.3 chinese

This is an implementation of the chinese remainder theorem. Let x=Mod(a,m) and y=Mod(b,n) be integers mod m, n respectively. Then chinese(x,y) gives an integer modulo mn whose reduction mod m is a and mod n is b (if the system can be solved). Let’s try an example. Input

            x=Mod(7,23)
            y=Mod(9,54)
            chinese(x,y).

PARI gives Mod(927,1242). We can check ourselves. Observe that fcator(1242) gives

        [23 1]
        [54 1],

as expected, and Mod(927,23) and Mod(927,54) give 7 and 9 respectively. What happens if the moduli are not relatively prime?

3.4 gcd

A simple function that gives the greatest common divisor of two integers. For example, input

            gcd(123456789,987654321)

to find that these two integers have GCD 9 (of course, you knew a priori that the two were divisible by 9, right?)

3.5 bezout

This function is an extension of gcd in that given integers m, n it solves the linear diophantine equation

mx + ny = gcd(m, n)

for integers x, y. Moreover, it gives the “smallest” (in an appropriate sense) solution. For example, input

            bezout(25678,166784)

PARI gives us

            [-37289, 5741, 2].

The first number is the coefficient of 25678, the second is the coefficient of 166784 and the last is gcd(25678, 166784). PARI is telling us that

25678(−37289) + 166784(5741) = 2.

Can you find all solutions to this linear diophantine equation? The function bezout also works with polynomials. Play with this if you like, and make sure you think about why the gcd of two polynomials is a well defined object.

3.6 lcm

I’ll bet you can guess what this does. Try it.

3.7 isprime

This function takes an integer input and returns 0 if the input is composite and 1 if it is prime. For example,

            isprime(257)
            isprime(256)

return 1, 0 respectively.

3.8 nextprime,prime

The command nextprime(x) gives the smallest prime bigger than or equal to x. For example,

            nextprime(Pi)
            nextprime(2346)

give 5, 2347 respectively. In a similar manner, prime(n) gives the n^{th} prime number.

3.9 eulerphi

This is just Euler’s “totient” function. Incidentially, one reference [1] explains

Euler’s generalization of the Fermat’s Little Theorem depends on a function which indeed was invented by Euler (1707-1783) but named by J.J.Sylvester (1814-1897) in 1883. I never saw an authoritative explanation for the name totient he has given the function. In Sylvestor’s opinion mathematics is essentially about seeing “differences in similarity, similarity in difference.” The word totient rhymes with quotient and the function has to do with division but in an unusual way. (Scott Brodie brougt to my attention that the Oxford English Dictionary brings up the latin root tot for adding up, total. Tot has also the meaning (of unkown origin) of a small child or a small drink. It would be very much in the spirit of the above maxim to use the word with two so different meanings. You could expect this from Sylvester whom E.T. Bell has christened hothead.)

Hm, OK. At any rate, we input

            eulerphi(345)

to find that there are 176 units modulo 345.

3.10 znprimroot

It is an important theorem that the group of units modulo n is cyclic if and only if n = p, p^{n} or 2p^{n} where p is prime. This means that there is one unit g called a primitive root whose successive powers give all the units modulo n. Let’s try it for n = 13. We input

            znprimroot(13)

and get

            Mod(2,13)

Now lets make a little table. Type ?vector to fin out how the vector command works and what it does. Once you understand this, input

            vector(12,n,Mod(2^n,13)).

PARI gives some expression with lots of Mod’s in it, so lift the result (using lift) to obtain the vector

            [2, 4, 8, 3, 6, 12, 11, 9, 5, 10, 7, 1],

which clearly contains all the units modulo 13.

3.11 znlog

Continuing with our example above, suppose that you know that g is a primitive root modulo n and you have some unit u modulo n and you want to know what power of g gives you u. This is the purpose of the function znlog. For example, input

            g=znprimroot(13)
            znlog(5,g)

to find that

2^{9} ≡ 5 mod 13

as you might as well check using Mod(2^9,5).

3.12 znorder

Still working in the unit group modulo n, this function gives you the order of u modulo n for any unit u. This is the smallest integer k such that

u^{k} ≡ 1 mod n.

(Why should such an integer always exist?) Be sure to input u in the form Mod(u,n) to tell PARI that we are working modulo n. For example,

            znorder(Mod(2,29))

returns 28, which incidentally tells us that 2 is a primitive root modulo 29.

4 Obtaining PARI

You can get PARI as an executable to install on your MAC or PC. The website is:

http://www.parigp-home.de/

Make sure to obtain the manuals, tutorial, and user’s guide as well as these are invaluable references.

5 Exercises

1. Factor 50260389661.

2. Factor x^{p} − x modulo p. Is this surprising?

3. Find all solutions to the linear diophantine equation

56752875x + 91796391y = gcd(56752875, 91796391).

4. Solve the system of equations

x ≡ i^{2} mod p_{i}

where p_{i} is the i^{th} prime and 1 ≤ i ≤ 20. Hint: Lookup the function if and write a recursive function. If you write any more than one line of code for this, you’re thinking along the wrong lines.

5. Make a table of values of φ(n), Euler’s function, for n < 50. Hint: use the vector function. Do you notice any patterns? If not, make a list of the factorization of n and of φ(n) for n < 50. Now do you notice anything?

6. Can you figure out what

        vector(20,n,sum(j=1,n,if(lift(Mod(n,j))==0,eulerphi(j),0)))

does? It might help to look up the syntax of if and sum. Any conjectures? Proofs?

7. What is the sum of the reciprocals of the first 40000 primes? Do you think this sum will converge or diverge? Hint: it might help to look at the syntax for the sum above, and to use the isprime function. If you want to be really slick and efficient, you shouldn’t do this, of course, and there is a better way using prime. Can you figure out how?

8. What proportion of the first 100 primes have 2 as a primitive root? How about the first 1000, 10000, 40000 primes? Can you make a bold conjecture?

References

[1] http://www.cut-the-knot.com/blue/Euler.shtml

[2]