I always use PARI when I need to do computations and I am a big fan of this little program. I believe that it is possible to do in PARI everything you can do with such big programs as Maple and Mathematica. Well… almost everything. Here I’d like to present some tricks to do things in PARI that seem impossible from first sight, or just convenient hints. Readers are very welcome to publish their own tricks in comments. This way we may create something like a library of tricks.

#### Weird functions

Ok, which kind of functions you can use in PARI? Polynomials in several variables, rational functions in several variables are o.k., power series (in several variables also!). There is one thing you can do in one variable and cannot do with several. It is

#### Factoring polynomials.

Suppose you want to factor something like this:

Here is a solution, which I heard from Don Zagier. You substitute in place of y some big number. I think he prefers to use some big prime numbers, but to me any big number will do the job. Then factor the resulting polynomial as a polynomial in x, then if it does not factor you know it is irreducible. If it factors you try to guess the factors as polynomials in x and y:

(11:22) gp > x^2+2*x*y+y^2 %1 = x^2 + 2*y*x + y^2 (11:23) gp > factor(%) *** factor: sorry, factor for general polynomials is not yet implemented. (11:23) gp > subst(%,y,100000000) %2 = x^2 + 200000000*x + 10000000000000000 (11:23) gp > factor(%) %3 = [x + 100000000 2] (11:23) gp > %1/(x+y) %4 = x + y

#### Algebraic functions.

Sometimes you need to use functions like . In PARI you have Mods. Mod is an object in a finite extension of something. Their primary use, I guess, is for number fields, so that Mod(x, x^2+x+2) means . But nobody stops us from using it like Mod(y, y^2-x^2-1). So let us try:

(11:24) gp > y0=Mod(y,y^2-x^2-1) %5 = Mod(y, -x^2 + (y^2 - 1)) (12:12) gp > y0^2 %6 = Mod(y^2, -x^2 + (y^2 - 1))

Oops. We expected to get . This is the problem which everyone working with PARI should know about.

#### Side note: priority of variables

It is variable order. All variables are arranged in the order according to the time of first usage. Since x was used before y, it has ‘bigger priority’ than y. Therefore every expression in x and y is considered in first place as an expression in x, so the equation is treated like where is a kind of parameter. So this corresponds to . Therefore we should do it in a slightly different way:

(12:12) gp > y0=Mod(y,y^2-x0^2-1) %7 = Mod(y, y^2 + (-x0^2 - 1)) (12:19) gp > y0^2 %8 = Mod(x0^2 + 1, y^2 + (-x0^2 - 1))

Don’t forget to use ‘*lift*‘ when you want to get your final answer in a readable form.

#### Transcendental functions

In general: what can we do with transcendental functions? Since they are transcendental you cannot get any algebraic statements about them, so there is nothing to ask. There are, of course, exceptions to this. For example sometimes one transcendental function algebraically depends on another transcendental function. Like sine and cosine. Therefore the previous approach works perfectly well. You should encode sine and cosine in the following way:

(12:19) gp > Cos %9 = Cos (12:26) gp > Sin %10 = Sin (12:26) gp > Cos0=Mod(Cos, Cos^2+Sin^2-1) %11 = Mod(Cos, Cos^2 + (Sin^2 - 1))

Another exception is when you are expanding a transcendental function into a power series. This goes without problems:

(12:26) gp > sin(x) %12 = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9 - 1/39916800*x^11 + 1/6227020800*x^13 - 1/307674368000*x^15 + O(x^17)

Another thing you may want to do with a transcendental function is to differentiate it. Well, we come back to our example with Sin and Cos.

(12:28) gp > D(f)=subst(deriv(lift(f),Sin)*Cos-deriv(lift(f),Cos)*Sin, Cos, Cos0) (12:30) gp > D(Sin) %13 = Mod(Cos, Cos^2 + (Sin^2 - 1)) (12:31) gp > D(Cos0) %14 = -Sin

Now you can differentiate any combination of sine and cosine.

#### Differential equations

If you have a differential equation, like , you can encode it using variables y, dy and defining a differentiation operation like above which sends y to dy and dy to y:

(12:31) gp > D(f)=deriv(f,y)*dy+deriv(f,dy)*y

#### Numbers

There are several ways of dealing with algebraic numbers.

1. Using ‘Mod’. Just type Mod(x, x^2+x+2) when you need .

2. Using approximation. Simple approach: use (-1+sqrt(-7))/2 and in the end use very powerful algdep or lindep functions if you need the minimal equation:

(12:35) gp > (-1+sqrt(-7))/2 %15 = -1/2 + 1.322875655532295295250807877*I (12:41) gp > algdep(%,2) %16 = x^2 + x + 2

3. Using approximation, but with several embeddings of your number field in C.

#### Intersection theory

I have some experience computing some intersections of algebraic varieties. The approach is to use Mods to define algebraic varieties. Say, elliptic curve is Mod(y, x0^3+a*x0+b-y^2). For higher dimensional varieties one can use Mods with Mods inside. Then if you need to intersect something you simply get more equations. In the end you probably want to get points. Then the coordinates of these points will be some mods which give them as algebraic numbers. To find multiplicities solve all the equations in power series and look at the exponent of the main term. In the process of writing my thesis I was finding some points which are defined over some field of high degree (I think it was 12). If you try to use this approach don’t forget about the function ‘*polcompositum*‘, which helps, if you have numbers in different fields, to pass to the composite field.

#### Final remarks

The only thing that I not-so-like in PARI is its programming abilities. I think it is better to use some carefully designed standard scripting language for programming. That’s why I am working on integration of PARI with Python. This is still work in progress (look at some screenshots here).

If you know some other not-so-obvious tricks for PARI, you are very welcome to post them below.

## 8 comments

Comments feed for this article

November 21, 2007 at 12:06 am

The interesting bits of today — last week… « It’s Equal, but It’s Different…[...] Tricks using PARI [...]

November 25, 2007 at 2:54 am

PARI as good as Mathematica? « The Lyceum Mathematikoi[...] Posted in Computer Science at 2:54 am by saij I have yet to get into it, but Vivatsgas believes in it and lays down some tricks. [...]

July 5, 2008 at 8:51 pm

gcd calculator!Hi! I programmed an attractive online calculator that find the greatest common divisor(GCD) between two numbers. I will be happy, if you add the link in your blog. I hope that you and your visitors will enjoy!

—

http://gcd.awardspace.com

—

bye!

January 5, 2009 at 7:44 am

Jaime MontuertoHi,

I have this expression that’s beyond me, but relates to Fermat’s Last theorem ‘s diophantine equation, a^n + b^n = c^n, a,b,c,n all positive integers and n > 2. We all knew that it is now a theorem. I seems to me that it is related to this equation,

gcd(c^n – 1, a^n + b^n – 1) = 1

can you show or give a counterexample that will give a common factor larger than 1 that is positive ODD integer?

Thanks

September 29, 2009 at 7:02 pm

JohnHi Jaime,

Pick “a” and ‘c” odd, b even (for example a=c=3, b=2). Then trivially both sides will be even, and for all n you’ll have a GCD greater than 1.

gcd(3^3 – 1, 3^3 + 2^3 – 1) = gcd(26, 34) = 2.

September 29, 2009 at 7:03 pm

JohnOops, now I see the “odd” restriction, sorry!

November 23, 2009 at 4:07 pm

TempI require help with this :

If m1 & m2 are odd natural numbers and k1 & k2 are natural numbers, then does gcd(4^k1 + m1^2,4^k2+m2^2) always equal 1?

BTW, I arrived at that answer using the PolynomialGCD[] function in Mathematica 7 but I don’t really trust the program or my skills with it.(I’m not even sure if PolynomialGCD[] was the correct function to use for this kind of thing).Please verify this result for me.

Regards,

T.

December 14, 2011 at 8:55 am

Immobilier Djerba MidounJ’aime <3