Introduction to PARI/GP

  • PARI is a C library, allowing fast computations.
  • GP is the name of the scripting language giving access to the PARI functions.

There are many ways to use PARI/GP:

  • Jupyter interface to GP (the application that you are currently using)
  • In your browser using JavaScript
  • GP interactive shell (command-line application)
  • Python package cypari2
  • C code calling the PARI functions

Basic objects

In [1]:
57!
Out[1]:
40526919504877216755680601905432322134980384796226602145184481280000000000000
In [2]:
2 / 6
Out[2]:
1/3
In [3]:
(1+I)^2
Out[3]:
2*I
In [4]:
(x+1)^(-2)
Out[4]:
1/(x^2 + 2*x + 1)
In [5]:
Mod(2,5)^3
Out[5]:
Mod(3, 5)
In [6]:
Mod(x, x^2+x+1)^3
Out[6]:
Mod(1, x^2 + x + 1)
In [7]:
a = ffgen([3,5],'a) 
a^12  \\ in F_3^5
Out[7]:
2*a^4 + 2*a^3 + 2
In [8]:
Pi
Out[8]:
3.1415926535897932384626433832795028842
In [9]:
log(2)
Out[9]:
0.69314718055994530941723212145817656807
In [9]:
\p100
   realprecision = 115 significant digits (100 digits displayed)
In [10]:
log(2)
Out[10]:
0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875
In [11]:
exp(%)
Out[11]:
2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
In [12]:
log(1+x)
Out[12]:
x - 1/2*x^2 + 1/3*x^3 - 1/4*x^4 + 1/5*x^5 - 1/6*x^6 + 1/7*x^7 - 1/8*x^8 + 1/9*x^9 - 1/10*x^10 + 1/11*x^11 - 1/12*x^12 + 1/13*x^13 - 1/14*x^14 + 1/15*x^15 + O(x^16)
In [13]:
exp(%)
Out[13]:
1 + x + O(x^16)

Functions

In [13]:
?
Help topics: for a list of relevant subtopics, type ?n for n in
   0: user-defined functions (aliases, installed and user functions)
   1: PROGRAMMING under GP
   2: Standard monadic or dyadic OPERATORS
   3: CONVERSIONS and similar elementary functions
   4: functions related to COMBINATORICS
   5: NUMBER THEORETICAL functions
   6: POLYNOMIALS and power series
   7: Vectors, matrices, LINEAR ALGEBRA and sets
   8: TRANSCENDENTAL functions
   9: SUMS, products, integrals and similar functions
  10: General NUMBER FIELDS
  11: Associative and central simple ALGEBRAS
  12: ELLIPTIC CURVES
  13: L-FUNCTIONS
  14: MODULAR FORMS
  15: MODULAR SYMBOLS
  16: GRAPHIC functions
  17: The PARI community
Also:
  ? functionname (short on-line help)
  ?\             (keyboard shortcuts)
  ?.             (member functions)
Extended help (if available):
  ??             (opens the full user's manual in a dvi previewer)
  ??  tutorial / refcard / libpari (tutorial/reference card/libpari manual)
  ??  refcard-ell (or -lfun/-mf/-nf: specialized reference card)
  ??  keyword    (long help text about "keyword" from the user's manual)
  ??? keyword    (a propos: list of related functions).

Help

In [13]:
?4
binomial      fibonacci     hammingweight numbpart      numtoperm
partitions    permorder     permsign      permtonum     stirling

In [13]:
?atan
atan(x): arc tangent of x.

In [13]:
???determinant
algdisc           bnfsunit          charker           charpoly
ellheightmatrix   ellpadicregulator forsubgroup       matdet
matdetint         matdetmod         mathnfmod         matrixqz
mspolygon         nfdetint          nfhnfmod          polresultant
rnfdet            

See also:
  Finite abelian groups
  Pseudo-bases, determinant

Vectors and matrices

In [14]:
V = [1,2,3]
Out[14]:
[1, 2, 3]
In [15]:
W = [4,5,6]~
Out[15]:
[4, 5, 6]~
In [16]:
M = [1,2,3;4,5,6]
Out[16]:
[1 2 3]

[4 5 6]
In [17]:
V*W
Out[17]:
32
In [18]:
M*W
Out[18]:
[32, 77]~
In [19]:
U = [1..10]
Out[19]:
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

Components

In [20]:
V[2]
Out[20]:
2
In [21]:
W[1..2]
Out[21]:
[4, 5]~
In [22]:
M[2,2]
Out[22]:
5
In [23]:
M[1,]
Out[23]:
[1, 2, 3]
In [24]:
M[,2]
Out[24]:
[2, 5]~
In [25]:
M[1..2,1..2]
Out[25]:
[1 2]

[4 5]

Polymorphism

In [26]:
factor(91)
Out[26]:
[ 7 1]

[13 1]
In [27]:
factor(91+I)
Out[27]:
[      -1 1]

[   1 + I 1]

[ 4 + 5*I 1]

[1 + 10*I 1]
In [28]:
factor(x^4+4)
Out[28]:
[x^2 - 2*x + 2 1]

[x^2 + 2*x + 2 1]
In [29]:
factor((x^4+4)*I)
Out[29]:
[x + (-1 - I) 1]

[ x + (1 - I) 1]

[x + (-1 + I) 1]

[ x + (1 + I) 1]
In [30]:
factor((x^4+1)*Mod(1,y^2-2))
Out[30]:
[x^2 + Mod(-y, y^2 - 2)*x + 1 1]

[ x^2 + Mod(y, y^2 - 2)*x + 1 1]
In [31]:
factor((x^4+4)*Mod(1,13))
Out[31]:
[Mod(1, 13)*x + Mod(4, 13) 1]

[Mod(1, 13)*x + Mod(6, 13) 1]

[Mod(1, 13)*x + Mod(7, 13) 1]

[Mod(1, 13)*x + Mod(9, 13) 1]

Numerical integration

In [31]:
\p38
   realprecision = 38 significant digits
In [32]:
intnum(x=0,1,1/(1+x^2))/Pi
Out[32]:
0.25000000000000000000000000000000000000
In [33]:
sumnum(n=1,1/n^2)/Pi^2
Out[33]:
0.16666666666666666666666666666666666667
In [34]:
sumalt(n=1,(-1)^n*log(n))
Out[34]:
0.22579135264472743236309761494744107198
In [35]:
exp(2*%)
Out[35]:
1.5707963267948966192313216916397514427

Comprehension

In [36]:
[n^2|n<-[1..10]]
Out[36]:
[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
In [37]:
[n^2|n<-[1..10],isprime(n)]
Out[37]:
[4, 9, 25, 49]
In [38]:
[a,b] = [1,2];
In [38]:
print("a=",a," b=",b)
a=1 b=2

Programming in GP

In [39]:
dist(a,b) = sqrt(a^2+b^2);
In [40]:
dist(1,2)
Out[40]:
2.2360679774997896964091736687312762354

In a GP file, line ending terminates the line unless the they are preceded by = or \ or are in a section delimited by braces:

In [41]:
cmpz(a,b) = if (abs(a)>abs(b), print(a), print(b))
Out[41]:
(a,b)->if(abs(a)>abs(b),print(a),print(b))
In [42]:
cmpz(a,b) =
  if (abs(a)>abs(b), print(a), print(b))
Out[42]:
(a,b)->if(abs(a)>abs(b),print(a),print(b))
In [43]:
cmpz(a,b) = if (abs(a)>abs(b), \
              print(a), \
              print(b))
Out[43]:
(a,b)->if(abs(a)>abs(b),print(a),print(b))
In [44]:
cmpz(a,b) =
{
  if (abs(a)>abs(b),
    print(a)
   ,print(b))
}
Out[44]:
(a,b)->if(abs(a)>abs(b),print(a),print(b))

Example of a function

  • Put the opening brace on the line after the = sign.
  • End the function by a semicolon.
  • Declare any local variables with my().
  • Do not declare the loop index: it is local to the loop anyway.
  • The function return value is the last computed value.
  • Indent you code following the example.
In [45]:
fibo(n)=
{
  my(u0=0,u1=1);
  for(i=2,n,
    [u0,u1]=[u1,u0+u1]);
  u1;
}
In [46]:
fibo(100)
Out[46]:
354224848179261915075

While loop

In [47]:
rho(n)=
{
  my(x=2,y=5);
  while(gcd(y-x,n)==1,
      x=(x^2+1)%n;
      y=(y^2+1)%n; y=(y^2+1)%n
      );
  gcd(n,y-x);
}
In [48]:
rho(2^64+1)
Out[48]:
274177

Control flow

In [49]:
wieferich(n)=
{
  forprime(p=2,n,
    if(Mod(2,p^2)^(p-1)==1,
      return(p)));
}
In [50]:
wieferich2(n)=
{
  my(r);
  forprime(p=2,n,
    if(Mod(2,p^2)^(p-1)==1,r=p;break));
  r;
}
In [51]:
wieferich(10000)
Out[51]:
1093
In [52]:
wieferich2(10000)
Out[52]:
1093

Constructors

In [53]:
V=vector(10,i,1/i)
Out[53]:
[1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10]
In [54]:
[1/i|i<-[1..10]]
Out[54]:
[1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10]
In [55]:
M=matrix(4,4,i,j,i*j)
Out[55]:
[1 2  3  4]

[2 4  6  8]

[3 6  9 12]

[4 8 12 16]

forvec

Instead of

In [56]:
s3(n)=
{
  my(m=sqrtint(n));
  for (i=1,m,
    for (j=1,m,
      for (k=1,m,
        if (i^2+j^2+k^2==n,
          return([i,j,k])))));
}
s3(12345)
Out[56]:
[4, 77, 80]

use forvec:

In [57]:
s3(n)=
{
  my(m=sqrtint(n));
  forvec(v=vector(3,i,[1,m]),
    if (v*v~==n,
      return(v)));
}
s3(12345)
Out[57]:
[4, 77, 80]

For a better algorithm, see qfsolve.

Associative arrays

In [58]:
birthday(n)=
{
  my(M = Map());
  for(i=1,oo,
    my(x=random(n), j);
    if(mapisdefined(M,x,&j),
      return([i,j]));
    mapput(M,x,i));
}
In [59]:
birthday(2^30)
Out[59]:
[23001, 18881]

Algebraic Number Theory

Irreducibility

In GP, we describe a number field $K$ as $K = Q[x]/f(x)$ where $f \in \mathbb{Z}[x]$ is a monic irreducible polynomial.

In [60]:
f = x^4 - 2*x^3 + x^2 - 5
Out[60]:
x^4 - 2*x^3 + x^2 - 5
In [61]:
polisirreducible(f)
Out[61]:
1

GP knows cyclotomic polynomials:

In [62]:
g = polcyclo(30)
Out[62]:
x^8 + x^7 - x^5 - x^4 - x^3 + x + 1

Polmod

To perform simple operations in $K = \mathbb{Q}[x]/f(x) = \mathbb{Q}(\alpha)$ where $f(\alpha) = 0$, we can use Mod:

In [63]:
Mod(x,f)^5
Out[63]:
Mod(3*x^3 - 2*x^2 + 5*x + 10, x^4 - 2*x^3 + x^2 - 5)

Interpretation: $\alpha^5 = 3 \alpha^3 - 2 \alpha^2 + 5 \alpha + 10$. We check that the roots of $g$ are 30th roots of unity:

In [64]:
lift(Mod(x,g)^15)
Out[64]:
-1

We used lift to make the output more readable.

polredbest

Sometimes we can find a simpler defining polynomial for the same number field by using polredbest:

In [65]:
h = x^5 + 7*x^4 + 22550*x^3 - 281686*x^2 - 85911*x + 3821551;
In [66]:
polredbest(h)
Out[66]:
x^5 - x^3 - 2*x^2 + 1

Interpretation: $\mathbb{Q}[x]/h(x) \cong \mathbb{Q}[x]/(x^5 - x^3 - 2 x^2 + 1)$.

nfinit and precomputed information

Most operations on number fields require a precomputation, which is performed by the initialisation function nfinit.

In [67]:
K = nfinit(f);

K contains the precomputation of the number field K = \mathbb{Q}[x]/f(x).

In [68]:
 K.pol
Out[68]:
x^4 - 2*x^3 + x^2 - 5
In [69]:
K.sign
Out[69]:
[2, 1]

K has signature (2, 1): it has two real embeddings and one pair of conjugate complex embeddings.

In [70]:
K.disc
Out[70]:
-1975
In [71]:
K.zk
Out[71]:
[1, 1/2*x^2 - 1/2*x - 1/2, x, 1/2*x^3 - 1/2*x^2 - 1/2*x]
In [72]:
w = K.zk[2]
Out[72]:
1/2*x^2 - 1/2*x - 1/2

K has discriminant $-1975$, and its ring of integers is $$ \mathbb{Z}_K = \mathbb{Z} + \mathbb{Z}\frac{\alpha^2 - \alpha - 1}{2} + \mathbb{Z}\alpha + \mathbb{Z}\frac{\alpha^3 -\alpha^2 - \alpha}{2} $$

In [ ]: