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:

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 2n+O(log n) arithmetic operations.

The underlying method is described in:
M. A. Alekseyev and G. P. Michon. "Making Walks Count: From Silent Circles to Hamiltonian Cycles". In: J. Beineke, J. Rosenhouse (eds.) The Mathematics of Various Entertaining Subjects: Research in Recreational Math, Volume 2, Princeton University Press, 2017, pp. 157–168. ISBN 978-0-691-17192-0. doi:10.1515/9781400889136-012 doi:10.2307/j.ctt1s4773q.14 arXiv:1602.01396

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)); }

{ A137886(n) = nhp( matrix(2*n,2*n,i,j, min(i,j)<=n && max(i,j)>n && abs(j-i)!=n ) ); }

{ A137884(n) = nhp( matrix(3*n-1,3*n-1,i,j, abs(i-j)%3==1 ) ); }

II. Inversion of Multiplicative Functions

invphi.gp provides the following functions:

invphi(n) computes all solutions to the equation eulerphi(x)=n with respect to x
invphiNum(n) equivalent to but faster than #invphi(n)
invphiMin(n) equivalent to but faster than vecmin(invphi(n))
invphiMax(n) equivalent to but faster than vecmax(invphi(n))
invsigma(n,k) computes all solutions to the equation sigma(x,k)=n with respect to x. If k is omitted, it's assumed that k=1
invsigmaNum(n,k) equivalent to but faster than #invsigma(n,k)
invsigmaMin(n,k) equivalent to but faster than vecmin(invsigma(n,k))
invsigmaMax(n,k) equivalent to but faster than vecmax(invsigma(n,k))
invphitau(n,m) computes all solutions to the system { eulerphi(x)=n, numdiv(x)=m } with respect to x
                   
The implementation is based on the dynamic programming approach described in my paper:
M. A. Alekseyev, Computing the Inverses, their Power Sums, and Extrema for Euler's Totient and Other Multiplicative Functions. Journal of Integer Sequences 19 (2016), Article 16.5.2

Usage examples:

{ A014197(n) = invphiNum(n); }
{ A057635(n) = invphiMax(n); }

{ A072074(n) = invphiNum(10^n); }
{ A072075(n) = invphiMin(10^n); }
{ A072076(n) = invphiMax(10^n); }

{ A110078(n) = invsigmaNum(10^n); }
{ A110077(n) = invsigmaMin(10^n); }
{ A110076(n) = invsigmaMax(10^n); }

III. Binomial coefficients modulo integers

binomod.gp provides the following functions:

binomod(n,k,m) computes binomial(n,k) modulo integer m
binocoprime(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 input m is factored inside binocoprime(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.

The implementation is based on Lucas' Theorem and its generalization given in the paper:
Andrew Granville, The Arithmetic Properties of Binomial Coefficients. In: Proceedings of the Organic Mathematics Workshop, Simon Fraser University, December 12-14, 1995.

Usage example:

{ A080469() = for(n=2,10^9, if( !isprime(n) && binomod(3*n,n,n)==Mod(3,n)^n, print(n); ); ); }

{ A109760() = for(n=2,10^9, if( !isprime(n) && binomod(5*n,n,n)==Mod(5,n)^n, print(n); ); ); }

{ A109769() = for(n=2,10^9, if( !isprime(n) && binomod(7*n,n,n)==Mod(7,n)^n, print(n); ); ); }


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 = [ a1, a2, ..., ak ] computes the number of subgroups of the direct product of cyclic groups C(pa1) x C(pa2) x ... x C(pak). It implements the formula given in the paper:
G. A. Miller, On the Subgroups of an Abelian Group. The Annals of Mathematics, 2nd Ser. 6:1 (1904), 1-6.

Usage examples:

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

{ A061034(n) = my(f=factorint(n)); prod(i=1,#f~, vecmax(apply(x->numsubgrp(f[i,1],x),partitions(f[i,2]))) ); }

V. The Number of Self-Dual (Normal) Bases of GF(qm) 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(qm) 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(qm) Over GF(q). Proceedings of the American Mathematical Society 109:1 (1990), 23-29.

Usage examples:

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

VI. Sequences with Distinct Adjacent Elements

nseqadj.gp provides a function M(s) which, for a given k-dimensional vector s, computes the number of linear sequences, consisting elements from k classes with s[i] elements in the i-th class, where every pair of adjacent elements are from distinct classes. This script implements the formula given in the paper:
L. Q. Eifler, K. B. Reid Jr., D. P. Roselle,  Sequences with adjacent elements unequal. Aequationes Mathematicae 6:2-3 (1971), 256-262.

Usage example:

{ A110706(n) = M([n,n,n]); }

VII. Number of Monic Irreducible Multivariate Polynomials over Finite Fields

numirrpol.gp provides a function numirrpol(q,n,u) that counts the number of monic irreducible polynomials in n variables over GF(q) of degree at most u. Namely, it returns a vector of size u with the j-th component (j=1,2,...,u) equal to the number of such polynomials of degree j.
The implementation is based on the formula that I derived and posted in Russian forum in 2006. I also used it to contribute a whole bunch of related sequences (A115457 .. A115505) to OEIS. While I viewed this formula rather trivial and/or well-known, to my surprise the same formula was recently published in:
A. Bodin, Number of irreducible polynomials in several variables over finite fields. Amer. Math. Monthly 115 (2008), 653-660.
Interestingly, this publication even cites the sequence A115457 that I added to OEIS back in 2006.

Usage example:

{ A115457(n) = numirrpol(2,2,n); }

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

permod.gp and linrec.gp to come

IX. Empirical Recurrent Formulas with Polynomial Coefficients

isreccur.gp to come

X. Continued Fraction of Square Roots

powerful.gp to come


Back to Max's homepage.