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 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 at the official PARI/GP Resources page.

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.)

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

invphi(n) | computes all solutions to the
equation with respect to eulerphi(x)=nx |

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
with respect to sigma(x,k)=n. If x is
omitted, it's assumed that kk=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 with respect to {
eulerphi(x)=n, numdiv(x)=m }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.

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

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

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

Dieter Jungnickel, Alfred J. Menezes, Scott A. Vanstone, On the Number of Self-Dual Bases of GF(q

Usage examples:

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

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

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

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

Back to Max's homepage.