Pari/GP Reference Documentation | Contents - Global index - GP keyboard shortcuts |
break breakpoint dbg_down dbg_err dbg_up dbg_x for forcomposite fordiv fordivfactored foreach forell forfactored forpart forperm forprime forprimestep forsquarefree forstep forsubgroup forsubset forvec if iferr next return setdebug until while | |
A number of control statements are available in GP. They are simpler and have a syntax slightly different from their C counterparts, but are quite powerful enough to write any kind of program. Some of them are specific to GP, since they are made for number theorists. As usual, X will denote any simple variable name, and seq will always denote a sequence of expressions, including the empty sequence. Caveat. In constructs like
for (X = a,b, seq)
the variable
n = 5;
for (n = 1, 10,
if (something_nice(), break);
);
\\ at this point
If the sequence
? for (n = 1, 10, n += 2; print(n)) 3 6 9 12
| |
break({n = 1}) | |
Interrupts execution of current seq, and immediately exits from the n innermost enclosing loops, within the current function call (or the top level loop); the integer n must be positive. If n is greater than the number of enclosing loops, all enclosing loops are exited.
| |
breakpoint() | |
Interrupt the program and enter the breakloop. The program continues when the breakloop is exited.
? f(N,x)=my(z=x^2+1);breakpoint();gcd(N,z^2+1-z); ? f(221,3) *** at top-level: f(221,3) *** ^ — — -- *** in function f: my(z=x^2+1);breakpoint();gcd(N,z *** ^ — — — — — — -- *** Break loop: type <Return> to continue; 'break' to go back to GP break> z 10 break> %2 = 13
| |
dbg_down({n = 1}) | |
(In the break loop) go down n frames. This allows to cancel a previous
call to
? x = 0; ? g(x) = x-3; ? f(x) = 1 / g(x+1); ? for (x = 1, 5, f(x+1)) *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - *** in function f: 1/g(x+1) *** ^ — — - *** _/_: impossible inverse in gdiv: 0. *** Break loop: type 'break' to go back to GP prompt break> dbg_up(3) \\ go up 3 frames *** at top-level: for(x=1,5,f(x+1)) *** ^ — — — — — -- break> x 0 break> dbg_down() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - break> x 1 break> dbg_down() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - break> x 1 break> dbg_down() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - *** in function f: 1/g(x+1) *** ^ — — - break> x 2
The above example shows that the notion of GP frame is
finer than the usual stack of function calls (as given for instance by the
GDB
| |
dbg_err() | |
In the break loop, return the error data of the current error, if any.
See
? iferr(1/(Mod(2,12019)^(6!)-1),E,Vec(E)) %1 = ["e_INV", "Fp_inv", Mod(119, 12019)] ? 1/(Mod(2,12019)^(6!)-1) *** at top-level: 1/(Mod(2,12019)^(6!)- *** ^ — — — — — — -- *** _/_: impossible inverse in Fp_inv: Mod(119, 12019). *** Break loop: type 'break' to go back to GP prompt break> Vec(dbg_err()) ["e_INV", "Fp_inv", Mod(119, 12019)]
| |
dbg_up({n = 1}) | |
(In the break loop) go up n frames, which allows to inspect data of the
parent function. To cancel a
? x = 0; ? g(x) = x-3; ? f(x) = 1 / g(x+1); ? for (x = 1, 5, f(x+1)) *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - *** in function f: 1/g(x+1) *** ^ — — - *** _/_: impossible inverse in gdiv: 0. *** Break loop: type 'break' to go back to GP prompt break> x 2 break> dbg_up() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - break> x 1 break> dbg_up() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - break> x 1 break> dbg_up() *** at top-level: for(x=1,5,f(x+1)) *** ^ — — — — — -- break> x 0 break> dbg_down() \\ back up once *** at top-level: for(x=1,5,f(x+1)) *** ^ — — - break> x 1
The above example shows that the notion of GP frame is
finer than the usual stack of function calls (as given for instance by the
GDB
| |
dbgx(A, {n}) | |
Print the inner structure of A, complete if n is omitted, up
to level n otherwise. This function is useful for debugging. It is similar
to
| |
for(X = a, b, seq) | |
Evaluates seq, where
the formal variable X goes from a to b, where a and b must be in
ℝ. Nothing is done if a > b. If b is set to
| |
forcomposite(n = a, {b}, seq) | |
Evaluates seq, where the formal variable n ranges over the composite numbers between the nonnegative real numbers a to b, including a and b if they are composite. Nothing is done if a > b.
? forcomposite(n = 0, 10, print(n)) 4 6 8 9 10
Omitting b means we will run through all composites ≥ a,
starting an infinite loop; it is expected that the user will break out of
the loop himself at some point, using Note that the value of n cannot be modified within seq:
? forcomposite(n = 2, 10, n = []) *** at top-level: forcomposite(n=2,10,n=[]) *** ^ — *** index read-only: was changed to [].
| |
fordiv(n, X, seq) | |
Evaluates seq, where
the formal variable X ranges through the divisors of n
(see
This routine uses To avoid storing all divisors, possibly using a lot of memory, the following (slower) routine loops over the divisors using essentially constant space:
FORDIV(N)= { my(F = factor(N), P = F[,1], E = F[,2]); forvec(v = vector(#E, i, [0,E[i]]), X = factorback(P, v)); } ? for(i=1, 10^6, FORDIV(i)) time = 11,180 ms. ? for(i=1, 10^6, fordiv(i, d, )) time = 2,667 ms. Of course, the divisors are no longer sorted by inreasing size.
| |
fordivfactored(n, X, seq) | |
Evaluates seq, where
the formal variable X ranges through [d,
It is assumed that
This function is particularly useful when n is hard to factor and one must evaluate multiplicative function on its divisors: we avoid refactoring each divisor in turn. It also provides a small speedup when n is easy to factor; compare
? A = 10^8; B = A + 10^5; ? for (n = A, B, fordiv(n, d, eulerphi(d))); time = 2,091 ms. ? for (n = A, B, fordivfactored(n, d, eulerphi(d))); time = 1,298 ms. \\ avoid refactoring the divisors ? forfactored (n = A, B, fordivfactored(n, d, eulerphi(d))); time = 1,270 ms. \\ also avoid factoring the consecutive n's !
| |
foreach(V, X, seq) | |
Evaluates seq, where the formal variable X ranges through the
components of V (
| |
forell(E, a, b, seq, {flag = 0}) | |
Evaluates seq, where the formal variable E = [name, M, G] ranges through all elliptic curves of conductors from a to b. In this notation name is the curve name in Cremona's elliptic curve database, M is the minimal model, G is a ℤ-basis of the free part of the Mordell-Weil group E(ℚ). If flag is nonzero, select only the first curve in each isogeny class.
? forell(E, 1, 500, my([name,M,G] = E); \ if (#G > 1, print(name))) 389a1 433a1 446d1 ? c = 0; forell(E, 1, 500, c++); c \\ number of curves %2 = 2214 ? c = 0; forell(E, 1, 500, c++, 1); c \\ number of isogeny classes %3 = 971
The
The library syntax is
| |
forfactored(N = a, b, seq) | |
Evaluates seq, where
the formal variable N is [n,
This function is only implemented for |a|, |b| < 264 (232 on a 32-bit
machine). It uses a sieve and runs in time O(sqrt{b} + b-a). It should
be at least 3 times faster than regular factorization as long as the interval
length b-a is much larger than sqrt{b} and get relatively faster as
the bounds increase. The function slows down dramatically
if
? B = 10^9; ? for (N = B, B+10^6, factor(N)) time = 4,538 ms. ? forfactored (N = B, B+10^6, [n,fan] = N) time = 1,031 ms. ? B = 10^11; ? for (N = B, B+10^6, factor(N)) time = 15,575 ms. ? forfactored (N = B, B+10^6, [n,fan] = N) time = 2,375 ms. ? B = 10^14; ? for (N = B, B+10^6, factor(N)) time = 1min, 4,948 ms. ? forfactored (N = B, B+10^6, [n,fan] = N) time = 58,601 ms.
The last timing is with the default
Note that all PARI multiplicative functions accept the
? s = 0; forfactored(N = 1, 10^7, s += moebius(N)*eulerphi(N)); s time = 6,001 ms. %1 = 6393738650 ? s = 0; for(N = 1, 10^7, s += moebius(N)*eulerphi(N)); s time = 28,398 ms. \\ slower, we must factor N. Twice. %2 = 6393738650 The following loops over the fundamental dicriminants less than X:
? X = 10^8; ? forfactored(d=1,X, if (isfundamental(d),)); time = 34,030 ms. ? for(d=1,X, if (isfundamental(d),)) time = 1min, 24,225 ms.
| |
forpart(X = k, seq, {a = k}, {n = k}) | |
Evaluate seq over the partitions X = [x1,...xn] of the
integer k, i.e. increasing sequences x1 ≤ x2... ≤ xn of sum
x1+...+ xn = k. By convention, 0 admits only the empty partition and
negative numbers have no partitions. A partition is given by a
? forpart(X=4, print(X)) Vecsmall([4]) Vecsmall([1, 3]) Vecsmall([2, 2]) Vecsmall([1, 1, 2]) Vecsmall([1, 1, 1, 1]) Optional parameters n and a are as follows: * n = nmax (resp. n = [nmin,nmax]) restricts partitions to length less than nmax (resp. length between nmin and nmax), where the length is the number of nonzero entries. * a = amax (resp. a = [amin,amax]) restricts the parts to integers less than amax (resp. between amin and amax). By default, parts are positive and we remove zero entries unless amin ≤ 0, in which case we fix the size #X = nmax:
\\ at most 3 nonzero parts, all <= 4 ? forpart(v=5,print(Vec(v)), 4, 3) [1, 4] [2, 3] [1, 1, 3] [1, 2, 2] \\ between 2 and 4 parts less than 5, fill with zeros ? forpart(v=5,print(Vec(v)),[0,5],[2,4]) [0, 0, 1, 4] [0, 0, 2, 3] [0, 1, 1, 3] [0, 1, 2, 2] [1, 1, 1, 2] \\ no partitions of 1 with 2 to 4 nonzero parts ? forpart(v=1,print(v),[0,5],[2,4]) ? The behavior is unspecified if X is modified inside the loop.
The library syntax is
| |
forperm(a, p, seq) | |
Evaluates seq, where the formal variable p goes through some
permutations given by a
? forperm(3, p, print(p)) Vecsmall([1, 2, 3]) Vecsmall([1, 3, 2]) Vecsmall([2, 1, 3]) Vecsmall([2, 3, 1]) Vecsmall([3, 1, 2]) Vecsmall([3, 2, 1])
When a is itself a
? forperm([2,1,1,3], p, print(p)) Vecsmall([2, 1, 1, 3]) Vecsmall([2, 1, 3, 1]) Vecsmall([2, 3, 1, 1]) Vecsmall([3, 1, 1, 2]) Vecsmall([3, 1, 2, 1]) Vecsmall([3, 2, 1, 1])
| |
forprime(p = a, {b}, seq) | |
Evaluates seq,
where the formal variable p ranges over the prime numbers between the real
numbers a to b, including a and b if they are prime. More precisely,
the value of
p is incremented to
? forprime(p = 4, 10, print(p)) 5 7
Setting b to Note that the value of p cannot be modified within seq:
? forprime(p = 2, 10, p = []) *** at top-level: forprime(p=2,10,p=[]) *** ^ — *** prime index read-only: was changed to [].
| |
forprimestep(p = a, b, q, seq) | |
Evaluates seq, where the formal variable p ranges over the prime
numbers in an arithmetic progression in [a,b]: q is either an integer
(p = a (mod q)) or an intmod
? forprimestep(p = 4, 30, 5, print(p)) 19 29 ? forprimestep(p = 4, 30, Mod(1,5), print(p)) 11
Setting b to Note that the value of p cannot be modified within seq:
? forprimestep(p = 2, 10, 3, p = []) *** at top-level: forprimestep(p=2,10,3,p=[]) *** ^ — *** prime index read-only: was changed to [].
| |
forsquarefree(N = a, b, seq) | |
Evaluates seq, where the formal variable N is [n,
? forsquarefree(N=-3,9,print(N)) [-3, [-1, 1; 3, 1]] [-2, [-1, 1; 2, 1]] [-1, Mat([-1, 1])] [1, matrix(0,2)] [2, Mat([2, 1])] [3, Mat([3, 1])] [5, Mat([5, 1])] [6, [2, 1; 3, 1]] [7, Mat([7, 1])]
This function is only implemented for |a|, |b| < 264 (232 on a 32-bit
machine). It uses a sieve and runs in time O(sqrt{b} + b-a). It should
be at least 5 times faster than regular factorization as long as the interval
length b-a is much larger than sqrt{b} and get relatively faster as
the bounds increase. The function slows down dramatically
if
? B = 10^9; ? for (N = B, B+10^6, factor(N)) time = 2,463 ms. ? forfactored (N = B, B+10^6, [n,fan] = N) time = 567 ms. ? forsquarefree (N = B, B+10^6, [n,fan] = N) time = 343 ms. ? B = 10^11; ? for (N = B, B+10^6, factor(N)) time = 8,012 ms. ? forfactored (N = B, B+10^6, [n,fan] = N) time = 1,293 ms. ? forsquarefree (N = B, B+10^6, [n,fan] = N) time = 713 ms. ? B = 10^14; ? for (N = B, B+10^6, factor(N)) time = 41,283 ms. ? forsquarefree (N = B, B+10^6, [n,fan] = N) time = 33,399 ms.
The last timing is with the default
Note that all PARI multiplicative functions accept the
? s = 0; forsquarefree(N = 1, 10^7, s += moebius(N)*eulerphi(N)); s time = 2,003 ms. %1 = 6393738650 ? s = 0; for(N = 1, 10^7, s += moebius(N)*eulerphi(N)); s time = 18,024 ms. \\ slower, we must factor N. Twice. %2 = 6393738650 The following loops over the fundamental dicriminants less than X:
? X = 10^8; ? for(d=1,X, if (isfundamental(d),)) time = 53,387 ms. ? forfactored(d=1,X, if (isfundamental(d),)); time = 13,861 ms. ? forsquarefree(d=1,X, D = quaddisc(d); if (D <= X, )); time = 14,341 ms.
Note that in the last loop, the fundamental discriminants
D are not evaluated in order (since
? forsquarefree(d=1,X/4, D = quaddisc(d)); time = 3,642 ms. ? forsquarefree(d=X/4+1,X, if (d[1] % 4 == 1,)); time = 7,772 ms. This is the price we pay for a faster evaluation, We can run through negative fundamental discriminants in the same way:
? forfactored(d=-X,-1, if (isfundamental(d),));
| |
forstep(X = a, b, s, seq) | |
Evaluates seq, where the formal variable X goes from a to b in increments of s. Nothing is done if s > 0 and a > b or if s < 0 and a < b. The s can be * a positive real number, preferably an integer: X = a, a+s, a+2s...
* an intmod * a vector of steps [s1,...,sn] (the successive steps in ℝ* are used in the order they appear in s): X = a, a+s1, a+s1+s2, ...
? forstep(x=5, 10, 2, print(x)) 5 7 9 ? forstep(x=5, 10, Mod(1,3), print(x)) 7 10 ? forstep(x=5, 10, [1,2], print(x)) 5 6 8 9
Setting b to
| |
forsubgroup(H = G, {bound}, seq) | |
Evaluates seq for each subgroup H of the abelian group G (given in SNF form or as a vector of elementary divisors). If bound is present, and is a positive integer, restrict the output to subgroups of index less than bound. If bound is a vector containing a single positive integer B, then only subgroups of index exactly equal to B are computed The subgroups are not ordered in any obvious way, unless G is a p-group in which case Birkhoff's algorithm produces them by decreasing index. A subgroup is given as a matrix whose columns give its generators on the implicit generators of G. For example, the following prints all subgroups of index less than 2 in G = ℤ/2ℤ g1 x ℤ/2ℤ g2:
? G = [2,2]; forsubgroup(H=G, 2, print(H)) [1; 1] [1; 2] [2; 1] [1, 0; 1, 1]
The last one, for instance is generated by (g1, g1 + g2). This
routine is intended to treat huge groups, when For maximal speed the subgroups have been left as produced by the algorithm. To print them in canonical form (as left divisors of G in HNF form), one can for instance use
? G = matdiagonal([2,2]); forsubgroup(H=G, 2, print(mathnf(concat(G,H)))) [2, 1; 0, 1] [1, 0; 0, 2] [2, 0; 0, 1] [1, 0; 0, 1]
Note that in this last representation, the index [G:H] is given by the
determinant. See
The library syntax is
| |
forsubset(nk, s, seq) | |
If nk is a nonnegative integer n, evaluates
? forsubset([5,3], s, print(s)) Vecsmall([1, 2, 3]) Vecsmall([1, 2, 4]) Vecsmall([1, 2, 5]) Vecsmall([1, 3, 4]) Vecsmall([1, 3, 5]) Vecsmall([1, 4, 5]) Vecsmall([2, 3, 4]) Vecsmall([2, 3, 5]) Vecsmall([2, 4, 5]) Vecsmall([3, 4, 5])
? forsubset(3, s, print(s)) Vecsmall([]) Vecsmall([1]) Vecsmall([2]) Vecsmall([3]) Vecsmall([1, 2]) Vecsmall([1, 3]) Vecsmall([2, 3]) Vecsmall([1, 2, 3])
The running time is proportional to the number
of subsets enumerated, respectively 2n and
? c = 0; forsubset([40,35],s,c++); c time = 128 ms. %4 = 658008 ? binomial(40,35) %5 = 658008
| |
forvec(X = v, seq, {flag = 0}) | |
Let v be an n-component vector (where n is arbitrary) of
two-component vectors [ai,bi] for 1 ≤ i ≤ n, where all entries
ai, bi are real numbers.
This routine lets X vary over the n-dimensional
box given by v with unit steps: X is an n-dimensional vector whose i-th
entry X[i] runs through ai, ai+1, ai+2,... stopping with the
first value greater than bi (note that neither ai nor
bi - ai are required to be integers). The values of X are ordered
lexicographically, like embedded If flag = 1, generate only nondecreasing vectors X, and if flag = 2, generate only strictly increasing vectors X.
? forvec (X=[[0,1],[-1,1]], print(X)); [0, -1] [0, 0] [0, 1] [1, -1] [1, 0] [1, 1] ? forvec (X=[[0,1],[-1,1]], print(X), 1); [0, 0] [0, 1] [1, 1] ? forvec (X=[[0,1],[-1,1]], print(X), 2) [0, 1] As a shortcut, a vector of the form v = [[0,c1-1],...[0,cn-1]] can be abbreviated as v = [c1,...cn] and flag is ignored in this case. More generally, if v is a vector of nonnegative integers ci the loop runs over representatives of ℤn/vℤn; and flag is again ignored. The vector v may contain zero entries, in which case the loop spans an infinite lattice. The values are ordered lexicographically, graded by increasing L1-norm on free (ci = 0) components.
This allows to iterate over elements of abelian groups using their
? forvec (X=[2,3], print(X)); [0, 0] [0, 1] [0, 2] [1, 0] [1, 1] [1, 2] ? my(i);forvec (X=[0,0], print(X); if (i++ > 10, break)); [0, 0] [-1, 0] [0, -1] [0, 1] [1, 0] [-2, 0] [-1, -1] [-1, 1] [0, -2] [0, 2] [1, -1] ? zn = znstar(36,1); ? forvec (chi = zn.cyc, if (chareval(zn,chi,5) == 5/6, print(chi))); [1, 0] [1, 1] ? bnrchar(zn, [5], [5/6]) \\ much more efficient in general %5 = [[1, 1], [1, 0]]
| |
if(a, {seq1}, {seq2}) | |
Evaluates the expression sequence seq1 if a is nonzero, otherwise the expression seq2. Of course, seq1 or seq2 may be empty:
The value of an
x = if(n % 4 == 1, y, z); sets x to y if n is 1 modulo 4, and to z otherwise.
Successive 'else' blocks can be abbreviated in a single compound
if (test1, seq1, test2, seq2, ... testn, seqn, seqdefault); is equivalent to
if (test1, seq1 , if (test2, seq2 , ... if (testn, seqn, seqdefault)...)); For instance, this allows to write traditional switch / case constructions:
if (x == 0, do0(), x == 1, do1(), x == 2, do2(), dodefault());
Remark.
The boolean operators
if (x != 0 && f(1/x), ...) is a perfectly safe statement.
Remark. Functions such as
| |
iferr(seq1, E, seq2, {pred}) | |
Evaluates the expression sequence seq1. If an error occurs,
set the formal parameter E set to the error data.
If pred is not present or evaluates to true, catch the error
and evaluate seq2. Both pred and seq2 can reference E.
The error type is given by
? ecm(N, B = 1000!, nb = 100)= { for(a = 1, nb, iferr(ellmul(ellinit([a,1]*Mod(1,N)), [0,1]*Mod(1,N), B), E, return(gcd(lift(component(E,2)),N)), errname(E)=="e_INV" && type(component(E,2)) == "t_INTMOD")) } ? ecm(2^101-1) %2 = 7432339208719
The return value of Internal errors, "system" errors.
*
*
*
*
* Syntax errors, type errors.
*
*
*
*
*
*
*
* Overflows.
*
*
*
*
*
*
* Errors triggered intentionally.
*
* Mathematical errors.
*
*
*
*
*
*
nfalgtobasis(nfinit(t^3-2), Mod(t,t^2+1))
E has three component, 1 (
*
*
*
| |
next({n = 1}) | |
Interrupts execution of current seq, resume the next iteration of the innermost enclosing loop, within the current function call (or top level loop). If n is specified, resume at the n-th enclosing loop. If n is bigger than the number of enclosing loops, all enclosing loops are exited.
| |
return({x = 0}) | |
Returns from current subroutine, with
result x. If x is omitted, return the
| |
setdebug({D}, {n}) | |
Sets debug level for domain D to n (0 ≤ n ≤ 20).
The domain D is a character string describing a Pari feature or code
module, such as
? setdebug()[,1] \\ list of all domains ["alg", "arith", "bern", "bnf", "bnr", "bnrclassfield", ..., "zetamult"] ? \g 1 \\ sets all debug levels to 1 debug = 1 ? setdebug("bnf", 0); \\ kills messages related to bnfinit and bnfisrincipal
| |
until(a, seq) | |
Evaluates seq until a is not
equal to 0 (i.e. until a is true). If a is initially not equal to 0,
seq is evaluated once (more generally, the condition on a is tested
after execution of the seq, not before as in
| |
while(a, seq) | |
While a is nonzero, evaluates the expression sequence seq. The test is made before evaluating the seq, hence in particular if a is initially equal to zero the seq will not be evaluated at all.
| |