This site is supported by donations to The OEIS Foundation.
User:Max Alekseyev/gpscripts
This page was previously hosted at home.gwu.edu/~maxal/gpscripts/
Contents
 1 PARI/GP Scripts for Miscellaneous Math Problems
 1.1 Number of Hamiltonian paths and cycles in graphs
 1.2 Inversion of Multiplicative Functions
 1.3 Binomial coefficients modulo integers
 1.4 Number of subgroups of an abelian group
 1.5 The Number of SelfDual (Normal) Bases of GF(q^{m}) over GF(q)
 1.6 Sequences with Distinct Adjacent Elements
 1.7 Number of Monic Irreducible Multivariate Polynomials over Finite Fields
 1.8 Period of Linear Recurrent Sequences (e.g., Fibonacci Numbers) Modulo Primes
 1.9 Empirical Recurrent Formulas with Polynomial Coefficients
 1.10 Continued Fraction of Square Roots
PARI/GP Scripts for Miscellaneous Math Problems
PARI/GP is easy to use, fast, and powerful (freeware!) tool for numbertheoretical 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 inspired by certain computationally hard sequences in the Online 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 the 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 the 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 the 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.
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 9780691171920. doi:10.1515/9781400889136012 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*i1,2*j1)==1)); } { A107763(n) = nhc(matrix(n,n,i,j,gcd(2*i1,2*j1)==1)); } { A137886(n) = nhp( matrix(2*n,2*n,i,j, min(i,j)<=n && max(i,j)>n && abs(ji)!=n ) ); } { A137884(n) = nhp( matrix(3*n1,3*n1,i,j, abs(ij)%3==1 ) ); }
Inversion of Multiplicative Functions
invphi.gp provides the following functions:

invphi(n)
computes all solutions to the equationeulerphi(x)=n
with respect to x 
invphiNum(n)
equivalent to but faster than#invphi(n)

invphiMin(n)
equivalent to but faster thanvecmin(invphi(n))

invphiMax(n)
equivalent to but faster thanvecmax(invphi(n))

invsigma(n,k)
computes all solutions to the equationsigma(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 thanvecmin(invsigma(n,k))

invsigmaMax(n,k)
equivalent to but faster thanvecmax(invsigma(n,k))

invphitau(n,m)
computes all solutions to the system{ eulerphi(x)=n, numdiv(x)=m }
with respect to x 
invphitauNum(n,m)
equivalent to but faster than#invphitau(n,m)

invphitauMin(n,m)
equivalent to but faster thanvecmin(invphitau(n,m))

invphitauMax(n,m)
equivalent to but faster thanvecmax(invphitau(n,m))
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); }
Binomial coefficients modulo integers
binomod.gp provides the following functions:

binomod(n,k,m)
computesbinomial(n,k)
modulo integer m 
binocoprime(n,k,m)
tests ifbinomial(n,k)
is coprime to integer m 
binoval(n,k,p)
valuation ofbinomial(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 1214, 1995.
 Usage example
{ A080469() = forcomposite(n=2,10^9, if( binomod(3*n,n,n)==Mod(3,n)^n, print(n); ); ); } { A109760() = forcomposite(n=2,10^9, if( binomod(5*n,n,n)==Mod(5,n)^n, print(n); ); ); } { A109769() = forcomposite(n=2,10^9, if( binomod(7*n,n,n)==Mod(7,n)^n, print(n); ); ); }
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(p^{a1}) x C(p^{a2}) x ... x C(p^{ak}). 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), 16.
 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]))) ); }
The Number of SelfDual (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 selfdual and selfdual 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:
D. Jungnickel, A. J. Menezes, S. A. Vanstone, On the Number of SelfDual Bases of GF(q^{m}) Over GF(q). Proceedings of the American Mathematical Society 109:1 (1990), 2329.
 Usage examples
{ A088437(n) = sd(n,2); } { A135488(n) = sdn(n,2); }
Sequences with Distinct Adjacent Elements
nseqadj.gp provides a function M(s) which, for a given kdimensional vector s, computes the number of linear sequences, consisting elements from k classes with s[i] elements in the ith 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:23 (1971), 256262.
 Usage example
{ A110706(n) = M([n,n,n]); }
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 jth 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 dxdy.ru 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 wellknown, 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), 653660.
Interestingly, this publication even cites the sequence A115457 that I added to OEIS back in 2006.
 Usage example
{ A115457(n) = numirrpol(2,2,n); }
Period of Linear Recurrent Sequences (e.g., Fibonacci Numbers) Modulo Primes
permod.gp and linrec.gp to come
Empirical Recurrent Formulas with Polynomial Coefficients
isreccur.gp to come
Continued Fraction of Square Roots
powerful.gp to come