PARI/GP Scripts by Max Alekseyev

PARI/GP Scripts for Miscellaneous Math Problems

by Max Alekseyev

PARI/GP is easy to use, fast, and powerful (freeware!) tool for number-theoretical computations. As a rather frequent user of PARI/GP, I have developed a number of advanced scripts that I would like to share with the world. The most useful of them (at my humble opinion) are presented at this page.

Development of these scripts in many cases was insprired by certain computationally hard sequences in the On-line Encyclopedia of Integer Sequences (OEIS) and occasionally resulted in extension of such sequences. Whenever appropriate I illustrate the scripts usage with simple programs computing sequences in OEIS.

I welcome comments and suggestions regarding the scripts presented at this page. In particular, please let me know about:

bugs and miscalculations (a must to report!).
efficiency / performance concerns, both algorithmic and practical. Please tell me if you know a way to speed up things or a faster algorithm than currently implemented.
PARI/GP implementation issues including programming style, tips'n'tricks, etc.
your experience with these scripts. I would be pleased to know if you found them applicable to a computational problem of your interest (e.g., if they helped to extend a sequence in OEIS).

N.B. Some of the scripts planned for publication are not yet in a publishable form, in which case only an announce is given. Please email me, if you have an urgent need to try out a script that is not yet present.

P.S. Please also notice many other useful PARI/GP Contributed Scripts at the official PARI/GP Resources page.

I. Number of Hamiltonian paths and cycles in graphs

hamiltonian.gp provides two functions nhp(A) and nhc(A) that compute the number of Hamiltonian paths and cycles respectively in the graph defined by the adjacency n × n matrix A. The running time complexity is 2^{n+O(log n)} arithmetic operations.

Usage examples:

{ A076220(n) = nhp(matrix(n,n,i,j,gcd(i,j)==1)) }

{ A086595(n) = nhc(matrix(n,n,i,j,gcd(i,j)==1)) }

{ A103839(n) = nhp(matrix(n,n,i,j,isprime(i+j))) }

{ A107761(n) = nhp(matrix(n,n,i,j,gcd(2*i-1,2*j-1)==1)) }

{ A107763(n) = nhc(matrix(n,n,i,j,gcd(2*i-1,2*j-1)==1)) }

II. Inversion of Euler Totient Function

invphi.gp provides two functions invphi(n) and numinvphi(n) computing all solutions x to eulerphi(x)=n and the number of such solutions respectively. While the functionality of numinvphi(n) is identical to length(invphi(n)), the former is a bit faster.

Usage example:

{ A014197(n) = numinvphi(n) }

{ A057635(n) = trap(,0,vecmax(invphi(n))) }

III. Binomial Coefficients Modulo Primes

binomod.gp provides the following functions:

binomod(n,k,p)	computes binomial(n,k) modulo prime p (using Lucas theorem)
binomod1(n,k,p)	tests if binomial(n,k) is co-prime to prime p
binomod2(n,k,m)	tests if binomial(n,k) is co-prime to integer m
binoval(n,k,p)	valuation of binomial(n,k) with respect to prime p

Notes. The functions binomod1(n,k,p) and binomod2(n,k,m) are variants of binomod(n,k,p) that are a bit faster for their specific tasks. The input p is not tested for primality; if it is not prime, the result may be incorrect. The input m is factored inside binomod2(n,k,m) and that may take time for large m. The functionality of binoval(n,k,p) is equivalent to valuation(binomial(n,k),p) but the former does not compute binomial(n,k) explicitly and takes only O(log(n)) arithmetic operations.

For background on computation of binomial coefficient modulo prime numbers (and possible directions for further development) see the woderful paper: "The Arithmetic Properties of Binomial Coefficients" by Andrew Granville.

IV. Number of subgroups of an abelian group

ngs.gp provides the function numsubgrp(p,a) that for a prime p and a vector a = [ a₁, a₂, ..., a_k ] computes the number of subgroups of the direct product of cyclic groups C(p^a₁) x C(p^a₂) x ... x C(p^a_k). It implements the formula given in the paper: G. A. Miller "On the Subgroups of an Abelian Group", The Annals of Mathematics, 2nd Ser. (1904), 6 (1), 1-6.

Usage examples:

{ A006116(n) = numsubgrp(2,vector(n,i,1)) }

{ A061034(n) = local(f=factorint(n)); prod(i=1, matsize(f)[1], A061034pp(f[i,1],f[i,2]) ) }
{ A061034pp(p,k) = res=0; for(i=1, k, aux_part(p, k-i, i, [])); res } \\ for prime power p^k
{ aux_part(p, n, m, v) = \\ iterate over all partitions
v = concat(v,m);
if(n,
for(i=1, min(m,n), aux_part(p, n-i, i, v))
,
res=max(res,numsubgrp(p,v));
); }

V. The Number of Self-Dual (Normal) Bases of GF(q^m) over GF(q)

nsdb.gp provides two functions sd(m,q) and sdn(m,q) that compute respectively the number of distinct self-dual and self-dual normal bases of the finite field GF(q^m) over GF(q). The number q is a power of prime in sd(m,q) and a prime in sdn(m,q). This script implements the formulae given in the paper: Dieter Jungnickel, Alfred J. Menezes, Scott A. Vanstone "On the Number of Self-Dual Bases of GF(q^m) Over GF(q)". Proceedings of the American Mathematical Society (1990), 109 (1), 23-29.

Usage examples:

{ A088437(n) = sd(n,2) }
{ A135488(n) = sdn(n,2) }

VI. Period of Linear Recurrent Sequences (e.g., Fibonacci Numbers) Modulo Primes

permod.gp and linrec.gp to come

VII. Empirical Recurrent Formulas with Polynomial Coefficients

isreccur.gp to come

VIII. Generation of Integer Partitions

partition.gp to come

Meanwhile, please see Richard Mathar's implementation at:
http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=671

IX. Continued Fraction of Square Roots

powerful.gp to come

X. Number of Monic Irreducible Multivariate Polynomials over Finite Fields

irreduc.gp to come

Back to Max's homepage.