Arithmetic functions

These functions are by definition functions whose natural domain of definition is either ℤ (or ℤ > 0). The way these functions are used is completely different from transcendental functions in that there are no automatic type conversions: in general only integers are accepted as arguments. An integer argument N can be given in the following alternate formats:

* t_MAT: its factorization fa = factor(N),

* t_VEC: a pair [N, fa] giving both the integer and its factorization.

This allows to compute different arithmetic functions at a given N while factoring the latter only once.

    ? N = 10!; faN = factor(N);
    ? eulerphi(N)
    %2 = 829440
    ? eulerphi(faN)
    %3 = 829440
    ? eulerphi(S = [N, faN])
    %4 = 829440
    ? sigma(S)
    %5 = 15334088


Arithmetic functions and the factoring engine

All arithmetic functions in the narrow sense of the word  — Euler's totient function, the Moebius function, the sums over divisors or powers of divisors etc. — call, after trial division by small primes, the same versatile factoring machinery described under factorint. It includes Shanks SQUFOF, Pollard Rho, ECM and MPQS stages, and has an early exit option for the functions moebius and (the integer function underlying) issquarefree. This machinery relies on a fairly strong probabilistic primality test, see ispseudoprime, but you may also set

    default(factor_proven, 1)

to ensure that all tentative factorizations are fully proven. This should not slow down PARI too much, unless prime numbers with hundreds of decimal digits occur frequently in your application.


Orders in finite groups and Discrete Logarithm functions

The following functions compute the order of an element in a finite group: ellorder (the rational points on an elliptic curve defined over a finite field), fforder (the multiplicative group of a finite field), znorder (the invertible elements in ℤ/nℤ). The following functions compute discrete logarithms in the same groups (whenever this is meaningful) elllog, fflog, znlog.

All such functions allow an optional argument specifying an integer N, representing the order of the group. (The order functions also allows any non-zero multiple of the order, with a minor loss of efficiency.) That optional argument follows the same format as given above:

* t_INT: the integer N,

* t_MAT: the factorization fa = factor(N),

* t_VEC: this is the preferred format and provides both the integer N and its factorization in a two-component vector [N, fa].

When the group is fixed and many orders or discrete logarithms will be computed, it is much more efficient to initialize this data once and for all and pass it to the relevant functions, as in

  ? p = nextprime(10^40);
  ? v = [p-1, factor(p-1)]; \\ data for discrete log & order computations
  ? znorder(Mod(2,p), v)
  %3 = 500000000000000000000000000028
  ? g = znprimroot(p);
  ? znlog(2, g, v)
  %5 = 543038070904014908801878611374


Dirichlet characters

The finite abelian group G = (ℤ/Nℤ)* can be written G = ⨁ i ≤ n (ℤ/diℤ) gi, with dn | ... | d2 | d1 (SNF condition), all di > 0, and ∏i di = φ(N).

The SNF condition makes the di unique, but the generators gi, of respective order di, are definitely not unique. The ⨁ notation means that all elements of G can be written uniquely as ∏i gini where ni ∈ ℤ/diℤ. The gi are the so-called SNF generators of G.

* a character on the abelian group ⨁ (ℤ/djℤ) gj is given by a row vector χ = [a1,...,an] of integers 0 ≤ ai < di such that χ(gj) = e(aj / dj) for all j, with the standard notation e(x) := exp(2iπ x). In other words, χ(∏ gjnj) = e(∑ aj nj / dj).

This will be generalized to more general abelian groups in later sections (Hecke characters), but in the present case of (ℤ/Nℤ)*, there is a useful alternate convention : namely, it is not necessary to impose the SNF condition and we can use Chinese reminders instead. If N = ∏ pep is the factorization of N into primes, the so-called Conrey generators of G are the generators of the (ℤ/pepℤ)* lifted to (ℤ/Nℤ)* by requesting that they be congruent to 1 modulo N/pep (for p odd we take the smallest positive primitive root mod p^2, and for p = 2 we take -1 if e2 > 1 and additionally 5 if e2 > 2). We can again write G = ⨁ i ≤ n (ℤ/Diℤ) Gi, where again ∏i Di = φ(N). These generators don't satisfy the SNF condition in general since their orders are now (p-1)pep-1 for p odd; for p = 2, the generator -1 has order 2 and 5 has order 2e2-2 (e2 > 2). Nevertheless, any m ∈ (ℤ/Nℤ)* can be uniquely decomposed as ∏ Gimi for some mi modulo Di and we can define a character by χ(Gj) = e(mj / Dj) for all j.

* The column vector of the mj, 0 ≤ mj < Dj is called the Conrey logarithm of m (discrete logarithm in terms of the Conrey generators). Note that discrete logarithms in PARI/GP are always expressed as t_COLs.

* The attached character is called the Conrey character attached to m.

To sum up a Dirichlet character can be defined by a t_INT (the Conrey label m), a t_COL (the Conrey logarithm of m, in terms of the Conrey generators) or a t_VEC (in terms of the SNF generators). The t_COL format, i.e. Conrey logarithms, is the preferred (fastest) representation.

Concretely, this works as follows:

G = znstar(N, 1) initializes (ℤ/Nℤ)*, which must be given as first arguments to all functions handling Dirichlet characters.

znconreychar transforms t_INT and t_COL to a SNF character.

znconreylog transforms t_INT and t_VEC to a Conrey logarithm.

znconreyexp transforms t_VEC and t_COL to a Conrey label.

Also available are charconj, chardiv, charmul, charker, chareval, charorder, zncharinduce, znconreyconductor (also computes the primitive character attached to the input character). The prefix char indicates that the function applies to all characters, the prefix znchar that it is specific to Dirichlet characters (on (ℤ/Nℤ)*) and the prefix znconrey that it is specific to Conrey representation.


addprimes({x = []})

Adds the integers contained in the vector x (or the single integer x) to a special table of "user-defined primes", and returns that table. Whenever factor is subsequently called, it will trial divide by the elements in this table. If x is empty or omitted, just returns the current list of extra primes.

The entries in x must be primes: there is no internal check, even if the factor_proven default is set. To remove primes from the list use removeprimes.

The library syntax is GEN addprimes(GEN x = NULL).


bestappr(x, {B})

Using variants of the extended Euclidean algorithm, returns a rational approximation a/b to x, whose denominator is limited by B, if present. If B is omitted, return the best approximation affordable given the input accuracy; if you are looking for true rational numbers, presumably approximated to sufficient accuracy, you should first try that option. Otherwise, B must be a positive real scalar (impose 0 < b ≤ B).

* If x is a t_REAL or a t_FRAC, this function uses continued fractions.

  ? bestappr(Pi, 100)
  %1 = 22/7
  ? bestappr(0.1428571428571428571428571429)
  %2 = 1/7
  ? bestappr([Pi, sqrt(2) + 'x], 10^3)
  %3 = [355/113, x + 1393/985]

By definition, a/b is the best rational approximation to x if |b x - a| < |v x - u| for all integers (u,v) with 0 < v ≤ B. (Which implies that n/d is a convergent of the continued fraction of x.)

* If x is a t_INTMOD modulo N or a t_PADIC of precision N = p^k, this function performs rational modular reconstruction modulo N. The routine then returns the unique rational number a/b in coprime integers |a| < N/2B and b ≤ B which is congruent to x modulo N. Omitting B amounts to choosing it of the order of sqrt{N/2}. If rational reconstruction is not possible (no suitable a/b exists), returns [].

  ? bestappr(Mod(18526731858, 11^10))
  %1 = 1/7
  ? bestappr(Mod(18526731858, 11^20))
  %2 = []
  ? bestappr(3 + 5 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + O(5^7))
  %2 = -1/3

In most concrete uses, B is a prime power and we performed Hensel lifting to obtain x.

The function applies recursively to components of complex objects (polynomials, vectors,...). If rational reconstruction fails for even a single entry, return [].

The library syntax is GEN bestappr(GEN x, GEN B = NULL).


bestapprPade(x, {B})

Using variants of the extended Euclidean algorithm, returns a rational function approximation a/b to x, whose denominator is limited by B, if present. If B is omitted, return the best approximation affordable given the input accuracy; if you are looking for true rational functions, presumably approximated to sufficient accuracy, you should first try that option. Otherwise, B must be a non-negative real (impose 0 ≤ degree(b) ≤ B).

* If x is a t_POLMOD modulo N this function performs rational modular reconstruction modulo N. The routine then returns the unique rational function a/b in coprime polynomials, with degree(b) ≤ B and degree(a) minimal, which is congruent to x modulo N. Omitting B amounts to choosing it equal to the floor of degree(N) / 2. If rational reconstruction is not possible (no suitable a/b exists), returns [].

  ? T = Mod(x^3 + x^2 + x + 3, x^4 - 2);
  ? bestapprPade(T)
  %2 = (2*x - 1)/(x - 1)
  ? U = Mod(1 + x + x^2 + x^3 + x^5, x^9);
  ? bestapprPade(U)  \\ internally chooses B = 4
  %3 = []
  ? bestapprPade(U, 5) \\ with B = 5, a solution exists
  %4 = (2*x^4 + x^3 - x - 1)/(-x^5 + x^3 + x^2 - 1)

* If x is a t_RFRAC or t_SER, this function implicitly converts the input to t_POLMOD modulo N = t^k fractions.

  ? T = 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + O(x^7);
  ? bestapprPade(T)
  %1 = 1/(-x + 1)

The function applies recursively to components of complex objects (polynomials, vectors,...). If rational reconstruction fails for even a single entry, return [].

The library syntax is GEN bestapprPade(GEN x, long B).


bezout(x, y)

Deprecated alias for gcdext

The library syntax is GEN gcdext0(GEN x, GEN y).


bigomega(x)

Number of prime divisors of the integer |x| counted with multiplicity:

  ? factor(392)
  %1 =
  [2 3]
  
  [7 2]
  
  ? bigomega(392)
  %2 = 5;  \\ = 3+2
  ? omega(392)
  %3 = 2;  \\ without multiplicity

The library syntax is long bigomega(GEN x).


charconj(cyc, chi)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector χ = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

This function returns the conjugate character.

  ? cyc = [15,5]; chi = [1,1];
  ? charconj(cyc, chi)
  %2 = [14, 4]
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? charconj(bnf, [1])
  %5 = [2]

For Dirichlet characters (when cyc is znstar(q,1)), characters in Conrey representation are available, see Section se:dirichletchar or ??character:

  ? G = znstar(8, 1);  \\ (Z/8Z)*
  ? charorder(G, 3)  \\ Conrey label
  %2 = 2
  ? chi = znconreylog(G, 3);
  ? charorder(G, chi)  \\ Conrey logarithm
  %4 = 2

The library syntax is GEN charconj0(GEN cyc, GEN chi). Also available is GEN charconj(GEN cyc, GEN chi), when cyc is known to be a vector of elementary divisors and chi a compatible character (no checks).


chardiv(cyc, a, b)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector a = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

Given two characters a and b, return the character a / b = a b.

  ? cyc = [15,5]; a = [1,1]; b =  [2,4];
  ? chardiv(cyc, a,b)
  %2 = [14, 2]
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? chardiv(bnf, [1], [2])
  %5 = [2]

For Dirichlet characters on (ℤ/Nℤ)*, additional representations are available (Conrey labels, Conrey logarithm), see Section se:dirichletchar or ??character. If the two characters are in the same format, the result is given in the same format, otherwise a Conrey logarithm is used.

  ? G = znstar(100, 1);
  ? G.cyc
  %2 = [20, 2]
  ? a = [10, 1]; \\ usual representation for characters
  ? b = 7; \\ Conrey label;
  ? c = znconreylog(G, 11); \\ Conrey log
  ? chardiv(G, b,b)
  %6 = 1   \\ Conrey label
  ? chardiv(G, a,b)
  %7 = [0, 5]~  \\ Conrey log
  ? chardiv(G, a,c)
  %7 = [0, 14]~   \\ Conrey log

The library syntax is GEN chardiv0(GEN cyc, GEN a, GEN b). Also available is GEN chardiv(GEN cyc, GEN a, GEN b), when cyc is known to be a vector of elementary divisors and a, b are compatible characters (no checks).


chareval(G, chi, x, {z})

Let G be an abelian group structure affording a discrete logarithm method, e.g G = znstar(N, 1) for (ℤ/Nℤ)* or a bnr structure, let x be an element of G and let chi be a character of G (see the note below for details). This function returns the value of chi at x.

Note on characters. Let K be some field. If G is an abelian group, let χ: G → K* be a character of finite order and let o be a multiple of the character order such that χ(n) = ζc(n) for some fixed ζ ∈ K* of multiplicative order o and a unique morphism c: G → (ℤ/oℤ,+). Our usual convention is to write G = (ℤ/o1ℤ) g1 ⨁ ...⨁ (ℤ/odℤ) gd for some generators (gi) of respective order di, where the group has exponent o := lcmi oi. Since ζ^o = 1, the vector (ci) in ∏ (ℤ/oiℤ) defines a character χ on G via χ(gi) = ζci (o/oi) for all i. Classical Dirichlet characters have values in K = ℂ and we can take ζ = exp(2iπ/o).

Note on Dirichlet characters. In the special case where bid is attached to G = (ℤ/qℤ)* (as per G = znstar(q,1)), the Dirichlet character chi can be written in one of the usual 3 formats: a t_VEC in terms of bid.gen as above, a t_COL in terms of the Conrey generators, or a t_INT (Conrey label); see Section se:dirichletchar or ??character.

The character value is encoded as follows, depending on the optional argument z:

* If z is omitted: return the rational number c(x)/o for x coprime to q, where we normalize 0 ≤ c(x) < o. If x can not be mapped to the group (e.g. x is not coprime to the conductor of a Dirichlet or Hecke character) we return the sentinel value -1.

* If z is an integer o, then we assume that o is a multiple of the character order and we return the integer c(x) when x belongs to the group, and the sentinel value -1 otherwise.

* z can be of the form [zeta, o], where zeta is an o-th root of 1 and o is a multiple of the character order. We return ζc(x) if x belongs to the group, and the sentinel value 0 otherwise. (Note that this coincides with the usual extension of Dirichlet characters to ℤ, or of Hecke characters to general ideals.)

* Finally, z can be of the form [vzeta, o], where vzeta is a vector of powers ζ^0,..., ζo-1 of some o-th root of 1 and o is a multiple of the character order. As above, we return ζc(x) after a table lookup. Or the sentinel value 0.

The library syntax is GEN chareval(GEN G, GEN chi, GEN x, GEN z = NULL).


chargalois(cyc, {ORD})

Let cyc represent a finite abelian group by its elementary divisors (any object which has a .cyc method is also allowed, i.e. the output of znstar or bnrinit). Return a list of representatives for the Galois orbits of complex characters of G. If ORD is present, select characters depending on their orders:

* if ORD is a t_INT, restrict to orders less than this bound;

* if ORD is a t_VEC or t_VECSMALL, restrict to orders in the list.

  ? G = znstar(96);
  ? #chargalois(G) \\ 16 orbits of characters mod 96
  %2 = 16
  ? #chargalois(G,4) \\ order less than 4
  %3 = 12
  ? chargalois(G,[1,4]) \\ order 1 or 4; 5 orbits
  %4 = [[0, 0, 0], [2, 0, 0], [2, 1, 0], [2, 0, 1], [2, 1, 1]]

Given a character χ, of order n (charorder(G,chi)), the elements in its orbit are the φ(n) characters χ^i, (i,n) = 1.

The library syntax is GEN chargalois(GEN cyc, GEN ORD = NULL).


charker(cyc, chi)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector χ = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

This function returns the kernel of χ, as a matrix K in HNF which is a left-divisor of matdiagonal(d). Its columns express in terms of the gj the generators of the subgroup. The determinant of K is the kernel index.

  ? cyc = [15,5]; chi = [1,1];
  ? charker(cyc, chi)
  %2 =
  [15 12]
  
  [ 0  1]
  
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? charker(bnf, [1])
  %5 =
  [3]

Note that for Dirichlet characters (when cyc is znstar(q, 1)), characters in Conrey representation are available, see Section se:dirichletchar or ??character.

  ? G = znstar(8, 1);  \\ (Z/8Z)*
  ? charker(G, 1) \\ Conrey label for trivial character
  %2 =
  [1 0]
  
  [0 1]

The library syntax is GEN charker0(GEN cyc, GEN chi). Also available is GEN charker(GEN cyc, GEN chi), when cyc is known to be a vector of elementary divisors and chi a compatible character (no checks).


charmul(cyc, a, b)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector a = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

Given two characters a and b, return the product character ab.

  ? cyc = [15,5]; a = [1,1]; b =  [2,4];
  ? charmul(cyc, a,b)
  %2 = [3, 0]
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? charmul(bnf, [1], [2])
  %5 = [0]

For Dirichlet characters on (ℤ/Nℤ)*, additional representations are available (Conrey labels, Conrey logarithm), see Section se:dirichletchar or ??character. If the two characters are in the same format, their product is given in the same format, otherwise a Conrey logarithm is used.

  ? G = znstar(100, 1);
  ? G.cyc
  %2 = [20, 2]
  ? a = [10, 1]; \\ usual representation for characters
  ? b = 7; \\ Conrey label;
  ? c = znconreylog(G, 11); \\ Conrey log
  ? charmul(G, b,b)
  %6 = 49   \\ Conrey label
  ? charmul(G, a,b)
  %7 = [0, 15]~  \\ Conrey log
  ? charmul(G, a,c)
  %7 = [0, 6]~   \\ Conrey log

The library syntax is GEN charmul0(GEN cyc, GEN a, GEN b). Also available is GEN charmul(GEN cyc, GEN a, GEN b), when cyc is known to be a vector of elementary divisors and a, b are compatible characters (no checks).


charorder(cyc, chi)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector χ = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

This function returns the order of the character chi.

  ? cyc = [15,5]; chi = [1,1];
  ? charorder(cyc, chi)
  %2 = 15
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? charorder(bnf, [1])
  %5 = 3

For Dirichlet characters (when cyc is znstar(q, 1)), characters in Conrey representation are available, see Section se:dirichletchar or ??character:

  ? G = znstar(100, 1); \\ (Z/100Z)*
  ? charorder(G, 7)   \\ Conrey label
  %2 = 4

The library syntax is GEN charorder0(GEN cyc, GEN chi). Also available is GEN charorder(GEN cyc, GEN chi), when cyc is known to be a vector of elementary divisors and chi a compatible character (no checks).


charpow(cyc, a, n)

Let cyc represent a finite abelian group by its elementary divisors, i.e. (dj) represents ∑j ≤ k ℤ/djℤ with dk | ... | d1; any object which has a .cyc method is also allowed, e.g. the output of znstar or bnrinit. A character on this group is given by a row vector a = [a1,...,an] such that χ(∏ gjnj) = exp(2π i∑ aj nj / dj), where gj denotes the generator (of order dj) of the j-th cyclic component.

Given n ∈ ℤ and a character a, return the character a^n.

  ? cyc = [15,5]; a = [1,1];
  ? charpow(cyc, a, 3)
  %2 = [3, 3]
  ? charpow(cyc, a, 5)
  %2 = [5, 0]
  ? bnf = bnfinit(x^2+23);
  ? bnf.cyc
  %4 = [3]
  ? charpow(bnf, [1], 3)
  %5 = [0]

For Dirichlet characters on (ℤ/Nℤ)*, additional representations are available (Conrey labels, Conrey logarithm), see Section se:dirichletchar or ??character and the output uses the same format as the input.

  ? G = znstar(100, 1);
  ? G.cyc
  %2 = [20, 2]
  ? a = [10, 1]; \\ standard representation for characters
  ? b = 7; \\ Conrey label;
  ? c = znconreylog(G, 11); \\ Conrey log
  ? charpow(G, a,3)
  %6 = [10, 1]   \\ standard representation
  ? charpow(G, b,3)
  %7 = 43   \\ Conrey label
  ? charpow(G, c,3)
  %8 = [1, 8]~  \\ Conrey log

The library syntax is GEN charpow0(GEN cyc, GEN a, GEN n). Also available is GEN charpow(GEN cyc, GEN a, GEN n), when cyc is known to be a vector of elementary divisors (no check).


chinese(x, {y})

If x and y are both intmods or both polmods, creates (with the same type) a z in the same residue class as x and in the same residue class as y, if it is possible.

  ? chinese(Mod(1,2), Mod(2,3))
  %1 = Mod(5, 6)
  ? chinese(Mod(x,x^2-1), Mod(x+1,x^2+1))
  %2 = Mod(-1/2*x^2 + x + 1/2, x^4 - 1)

This function also allows vector and matrix arguments, in which case the operation is recursively applied to each component of the vector or matrix.

  ? chinese([Mod(1,2),Mod(1,3)], [Mod(1,5),Mod(2,7)])
  %3 = [Mod(1, 10), Mod(16, 21)]

For polynomial arguments in the same variable, the function is applied to each coefficient; if the polynomials have different degrees, the high degree terms are copied verbatim in the result, as if the missing high degree terms in the polynomial of lowest degree had been Mod(0,1). Since the latter behavior is usually not the desired one, we propose to convert the polynomials to vectors of the same length first:

   ? P = x+1; Q = x^2+2*x+1;
   ? chinese(P*Mod(1,2), Q*Mod(1,3))
   %4 = Mod(1, 3)*x^2 + Mod(5, 6)*x + Mod(3, 6)
   ? chinese(Vec(P,3)*Mod(1,2), Vec(Q,3)*Mod(1,3))
   %5 = [Mod(1, 6), Mod(5, 6), Mod(4, 6)]
   ? Pol(%)
   %6 = Mod(1, 6)*x^2 + Mod(5, 6)*x + Mod(4, 6)

If y is omitted, and x is a vector, chinese is applied recursively to the components of x, yielding a residue belonging to the same class as all components of x.

Finally chinese(x,x) = x regardless of the type of x; this allows vector arguments to contain other data, so long as they are identical in both vectors.

The library syntax is GEN chinese(GEN x, GEN y = NULL). GEN chinese1(GEN x) is also available.


content(x, {D})

Computes the gcd of all the coefficients of x, when this gcd makes sense. This is the natural definition if x is a polynomial (and by extension a power series) or a vector/matrix. This is in general a weaker notion than the ideal generated by the coefficients:

  ? content(2*x+y)
  %1 = 1            \\ = gcd(2,y) over Q[y]

If x is a scalar, this simply returns the absolute value of x if x is rational (t_INT or t_FRAC), and either 1 (inexact input) or x (exact input) otherwise; the result should be identical to gcd(x, 0).

The content of a rational function is the ratio of the contents of the numerator and the denominator. In recursive structures, if a matrix or vector coefficient x appears, the gcd is taken not with x, but with its content:

  ? content([ [2], 4*matid(3) ])
  %1 = 2

The content of a t_VECSMALL is computed assuming the entries are signed integers.

The optional argument D allows to control over which ring we compute and get a more predictable behaviour:

* 1: we only consider the underlying ℚ-structure and the denominator is a (positive) rational number

* a simple variable, say 'x: all entries are considered as rational functions in K(x) for some field K and the content is an element of K.

  ? f = x + 1/y + 1/2;
  ? content(f) \\ as a t_POL in x
  %2 = 1/(2*y)
  ? content(f, 1) \\ Q-content
  %3 = 1/2
  ? content(f, y) \\ as a rational function in y
  %4 = 1/2
  ? g = x^2*y + y^2*x;
  ? content(g, x)
  %6 = y
  ? content(g, y)
  %7 = x

The library syntax is GEN content0(GEN x, GEN D = NULL).


contfrac(x, {b}, {nmax})

Returns the row vector whose components are the partial quotients of the continued fraction expansion of x. In other words, a result [a0,...,an] means that x ~ a0+1/(a1+...+1/an). The output is normalized so that an != 1 (unless we also have n = 0).

The number of partial quotients n+1 is limited by nmax. If nmax is omitted, the expansion stops at the last significant partial quotient.

  ? \p19
    realprecision = 19 significant digits
  ? contfrac(Pi)
  %1 = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2]
  ? contfrac(Pi,, 3)  \\ n = 2
  %2 = [3, 7, 15]

x can also be a rational function or a power series.

If a vector b is supplied, the numerators are equal to the coefficients of b, instead of all equal to 1 as above; more precisely, x ~ (1/b0)(a0+b1/(a1+...+bn/an)); for a numerical continued fraction (x real), the ai are integers, as large as possible; if x is a rational function, they are polynomials with deg ai = deg bi + 1. The length of the result is then equal to the length of b, unless the next partial quotient cannot be reliably computed, in which case the expansion stops. This happens when a partial remainder is equal to zero (or too small compared to the available significant digits for x a t_REAL).

A direct implementation of the numerical continued fraction contfrac(x,b) described above would be

  \\ "greedy" generalized continued fraction
  cf(x, b) =
  { my( a= vector(#b), t );
  
    x *= b[1];
    for (i = 1, #b,
      a[i] = floor(x);
      t = x - a[i]; if (!t || i == #b, break);
      x = b[i+1] / t;
    ); a;
  }

There is some degree of freedom when choosing the ai; the program above can easily be modified to derive variants of the standard algorithm. In the same vein, although no builtin function implements the related Engel expansion (a special kind of Egyptian fraction decomposition: x = 1/a1 + 1/(a_1a2) +... ), it can be obtained as follows:

  \\ n terms of the Engel expansion of x
  engel(x, n = 10) =
  { my( u = x, a = vector(n) );
    for (k = 1, n,
      a[k] = ceil(1/u);
      u = u*a[k] - 1;
      if (!u, break);
    ); a
  }

Obsolete hack. (don't use this): if b is an integer, nmax is ignored and the command is understood as contfrac(x,, b).

The library syntax is GEN contfrac0(GEN x, GEN b = NULL, long nmax). Also available are GEN gboundcf(GEN x, long nmax), GEN gcf(GEN x) and GEN gcf2(GEN b, GEN x).


contfracpnqn(x, {n = -1})

When x is a vector or a one-row matrix, x is considered as the list of partial quotients [a0,a1,...,an] of a rational number, and the result is the 2 by 2 matrix [pn,pn-1;qn,qn-1] in the standard notation of continued fractions, so pn/qn = a0+1/(a1+...+1/an). If x is a matrix with two rows [b0,b1,...,bn] and [a0,a1,...,an], this is then considered as a generalized continued fraction and we have similarly pn/qn = (1/b0)(a0+b1/(a1+...+bn/an)). Note that in this case one usually has b0 = 1.

If n ≥ 0 is present, returns all convergents from p0/q0 up to pn/qn. (All convergents if x is too small to compute the n+1 requested convergents.)

  ? a = contfrac(Pi,10)
  %1 = [3, 7, 15, 1, 292, 1, 1, 1, 3]
  ? allpnqn(x) = contfracpnqn(x,#x) \\ all convergents
  ? allpnqn(a)
  %3 =
  [3 22 333 355 103993 104348 208341 312689 1146408]
  
  [1  7 106 113  33102  33215  66317  99532  364913]
  ? contfracpnqn(a) \\ last two convergents
  %4 =
  [1146408 312689]
  
  [ 364913  99532]
  
  ? contfracpnqn(a,3) \\ first three convergents
  %5 =
  [3 22 333 355]
  
  [1  7 106 113]

The library syntax is GEN contfracpnqn(GEN x, long n). also available is GEN pnqn(GEN x) for n = -1.


core(n, {flag = 0})

If n is an integer written as n = df^2 with d squarefree, returns d. If flag is non-zero, returns the two-element row vector [d,f]. By convention, we write 0 = 0 x 1^2, so core(0, 1) returns [0,1].

The library syntax is GEN core0(GEN n, long flag). Also available are GEN core(GEN n) (flag = 0) and GEN core2(GEN n) (flag = 1)


coredisc(n, {flag = 0})

A fundamental discriminant is an integer of the form t = 1 mod 4 or 4t = 8,12 mod 16, with t squarefree (i.e. 1 or the discriminant of a quadratic number field). Given a non-zero integer n, this routine returns the (unique) fundamental discriminant d such that n = df^2, f a positive rational number. If flag is non-zero, returns the two-element row vector [d,f]. If n is congruent to 0 or 1 modulo 4, f is an integer, and a half-integer otherwise.

By convention, coredisc(0, 1)) returns [0,1].

Note that quaddisc(n) returns the same value as coredisc(n), and also works with rational inputs n ∈ ℚ*.

The library syntax is GEN coredisc0(GEN n, long flag). Also available are GEN coredisc(GEN n) (flag = 0) and GEN coredisc2(GEN n) (flag = 1)


dirdiv(x, y)

x and y being vectors of perhaps different lengths but with y[1] != 0 considered as Dirichlet series, computes the quotient of x by y, again as a vector.

The library syntax is GEN dirdiv(GEN x, GEN y).


direuler(p = a, b, expr, {c})

Computes the Dirichlet series attached to the Euler product of expression expr as p ranges through the primes from a to b. expr must be a polynomial or rational function in another variable than p (say X) and expr(X) is understood as the local factor expr(p-s).

The series is output as a vector of coefficients. If c is omitted, output the first b coefficients of the series; otherwise, output the first c coefficients. The following command computes the sigma function, attached to ζ(s)ζ(s-1):

  ? direuler(p=2, 10, 1/((1-X)*(1-p*X)))
  %1 = [1, 3, 4, 7, 6, 12, 8, 15, 13, 18]
  
  ? direuler(p=2, 10, 1/((1-X)*(1-p*X)), 5) \\ fewer terms
  %2 = [1, 3, 4, 7, 6]

Setting c < b is useless (the same effect would be achieved by setting b = c). If c > b, the computed coefficients are "missing" Euler factors:

  ? direuler(p=2, 10, 1/((1-X)*(1-p*X)), 15) \\ more terms, no longer = sigma !
  %3 = [1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 0, 28, 0, 24, 24]

The library syntax is direuler(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b)


dirmul(x, y)

x and y being vectors of perhaps different lengths representing the Dirichlet series ∑n xn n-s and ∑n yn n-s, computes the product of x by y, again as a vector.

  ? dirmul(vector(10,n,1), vector(10,n,moebius(n)))
  %1 = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]

The product length is the minimum of #x*v(y) and #y*v(x), where v(x) is the index of the first non-zero coefficient.

  ? dirmul([0,1], [0,1]);
  %2 = [0, 0, 0, 1]

The library syntax is GEN dirmul(GEN x, GEN y).


divisors(x, {flag = 0})

Creates a row vector whose components are the divisors of x. The factorization of x (as output by factor) can be used instead. If flag = 1, return pairs [d, factor(d)].

By definition, these divisors are the products of the irreducible factors of n, as produced by factor(n), raised to appropriate powers (no negative exponent may occur in the factorization). If n is an integer, they are the positive divisors, in increasing order.

  ? divisors(12)
  %1 = [1, 2, 3, 4, 6, 12]
  ? divisors(12, 1) \\ include their factorization
  %2 = [[1, matrix(0,2)], [2, Mat([2, 1])], [3, Mat([3, 1])],
        [4, Mat([2, 2])], [6, [2, 1; 3, 1]], [12, [2, 2; 3, 1]]]
  
  ? divisors(x^4 + 2*x^3 + x^2) \\ also works for polynomials
  %3 = [1, x, x^2, x + 1, x^2 + x, x^3 + x^2, x^2 + 2*x + 1,
        x^3 + 2*x^2 + x, x^4 + 2*x^3 + x^2]

The library syntax is GEN divisors0(GEN x, long flag). The functions GEN divisors(GEN N) (flag = 0) and GEN divisors_factored(GEN N) (flag = 1) are also available.


divisorslenstra(N, r, s)

Given three integers N > s > r ≥ 0 such that (r,s) = 1 and s^3 > N, find all divisors d of N such that d = r (mod s). There are at most 11 such divisors (Lenstra).

  ? N = 245784; r = 19; s = 65 ;
  ? divisorslenstra(N, r, s)
  %2 = [19, 84, 539, 1254, 3724, 245784]
  ? [ d | d <- divisors(N), d % s == r]
  %3 = [19, 84, 539, 1254, 3724, 245784]

When the preconditions are not met, the result is undefined:

  ? N = 4484075232; r = 7; s = 1303; s^3 > N
  %4 = 0
  ? divisorslenstra(N, r, s)
  ? [ d | d <- divisors(N), d % s == r ]
  %6 = [7, 2613, 9128, 19552, 264516, 3407352, 344928864]

(Divisors were missing but s^3 < N.)

The library syntax is GEN divisorslenstra(GEN N, GEN r, GEN s).


eulerphi(x)

Euler's φ (totient) function of the integer |x|, in other words |(ℤ/xℤ)*|.

  ? eulerphi(40)
  %1 = 16

According to this definition we let φ(0) := 2, since ℤ^ *= {-1,1}; this is consistent with znstar(0): we have znstar(n).no = eulerphi(n) for all n ∈ ℤ.

The library syntax is GEN eulerphi(GEN x).


factor(x, {D})

Factor x over domain D; if D is omitted, it is determined from x. For instance, if x is an integer, it is factored in ℤ, if it is a polynomial with rational coefficients, it is factored in ℚ[x], etc., see below for details. The result is a two-column matrix: the first contains the irreducibles dividing x (rational or Gaussian primes, irreducible polynomials), and the second the exponents. By convention, 0 is factored as 0^1.

x ∈ ℚ. See factorint for the algorithms used. The factorization includes the unit -1 when x < 0 and all other factors are positive; a denominator is factored with negative exponents. The factors are sorted in increasing order.

  ? factor(-7/106)
  %1 =
  [-1  1]
  
  [ 2 -1]
  
  [ 7  1]
  
  [53 -1]

By convention, 1 is factored as matrix(0,2) (the empty factorization, printed as [;]).

Large rational "primes" > 264 in the factorization are in fact pseudoprimes (see ispseudoprime), a priori not rigorously proven primes. Use isprime to prove primality of these factors, as in

  ? fa = factor(2^2^7 + 1)
  %1 =
  [59649589127497217 1]
  
  [5704689200685129054721 1]
  
  ? isprime( fa[,1] )
  %2 = [1, 1]~   \\ both entries are proven primes

Another possibility is to globally set the default factor_proven, which will perform a rigorous primality proof for each pseudoprime factor but will slow down PARI.

A t_INT argument D can be added, meaning that we look only for prime factors p < D. The limit D must be non-negative. In this case, all but the last factor are proven primes, but the remaining factor may actually be a proven composite! If the remaining factor is less than D^2, then it is prime.

  ? factor(2^2^7 +1, 10^5)
  %3 =
  [340282366920938463463374607431768211457 1]

Deprecated feature. Setting D = 0 is the same as setting it to primelimit + 1.

This routine uses trial division and perfect power tests, and should not be used for huge values of D (at most 10^9, say): factorint(, 1 + 8) will in general be faster. The latter does not guarantee that all small prime factors are found, but it also finds larger factors and in a more efficient way.

  ? F = (2^2^7 + 1) * 1009 * (10^5+3); factor(F, 10^5)  \\ fast, incomplete
  time = 0 ms.
  %4 =
  [1009 1]
  
  [34029257539194609161727850866999116450334371 1]
  
  ? factor(F, 10^9)    \\ slow
  time = 3,260 ms.
  %6 =
  [1009 1]
  
  [100003 1]
  
  [340282366920938463463374607431768211457 1]
  
  ? factorint(F, 1+8)  \\ much faster and all small primes were found
  time = 8 ms.
  %7 =
  [1009 1]
  
  [100003 1]
  
  [340282366920938463463374607431768211457 1]
  
  ? factor(F)   \\ complete factorization
  time = 60 ms.
  %8 =
  [1009 1]
  
  [100003 1]
  
  [59649589127497217 1]
  
  [5704689200685129054721 1]

Setting D = I will factor in the Gaussian integers ℤ[i]:

x ∈ ℚ(i). The factorization is performed with Gaussian primes in ℤ[i] and includes Gaussian units in {±1, ± i}; factors are sorted by increasing norm. Except for a possible leading unit, the Gaussian factors are normalized: rational factors are positive and irrational factors have positive imaginary part (a canonical represneta.

Unless factor_proven is set, large factors are actually pseudoprimes, not proven primes; a rational factor is prime if less than 264 and an irrational one if its norm is less than 264.

  ? factor(5*I)
  %1 =
  [  2 + I 1]
  
  [1 + 2*I 1]

One can force the factorization of a rational number by setting the domain D = I:

  ? factor(-5, I)
  %2 =
  [      I 1]
  
  [  2 + I 1]
  
  [1 + 2*I 1]
  ? factorback(%)
  %3 = -5

Univariate polynomials and rational functions. PARI can factor univariate polynomials in K[t]. The following base fields K are currently supported: ℚ, ℝ, ℂ, ℚp, finite fields and number fields. See factormod and factorff for the algorithms used over finite fields and nffactor for the algorithms over number fields. The irreducible factors are sorted by increasing degree and normalized: they are monic except when K = ℚ where they are primitive in ℤ[t].

The content is not included in the factorization, in particular factorback will in general recover the original x only up to multiplication by an element of K*: when K != ℚ, this scalar is pollead(x) (since irreducible factors are monic); and when K = ℚ you can either ask for the ℚ-content explicitly of use factorback:

  ? P = t^2 + 5*t/2 + 1; F = factor(P)
  %1 =
  [t + 2 1]
  
  [2*t + 1 1]
  
  ? content(P, 1) \\ Q-content
  %2 = 1/2
  
  ? pollead(factorback(F)) / pollead(P)
  %3 = 2

You can specify K using the optional "domain" argument D as follows

* K = ℚ : D a rational number (t_INT or t_FRAC),

* K = ℤ/pℤ with p prime : D a t_INTMOD modulo p; factoring modulo a non-prime number is not supported.

* K = 𝔽q : D a t_FFELT encoding the finite field; you can also use a t_POLMOD of t_INTMOD modulo a prime p but this is usualy less convenient;

* K = ℚ[X]/(T) a number field : D a t_POLMOD modulo T,

* K = ℚ(i) (alternate syntax for special case): D = I,

* K = ℚ(w) a quadratic number field (alternate syntax for special case): D a t_QUAD,

* K = ℝ : D a real number (t_REAL); truncate the factorization at accuracy precision(D). If x is inexact and precision(x) is less than precision(D), then the precision of x is used instead.

* K = ℂ : D a complex number with a t_REAL component, e.g. I * 1.; truncate the factorization as for K = ℝ,

* K = ℚp : D a t_PADIC; truncate the factorization at p-adic accuracy padicprec(D), possibly less if x is inexact with insufficient p-adic accuracy;

  ? T = x^2+1;
  ? factor(T, 1);                      \\ over Q
  ? factor(T, Mod(1,3))                \\ over F3
  ? factor(T, ffgen(ffinit(3,2,'t))^0) \\ over F3^2
  ? factor(T, Mod(Mod(1,3), t^2+t+2))  \\ over F3^2, again
  ? factor(T, O(3^6))                  \\ over Q3, precision 6
  ? factor(T, 1.)                      \\ over R, current precision
  ? factor(T, I*1.)                    \\ over C
  ? factor(T, Mod(1, y^3-2))           \\ over Q(21/3)

In most cases, it is possible and simpler to call a specialized variant rather than use the above scheme:

  ? factormod(T, 3)              \\ over F3
  ? factormod(T, [t^2+t+2, 3])   \\ over F3^2
  ? factormod(T, ffgen(3^2, 't)) \\ over F3^2
  ? factorpadic(T, 3,6)          \\ over Q3, precision 6
  ? nffactor(y^3-2, T)           \\ over Q(21/3)
  ? polroots(T)                  \\ over C
  ? polrootsreal(T)              \\ over R (real polynomial)

It is also possible to let the routine use the smallest field containing all coefficients, taking into account quotient structures induced by t_INTMODs and t_POLMODs (e.g. if a coefficient in ℤ/nℤ is known, all rational numbers encountered are first mapped to ℤ/nℤ; different moduli will produce an error):

  ? T = x^2+1;
  ? factor(T);                         \\ over Q
  ? factor(T*Mod(1,3))                 \\ over F3
  ? factor(T*ffgen(ffinit(3,2,'t))^0)  \\ over F3^2
  ? factor(T*Mod(Mod(1,3), t^2+t+2))   \\ over F3^2, again
  ? factor(T*(1 + O(3^6))              \\ over Q3, precision 6
  ? factor(T*1.)                       \\ over R, current precision
  ? factor(T*(1.+0.*I))                \\ over C
  ? factor(T*Mod(1, y^3-2))            \\ over Q(21/3)

Multiplying by a suitable field element equal to 1 ∈ K in this way is error-prone and is not recommanded. Factoring existing polynomials with obvious fields of coefficients is fine, the domain argument D should be used instead ad hoc conversions.

Note on inexact polynomials. Polynomials with inexact coefficients (e.g. floating point or p-adic numbers) are first rounded to an exact representation, then factored to (potentially) infinite accuracy and we return a truncated approximation of that virtual factorization. To avoid pitfalls, we advise to only factor exact polynomials:

  ? factor(x^2-1+O(2^2)) \\ rounded to x^2 + 3, irreducible in Q2
  %1 =
  [(1 + O(2^2))*x^2 + O(2^2)*x + (1 + 2 + O(2^2)) 1]
  
  ? factor(x^2-1+O(2^3)) \\ rounded to x^2 + 7, reducible !
  %2 =
  [  (1 + O(2^3))*x + (1 + 2 + O(2^3)) 1]
  
  [(1 + O(2^3))*x + (1 + 2^2 + O(2^3)) 1]
  
  ? factor(x^2-1, O(2^2)) \\ no ambiguity now
  %3 =
  [    (1 + O(2^2))*x + (1 + O(2^2)) 1]
  
  [(1 + O(2^2))*x + (1 + 2 + O(2^2)) 1]

Note about inseparable polynomials. Polynomials with inexact coefficients are considered to be squarefree: indeed, there exist a squarefree polynomial arbitrarily close to the input, and they cannot be distinguished at the input accuracy. This means that irreducible factors are repeated according to their apparent multiplicity. On the contrary, using a specialized function such as factorpadic with an exact rational input yields the correct multiplicity when the (now exact) input is not separable. Compare:

  ? factor(z^2 + O(5^2)))
  %1 =
  [(1 + O(5^2))*z + O(5^2) 1]
  
  [(1 + O(5^2))*z + O(5^2) 1]
  ? factor(z^2, O(5^2))
  %2 =
  [1 + O(5^2))*z + O(5^2) 2]

Multivariate polynomials and rational functions. PARI recursively factors multivariate polynomials in K[t1,..., td] for the same fields K as above and the argument D is used in the same way to specify K. The irreducible factors are sorted by their main variable (least priority first) then by increasing degree.

  ? factor(x^2 + y^2, Mod(1,5))
  %1 =
  [          x + Mod(2, 5)*y 1]
  
  [Mod(1, 5)*x + Mod(3, 5)*y 1]
  
  ? factor(x^2 + y^2, O(5^2))
  %2 =
  [  (1 + O(5^2))*x + (O(5^2)*y^2 + (2 + 5 + O(5^2))*y + O(5^2)) 1]
  
  [(1 + O(5^2))*x + (O(5^2)*y^2 + (3 + 3*5 + O(5^2))*y + O(5^2)) 1]
  
  ? lift(%)
  %3 =
  [ x + 7*y 1]
  
  [x + 18*y 1]

Note that the implementation does not really support inexact real fields (ℝ or ℂ) and usually misses factors even if the input is exact:

  ? factor(x^2 + y^2, I)  \\ over Q(i)
  %4 =
  [x - I*y 1]
  
  [x + I*y 1]
  
  ? factor(x^2 + y^2, I*1.) \\ over C
  %5 =
  [x^2 + y^2 1]

The library syntax is GEN factor0(GEN x, GEN D = NULL).

GEN factor(GEN x) GEN boundfact(GEN x, ulong lim).


factorback(f, {e})

Gives back the factored object corresponding to a factorization. The integer 1 corresponds to the empty factorization.

If e is present, e and f must be vectors of the same length (e being integral), and the corresponding factorization is the product of the f[i]e[i].

If not, and f is vector, it is understood as in the preceding case with e a vector of 1s: we return the product of the f[i]. Finally, f can be a regular factorization, as produced with any factor command. A few examples:

  ? factor(12)
  %1 =
  [2 2]
  
  [3 1]
  
  ? factorback(%)
  %2 = 12
  ? factorback([2,3], [2,1])   \\ 2^3 * 3^1
  %3 = 12
  ? factorback([5,2,3])
  %4 = 30

The library syntax is GEN factorback2(GEN f, GEN e = NULL). Also available is GEN factorback(GEN f) (case e = NULL).


factorcantor(x, p)

This function is obsolete, use factormod.

The library syntax is GEN factmod(GEN x, GEN p).


factorff(x, {p}, {a})

Obsolete, kept for backward compatibility: use factormod.

The library syntax is GEN factorff(GEN x, GEN p = NULL, GEN a = NULL).


factorial(x)

Factorial of x. The expression x! gives a result which is an integer, while factorial(x) gives a real number.

The library syntax is GEN mpfactr(long x, long prec). GEN mpfact(long x) returns x! as a t_INT.


factorint(x, {flag = 0})

Factors the integer n into a product of pseudoprimes (see ispseudoprime), using a combination of the Shanks SQUFOF and Pollard Rho method (with modifications due to Brent), Lenstra's ECM (with modifications by Montgomery), and MPQS (the latter adapted from the LiDIA code with the kind permission of the LiDIA maintainers), as well as a search for pure powers. The output is a two-column matrix as for factor: the first column contains the "prime" divisors of n, the second one contains the (positive) exponents.

By convention 0 is factored as 0^1, and 1 as the empty factorization; also the divisors are by default not proven primes if they are larger than 264, they only failed the BPSW compositeness test (see ispseudoprime). Use isprime on the result if you want to guarantee primality or set the factor_proven default to 1. Entries of the private prime tables (see addprimes) are also included as is.

This gives direct access to the integer factoring engine called by most arithmetical functions. flag is optional; its binary digits mean 1: avoid MPQS, 2: skip first stage ECM (we may still fall back to it later), 4: avoid Rho and SQUFOF, 8: don't run final ECM (as a result, a huge composite may be declared to be prime). Note that a (strong) probabilistic primality test is used; thus composites might not be detected, although no example is known.

You are invited to play with the flag settings and watch the internals at work by using gp's debug default parameter (level 3 shows just the outline, 4 turns on time keeping, 5 and above show an increasing amount of internal details).

The library syntax is GEN factorint(GEN x, long flag).


factormod(f, {D}, {flag = 0})

Factors the polynomial f over the finite field defined by the domain D as follows:

* D = p a prime: factor over 𝔽p;

* D = [T,p] for a prime p and T an irreducible polynomial over 𝔽p: factor over 𝔽p[x]/(T);

* D a t_FFELT: factor over the attached field;

* D omitted: factor over the field of definition of f, which must be a finite field.

The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix, the first column being the irreducible polynomials dividing f, and the second the exponents. By convention, the 0 polynomial factors as 0^1; a non-zero constant polynomial has empty factorization, a 0 x 2 matrix. The irreducible factors are ordered by increasing degree and the result is canonical: it will not change across multiple calls or sessions.

  ? factormod(x^2 + 1, 3)  \\ over F3
  %1 =
  [Mod(1, 3)*x^2 + Mod(1, 3) 1]
  ? liftall( factormod(x^2 + 1, [3, t^2+1]) ) \\ over F9
  %2 =
  [  x + t 1]
  
  [x + 2*t 1]
  ? t = ffgen(t^2+Mod(1,3)); factormod(x^2 + t^0) \\ same using t_FFELT
  %3 =
  [  x + t 1]
  
  [x + 2*t 1]
  ? factormod(x^2+Mod(1,3))
  %4 =
  [Mod(1, 3)*x^2 + Mod(1, 3) 1]
  ? liftall( factormod(x^2 + Mod(Mod(1,3),y^2+1)) )
  %5 =
  [  x + y 1]
  
  [x + 2*y 1]

If flag is non-zero, outputs only the degrees of the irreducible polynomials (for example to compute an L-function). By convention, a constant polynomial (including the 0 polynomial) has empty factorization. The degrees appear in increasing order but need not correspond to the ordering with flag = 0 when multiplicities are present.

  ? f = x^3 + 2*x^2 + x + 2;
  ? factormod(f, 5)  \\ (x+2)^2 * (x+3)
  %1 =
  [Mod(1, 5)*x + Mod(2, 5) 2]
  
  [Mod(1, 5)*x + Mod(3, 5) 1]
  ? factormod(f, 5, 1) \\ (deg 1) * (deg 1)^2
  %2 =
  [1 1]
  
  [1 2]

The library syntax is GEN factormod0(GEN f, GEN D = NULL, long flag).


factormodDDF(f, {D})

Distinct-degree factorization of the squarefree polynomial f over the finite field defined by the domain D as follows:

* D = p a prime: factor over 𝔽p;

* D = [T,p] for a prime p and T an irreducible polynomial over 𝔽p: factor over 𝔽p[x]/(T);

* D a t_FFELT: factor over the attached field;

* D omitted: factor over the field of definition of f, which must be a finite field.

This is somewhat faster than full factorization. The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix:

* the first column contains monic (squarefree) pairwise coprime polynomials dividing f, all of whose irreducible factors have degree d;

* the second column contains the degrees of the irreducible factors.

The factors are ordered by increasing degree and the result is canonical: it will not change across multiple calls or sessions.

  ? f = (x^2 + 1) * (x^2-1);
  ? factormodSQF(f,3) \\ squarefree over F3
  %2 =
  [Mod(1, 3)*x^4 + Mod(2, 3) 1]
  
  ? factormodDDF(f, 3)
  %3 =
  [Mod(1, 3)*x^2 + Mod(2, 3) 1]  \\ two degree 1 factors
  
  [Mod(1, 3)*x^2 + Mod(1, 3) 2]  \\ irred of degree 2
  
  ? for(i=1,10^5,factormodDDF(f,3))
  time = 424 ms.
  ? for(i=1,10^5,factormod(f,3))  \\ full factorization is slower
  time = 464 ms.
  
  ? liftall( factormodDDF(x^2 + 1, [3, t^2+1]) ) \\ over F9
  %6 =
  [x^2 + 1 1] \\ product of two degree 1 factors
  
  ? t = ffgen(t^2+Mod(1,3)); factormodDDF(x^2 + t^0) \\ same using t_FFELT
  %7 =
  [x^2 + 1 1]
  
  ? factormodDDF(x^2-Mod(1,3))
  %8 =
  [Mod(1, 3)*x^2 + Mod(2, 3) 1]
  

The library syntax is GEN factormodDDF(GEN f, GEN D = NULL).


factormodSQF(f, {D})

Squarefree factorization of the polynomial f over the finite field defined by the domain D as follows:

* D = p a prime: factor over 𝔽p;

* D = [T,p] for a prime p and T an irreducible polynomial over 𝔽p: factor over 𝔽p[x]/(T);

* D a t_FFELT: factor over the attached field;

* D omitted: factor over the field of definition of f, which must be a finite field.

This is somewhat faster than full factorization. The coefficients of f must be operation-compatible with the corresponding finite field. The result is a two-column matrix:

* the first column contains monic squarefree pairwise coprime polynomials dividing f;

* the second column contains the power to which the polynomial in column 1 divides f;

The factors are ordered by increasing degree and the result is canonical: it will not change across multiple calls or sessions.

  ? f = (x^2 + 1)^3 * (x^2-1)^2;
  ? factormodSQF(f, 3)  \\ over F3
  %1 =
  [Mod(1, 3)*x^2 + Mod(2, 3) 2]
  
  [Mod(1, 3)*x^2 + Mod(1, 3) 3]
  
  ? for(i=1,10^5,factormodSQF(f,3))
  time = 192 ms.
  ? for(i=1,10^5,factormod(f,3))  \\ full factorization is slower
  time = 409 ms.
  
  ? liftall( factormodSQF((x^2 + 1)^3, [3, t^2+1]) ) \\ over F9
  %4 =
  [x^2 + 1 3]
  
  ? t = ffgen(t^2+Mod(1,3)); factormodSQF((x^2 + t^0)^3) \\ same using t_FFELT
  %5 =
  [x^2 + 1 3]
  
  ? factormodSQF(x^8 + x^7 + x^6 + x^2 + x + Mod(1,2))
  %6 =
  [                Mod(1, 2)*x + Mod(1, 2) 2]
  
  [Mod(1, 2)*x^2 + Mod(1, 2)*x + Mod(1, 2) 3]

The library syntax is GEN factormodSQF(GEN f, GEN D = NULL).


ffcompomap(f, g)

Let k, l, m be three finite fields and f a (partial) map from k to l and g a (partial) map from l to m, return the (partial) map f o g from k to m.

  a = ffgen([3,5],'a); b = ffgen([3,10],'b); c = ffgen([3,20],'c);
  m = ffembed(a, b); n = ffembed(b, c);
  rm = ffinvmap(m); rn = ffinvmap(n);
  nm = ffcompomap(n,m);
  ffmap(n,ffmap(m,a)) == ffmap(nm, a)
  %5 = 1
  ffcompomap(rm, rn) == ffinvmap(nm)
  %6 = 1

The library syntax is GEN ffcompomap(GEN f, GEN g).


ffembed(a, b)

Given two finite fields elements a and b, return a map embedding the definition field of a to the definition field of b. Assume that the latter contains the former.

  ? a = ffgen([3,5],'a);
  ? b = ffgen([3,10],'b);
  ? m = ffembed(a, b);
  ? A = ffmap(m, a);
  ? minpoly(A) == minpoly(a)
  %5 = 1

The library syntax is GEN ffembed(GEN a, GEN b).


ffextend(a, P, {v})

Extend the field K of definition of a by a root of the polynomial P ∈ K[X] assumed to be irreducible over K. Return [r, m] where r is a root of P in the extension field L and m is a map from K to L, see ffmap. If v is given, the variable name is used to display the generator of L, else the name of the variable of P is used. A generator of L can be recovered using b = ffgen(r). The image of P in L[X] can be recovered using PL = ffmap(m,P).

  ? a = ffgen([3,5],'a);
  ? P = x^2-a; polisirreducible(P)
  %2 = 1
  ? [r,m] = ffextend(a, P, 'b);
  ? r
  %3 = b^9+2*b^8+b^7+2*b^6+b^4+1
  ? subst(ffmap(m, P), x, r)
  %4 = 0
  ? ffgen(r)
  %5 = b

The library syntax is GEN ffextend(GEN a, GEN P, long v = -1) where v is a variable number.


fffrobenius(m, {n = 1})

Return the n-th power of the Frobenius map over the field of definition of m.

  ? a = ffgen([3,5],'a);
  ? f = fffrobenius(a);
  ? ffmap(f,a) == a^3
  %3 = 1
  ? g = fffrobenius(a, 5);
  ? ffmap(g,a) == a
  %5 = 1
  ? h = fffrobenius(a, 2);
  ? h == ffcompomap(f,f)
  %7 = 1

The library syntax is GEN fffrobenius(GEN m, long n).


ffgen(k, {v = 'x})

Return a generator for the finite field k as a t_FFELT. The field k can be given by

* its order q

* the pair [p,f] where q = p^f

* a monic irreducible polynomial with t_INTMOD coefficients modulo a prime.

* a t_FFELT belonging to k.

If v is given, the variable name is used to display g, else the variable of the polynomial or the t_FFELT is used , else x is used.

When only the order is specified, the function uses the polynomial generated by ffinit and is deterministic: two calls to the function with the same parameters will always give the same generator.

For efficiency, the characteristic is not checked to be prime, nor If a polynomial is given it is not checked whether it is irreducible.

To obtain a multiplicative generator, call ffprimroot on the result.

  ? g = ffgen(16, 't);
  ? g.mod \\ recover the underlying polynomial.
  %2 = t^4+t^3+t^2+t+1
  ? g.p \\ recover the characteristic
  %3 = 2
  ? fforder(g) \\ g is not a multiplicative generator
  %4 = 5
  ? a = ffprimroot(g) \\ recover a multiplicative generator
  %5 = t^3+t^2+t
  ? fforder(a)
  %6 = 15

The library syntax is GEN ffgen(GEN k, long v = -1) where v is a variable number.

To create a generator for a prime finite field, the function GEN p_to_GEN(GEN p, long v) returns ffgen(p,v)^0.


ffinit(p, n, {v = 'x})

Computes a monic polynomial of degree n which is irreducible over 𝔽p, where p is assumed to be prime. This function uses a fast variant of Adleman and Lenstra's algorithm.

It is useful in conjunction with ffgen; for instance if P = ffinit(3,2), you can represent elements in 𝔽3^2 in term of g = ffgen(P,'t). This can be abbreviated as g = ffgen(3^2, 't), where the defining polynomial P can be later recovered as g.mod.

The library syntax is GEN ffinit(GEN p, long n, long v = -1) where v is a variable number.


ffinvmap(m)

m being a map from K to L two finite fields, return the partial map p from L to K such that for all k ∈ K, p(m(k)) = k.

  ? a = ffgen([3,5],'a);
  ? b = ffgen([3,10],'b);
  ? m = ffembed(a, b);
  ? p = ffinvmap(m);
  ? u = random(a);
  ? v = ffmap(m, u);
  ? ffmap(p, v^2+v+2) == u^2+u+2
  %7 = 1

The library syntax is GEN ffinvmap(GEN m).


fflog(x, g, {o})

Discrete logarithm of the finite field element x in base g, i.e.  an e in ℤ such that g^e = o. If present, o represents the multiplicative order of g, see Section se:DLfun; the preferred format for this parameter is [ord, factor(ord)], where ord is the order of g. It may be set as a side effect of calling ffprimroot. The result is undefined if e does not exist. This function uses

* a combination of generic discrete log algorithms (see znlog)

* a cubic sieve index calculus algorithm for large fields of degree at least 5.

* Coppersmith's algorithm for fields of characteristic at most 5.

  ? t = ffgen(ffinit(7,5));
  ? o = fforder(t)
  %2 = 5602   \\  not a primitive root.
  ? fflog(t^10,t)
  %3 = 10
  ? fflog(t^10,t, o)
  %4 = 10
  ? g = ffprimroot(t, &o);
  ? o   \\ order is 16806, bundled with its factorization matrix
  %6 = [16806, [2, 1; 3, 1; 2801, 1]]
  ? fforder(g, o)
  %7 = 16806
  ? fflog(g^10000, g, o)
  %8 = 10000

The library syntax is GEN fflog(GEN x, GEN g, GEN o = NULL).


ffmap(m, x)

Given a (partial) map m between two finite fields, return the image of x by m. The function is applied recursively to the component of vectors, matrices and polynomials. If m is a partial map that is not defined at x, return [].

  ? a = ffgen([3,5],'a);
  ? b = ffgen([3,10],'b);
  ? m = ffembed(a, b);
  ? P = x^2+a*x+1;
  ? Q = ffmap(m,P);
  ? ffmap(m,poldisc(P)) == poldisc(Q)
  %6 = 1

The library syntax is GEN ffmap(GEN m, GEN x).


ffnbirred(q, n, {fl = 0})

Computes the number of monic irreducible polynomials over 𝔽q of degree exactly n, (flag = 0 or omitted) or at most n (flag = 1).

The library syntax is GEN ffnbirred0(GEN q, long n, long fl). Also available are GEN ffnbirred(GEN q, long n) (for flag = 0) and GEN ffsumnbirred(GEN q, long n) (for flag = 1).


fforder(x, {o})

Multiplicative order of the finite field element x. If o is present, it represents a multiple of the order of the element, see Section se:DLfun; the preferred format for this parameter is [N, factor(N)], where N is the cardinality of the multiplicative group of the underlying finite field.

  ? t = ffgen(ffinit(nextprime(10^8), 5));
  ? g = ffprimroot(t, &o);  \\  o will be useful!
  ? fforder(g^1000000, o)
  time = 0 ms.
  %5 = 5000001750000245000017150000600250008403
  ? fforder(g^1000000)
  time = 16 ms. \\  noticeably slower, same result of course
  %6 = 5000001750000245000017150000600250008403

The library syntax is GEN fforder(GEN x, GEN o = NULL).


ffprimroot(x, {&o})

Return a primitive root of the multiplicative group of the definition field of the finite field element x (not necessarily the same as the field generated by x). If present, o is set to a vector [ord, fa], where ord is the order of the group and fa its factorization factor(ord). This last parameter is useful in fflog and fforder, see Section se:DLfun.

  ? t = ffgen(ffinit(nextprime(10^7), 5));
  ? g = ffprimroot(t, &o);
  ? o[1]
  %3 = 100000950003610006859006516052476098
  ? o[2]
  %4 =
  [2 1]
  
  [7 2]
  
  [31 1]
  
  [41 1]
  
  [67 1]
  
  [1523 1]
  
  [10498781 1]
  
  [15992881 1]
  
  [46858913131 1]
  
  ? fflog(g^1000000, g, o)
  time = 1,312 ms.
  %5 = 1000000

The library syntax is GEN ffprimroot(GEN x, GEN *o = NULL).


gcd(x, {y})

Creates the greatest common divisor of x and y. If you also need the u and v such that x*u + y*v = gcd(x,y), use the gcdext function. x and y can have rather quite general types, for instance both rational numbers. If y is omitted and x is a vector, returns the gcd of all components of x, i.e. this is equivalent to content(x).

When x and y are both given and one of them is a vector/matrix type, the GCD is again taken recursively on each component, but in a different way. If y is a vector, resp. matrix, then the result has the same type as y, and components equal to gcd(x, y[i]), resp. gcd(x, y[,i]). Else if x is a vector/matrix the result has the same type as x and an analogous definition. Note that for these types, gcd is not commutative.

The algorithm used is a naive Euclid except for the following inputs:

* integers: use modified right-shift binary ("plus-minus" variant).

* univariate polynomials with coefficients in the same number field (in particular rational): use modular gcd algorithm.

* general polynomials: use the subresultant algorithm if coefficient explosion is likely (non modular coefficients).

If u and v are polynomials in the same variable with inexact coefficients, their gcd is defined to be scalar, so that

  ? a = x + 0.0; gcd(a,a)
  %1 = 1
  ? b = y*x + O(y); gcd(b,b)
  %2 = y
  ? c = 4*x + O(2^3); gcd(c,c)
  %3 = 4

A good quantitative check to decide whether such a gcd "should be" non-trivial, is to use polresultant: a value close to 0 means that a small deformation of the inputs has non-trivial gcd. You may also use gcdext, which does try to compute an approximate gcd d and provides u, v to check whether u x + v y is close to d.

The library syntax is GEN ggcd0(GEN x, GEN y = NULL). Also available are GEN ggcd(GEN x, GEN y), if y is not NULL, and GEN content(GEN x), if y = NULL.


gcdext(x, y)

Returns [u,v,d] such that d is the gcd of x,y, x*u+y*v = gcd(x,y), and u and v minimal in a natural sense. The arguments must be integers or polynomials.

  ? [u, v, d] = gcdext(32,102)
  %1 = [16, -5, 2]
  ? d
  %2 = 2
  ? gcdext(x^2-x, x^2+x-2)
  %3 = [-1/2, 1/2, x - 1]

If x,y are polynomials in the same variable and inexact coefficients, then compute u,v,d such that x*u+y*v = d, where d approximately divides both and x and y; in particular, we do not obtain gcd(x,y) which is defined to be a scalar in this case:

  ? a = x + 0.0; gcd(a,a)
  %1 = 1
  
  ? gcdext(a,a)
  %2 = [0, 1, x + 0.E-28]
  
  ? gcdext(x-Pi, 6*x^2-zeta(2))
  %3 = [-6*x - 18.8495559, 1, 57.5726923]

For inexact inputs, the output is thus not well defined mathematically, but you obtain explicit polynomials to check whether the approximation is close enough for your needs.

The library syntax is GEN gcdext0(GEN x, GEN y).


hilbert(x, y, {p})

Hilbert symbol of x and y modulo the prime p, p = 0 meaning the place at infinity (the result is undefined if p != 0 is not prime).

It is possible to omit p, in which case we take p = 0 if both x and y are rational, or one of them is a real number. And take p = q if one of x, y is a t_INTMOD modulo q or a q-adic. (Incompatible types will raise an error.)

The library syntax is long hilbert(GEN x, GEN y, GEN p = NULL).


isfundamental(D)

True (1) if D is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise. D can be input in factored form as for arithmetic functions:

  ? isfundamental(factor(-8))
  %1 = 1
  \\ count fundamental discriminants up to 10^8
  ? c = 0; forfactored(d = 1, 10^8, if (isfundamental(d), c++)); c
  time = 40,840 ms.
  %2 = 30396325
  ? c = 0; for(d = 1, 10^8, if (isfundamental(d), c++)); c
  time = 1min, 33,593 ms. \\ slower !
  %3 = 30396325

The library syntax is long isfundamental(GEN D).


ispolygonal(x, s, {&N})

True (1) if the integer x is an s-gonal number, false (0) if not. The parameter s > 2 must be a t_INT. If N is given, set it to n if x is the n-th s-gonal number.

  ? ispolygonal(36, 3, &N)
  %1 = 1
  ? N

The library syntax is long ispolygonal(GEN x, GEN s, GEN *N = NULL).


ispower(x, {k}, {&n})

If k is given, returns true (1) if x is a k-th power, false (0) if not. What it means to be a k-th power depends on the type of x; see issquare for details.

If k is omitted, only integers and fractions are allowed for x and the function returns the maximal k ≥ 2 such that x = n^k is a perfect power, or 0 if no such k exist; in particular ispower(-1), ispower(0), and ispower(1) all return 0.

If a third argument &n is given and x is indeed a k-th power, sets n to a k-th root of x.

For a t_FFELT x, instead of omitting k (which is not allowed for this type), it may be natural to set

  k = (x.p ^ x.f - 1) / fforder(x)

The library syntax is long ispower(GEN x, GEN k = NULL, GEN *n = NULL). Also available is long gisanypower(GEN x, GEN *pty) (k omitted).


ispowerful(x)

True (1) if x is a powerful integer, false (0) if not; an integer is powerful if and only if its valuation at all primes dividing x is greater than 1.

  ? ispowerful(50)
  %1 = 0
  ? ispowerful(100)
  %2 = 1
  ? ispowerful(5^3*(10^1000+1)^2)
  %3 = 1

The library syntax is long ispowerful(GEN x).


isprime(x, {flag = 0})

True (1) if x is a prime number, false (0) otherwise. A prime number is a positive integer having exactly two distinct divisors among the natural numbers, namely 1 and itself.

This routine proves or disproves rigorously that a number is prime, which can be very slow when x is indeed a large prime integer. For instance a 1000 digits prime should require 15 to 30 minutes with default algorithms. Use ispseudoprime to quickly check for compositeness. Use primecert in order to obtain a primality proof instead of a yes/no answer; see also factor.

The function accepts vector/matrices arguments, and is then applied componentwise.

If flag = 0, use a combination of

* Baillie-Pomerance-Selfridge-Wagstaff compositeness test (see ispseudoprime),

* Selfridge "p-1" test if x-1 is smooth enough,

* Adleman-Pomerance-Rumely-Cohen-Lenstra (APRCL) for general medium-sized x (less than 1500 bits),

* Atkin-Morain's Elliptic Curve Primality Prover (ECPP) for general large x.

If flag = 1, use Selfridge-Pocklington-Lehmer "p-1" test; this requires partially factoring various auxilliary integers and is likely to be very slow.

If flag = 2, use APRCL only.

If flag = 3, use ECPP only.

The library syntax is GEN gisprime(GEN x, long flag).


isprimepower(x, {&n})

If x = p^k is a prime power (p prime, k > 0), return k, else return 0. If a second argument &n is given and x is indeed the k-th power of a prime p, sets n to p.

The library syntax is long isprimepower(GEN x, GEN *n = NULL).


ispseudoprime(x, {flag})

True (1) if x is a strong pseudo prime (see below), false (0) otherwise. If this function returns false, x is not prime; if, on the other hand it returns true, it is only highly likely that x is a prime number. Use isprime (which is of course much slower) to prove that x is indeed prime. The function accepts vector/matrices arguments, and is then applied componentwise.

If flag = 0, checks whether x has no small prime divisors (up to 101 included) and is a Baillie-Pomerance-Selfridge-Wagstaff pseudo prime. Such a pseudo prime passes a Rabin-Miller test for base 2, followed by a Lucas test for the sequence (P,1), where P ≥ 3 is the smallest odd integer such that P^2 - 4 is not a square mod x. (Technically, we are using an "almost extra strong Lucas test" that checks whether Vn is ± 2, without computing Un.)

There are no known composite numbers passing the above test, although it is expected that infinitely many such numbers exist. In particular, all composites ≤ 264 are correctly detected (checked using http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html).

If flag > 0, checks whether x is a strong Miller-Rabin pseudo prime for flag randomly chosen bases (with end-matching to catch square roots of -1).

The library syntax is GEN gispseudoprime(GEN x, long flag).


ispseudoprimepower(x, {&n})

If x = p^k is a pseudo-prime power (p pseudo-prime as per ispseudoprime, k > 0), return k, else return 0. If a second argument &n is given and x is indeed the k-th power of a prime p, sets n to p.

More precisely, k is always the largest integer such that x = n^k for some integer n and, when n ≤ 264 the function returns k > 0 if and only if n is indeed prime. When n > 264 is larger than the threshold, the function may return 1 even though n is composite: it only passed an ispseudoprime(n) test.

The library syntax is long ispseudoprimepower(GEN x, GEN *n = NULL).


issquare(x, {&n})

True (1) if x is a square, false (0) if not. What "being a square" means depends on the type of x: all t_COMPLEX are squares, as well as all non-negative t_REAL; for exact types such as t_INT, t_FRAC and t_INTMOD, squares are numbers of the form s^2 with s in ℤ, ℚ and ℤ/Nℤ respectively.

  ? issquare(3)          \\ as an integer
  %1 = 0
  ? issquare(3.)         \\ as a real number
  %2 = 1
  ? issquare(Mod(7, 8))  \\ in Z/8Z
  %3 = 0
  ? issquare( 5 + O(13^4) )  \\ in Q_13
  %4 = 0

If n is given, a square root of x is put into n.

  ? issquare(4, &n)
  %1 = 1
  ? n
  %2 = 2

For polynomials, either we detect that the characteristic is 2 (and check directly odd and even-power monomials) or we assume that 2 is invertible and check whether squaring the truncated power series for the square root yields the original input.

For t_POLMOD x, we only support t_POLMODs of t_INTMODs encoding finite fields, assuming without checking that the intmod modulus p is prime and that the polmod modulus is irreducible modulo p.

  ? issquare(Mod(Mod(2,3), x^2+1), &n)
  %1 = 1
  ? n
  %2 = Mod(Mod(2, 3)*x, Mod(1, 3)*x^2 + Mod(1, 3))

The library syntax is long issquareall(GEN x, GEN *n = NULL). Also available is long issquare(GEN x). Deprecated GP-specific functions GEN gissquare(GEN x) and GEN gissquareall(GEN x, GEN *pt) return gen0 and gen1 instead of a boolean value.


issquarefree(x)

True (1) if x is squarefree, false (0) if not. Here x can be an integer or a polynomial with coefficients in an integral domain.

  ? issquarefree(12)
  %1 = 0
  ? issquarefree(6)
  %2 = 1
  ? issquarefree(x^3+x^2)
  %3 = 0
  ? issquarefree(Mod(1,4)*(x^2+x+1))    \\ Z/4Z is not a domain !
   ***   at top-level: issquarefree(Mod(1,4)*(x^2+x+1))
   ***                 ^ —  —  —  —  —  —  —  —  —  — --
   *** issquarefree: impossible inverse in Fp_inv: Mod(2, 4).

A polynomial is declared squarefree if gcd(x,x') is 1. In particular a non-zero polynomial with inexact coefficients is considered to be squarefree. Note that this may be inconsistent with factor, which first rounds the input to some exact approximation before factoring in the apropriate domain; this is correct when the input is not close to an inseparable polynomial (the resultant of x and x' is not close to 0).

An integer can be input in factored form as in arithmetic functions.

  ? issquarefree(factor(6))
  %1 = 1
  \\ count squarefree integers up to 10^8
  ? c = 0; for(d = 1, 10^8, if (issquarefree(d), c++)); c
  time = 3min, 2,590 ms.
  %2 = 60792694
  ? c = 0; forfactored(d = 1, 10^8, if (issquarefree(d), c++)); c
  time = 45,348 ms. \\ faster !
  %3 = 60792694

The library syntax is long issquarefree(GEN x).


istotient(x, {&N})

True (1) if x = φ(n) for some integer n, false (0) if not.

  ? istotient(14)
  %1 = 0
  ? istotient(100)
  %2 = 0

If N is given, set N = n as well.

  ? istotient(4, &n)
  %1 = 1
  ? n
  %2 = 10

The library syntax is long istotient(GEN x, GEN *N = NULL).


kronecker(x, y)

Kronecker symbol (x|y), where x and y must be of type integer. By definition, this is the extension of Legendre symbol to ℤ x ℤ by total multiplicativity in both arguments with the following special rules for y = 0, -1 or 2:

* (x|0) = 1 if |x |= 1 and 0 otherwise.

* (x|-1) = 1 if x ≥ 0 and -1 otherwise.

* (x|2) = 0 if x is even and 1 if x = 1,-1 mod 8 and -1 if x = 3,-3 mod 8.

The library syntax is long kronecker(GEN x, GEN y).


lcm(x, {y})

Least common multiple of x and y, i.e. such that lcm(x,y)*gcd(x,y) = x*y, up to units. If y is omitted and x is a vector, returns the lcm of all components of x. For integer arguments, return the non-negative lcm.

When x and y are both given and one of them is a vector/matrix type, the LCM is again taken recursively on each component, but in a different way. If y is a vector, resp. matrix, then the result has the same type as y, and components equal to lcm(x, y[i]), resp. lcm(x, y[,i]). Else if x is a vector/matrix the result has the same type as x and an analogous definition. Note that for these types, lcm is not commutative.

Note that lcm(v) is quite different from

  l = v[1]; for (i = 1, #v, l = lcm(l, v[i]))

Indeed, lcm(v) is a scalar, but l may not be (if one of the v[i] is a vector/matrix). The computation uses a divide-conquer tree and should be much more efficient, especially when using the GMP multiprecision kernel (and more subquadratic algorithms become available):

  ? v = vector(10^5, i, random);
  ? lcm(v);
  time = 546 ms.
  ? l = v[1]; for (i = 1, #v, l = lcm(l, v[i]))
  time = 4,561 ms.

The library syntax is GEN glcm0(GEN x, GEN y = NULL).


logint(x, b, {&z})

Return the largest integer e so that b^e ≤ x, where the parameters b > 1 and x > 0 are both integers. If the parameter z is present, set it to b^e.

  ? logint(1000, 2)
  %1 = 9
  ? 2^9
  %2 = 512
  ? logint(1000, 2, &z)
  %3 = 9
  ? z
  %4 = 512

The number of digits used to write b in base x is 1 + logint(x,b):

  ? #digits(1000!, 10)
  %5 = 2568
  ? logint(1000!, 10)
  %6 = 2567

This function may conveniently replace

    floor( log(x) / log(b) )

which may not give the correct answer since PARI does not guarantee exact rounding.

The library syntax is long logint0(GEN x, GEN b, GEN *z = NULL).


moebius(x)

Moebius μ-function of |x|; x must be a non-zero integer.

The library syntax is long moebius(GEN x).


nextprime(x)

Finds the smallest pseudoprime (see ispseudoprime) greater than or equal to x. x can be of any real type. Note that if x is a pseudoprime, this function returns x and not the smallest pseudoprime strictly larger than x. To rigorously prove that the result is prime, use isprime.

The library syntax is GEN nextprime(GEN x).


numdiv(x)

Number of divisors of |x|. x must be of type integer.

The library syntax is GEN numdiv(GEN x).


omega(x)

Number of distinct prime divisors of |x|. x must be of type integer.

  ? factor(392)
  %1 =
  [2 3]
  
  [7 2]
  
  ? omega(392)
  %2 = 2;  \\ without multiplicity
  ? bigomega(392)
  %3 = 5;  \\ = 3+2, with multiplicity

The library syntax is long omega(GEN x).


precprime(x)

Finds the largest pseudoprime (see ispseudoprime) less than or equal to x. x can be of any real type. Returns 0 if x ≤ 1. Note that if x is a prime, this function returns x and not the largest prime strictly smaller than x. To rigorously prove that the result is prime, use isprime.

The library syntax is GEN precprime(GEN x).


prime(n)

The n-th prime number

  ? prime(10^9)
  %1 = 22801763489

Uses checkpointing and a naive O(n) algorithm. Will need about 30 minutes for n up to 1011; make sure to start gp with primelimit at least sqrt{pn}, e.g. the value sqrt{nlog (nlog n)} is guaranteed to be sufficient.

The library syntax is GEN prime(long n).


primecert(N, {flag = 0})

If N is a prime, return a PARI Primality Certificate for the prime N, as described below. Otherwise, return 0. A Primality Certificate c can be checked using primecertisvalid(c).

If flag = 0 (default), return an ECPP certificate (Atkin-Morain)

A PARI ECPP Primality Certificate for the prime N is either a prime integer N < 264 or a vector C of length ℓ whose ith component C[i] is a vector of length 5 containing Ni, ti, si, ai, Pi (in that order) where N1 = N. It is said to be valid if for each i = 1,..., ℓ, all of the following conditions are satisfied

* Ni is a positive integer

* ti is an integer such that ti < 4Ni^2

* si is a positive integer which divides mi where mi = Ni + 1 - ti

* If we set qi = (mi)/(si), then

  * qi > (Ni1/4+1)^2

  * qi = Ni+1 if 1 ≤ i < l

  * q_ℓ ≤ 264 is prime

* ai is an integer

  * P[i] is a vector of length 2 representing the affine point Pi = (xi, yi) on the elliptic curve E: y^2 = x^3 + a_ix + bi modulo Ni where bi = yi^2 - xi^3 - a_ixi satisfying the following:

  * mi Pi = oo

  * si Pi != oo

Theorem. If N is an integer and there exist positive integers m, q and a point P on the elliptic curve E: y^2 = x^3 + ax + b defined modulo N such that q > (N1/4 + 1)^2, q is a prime divisor of m, mP = oo and (m)/(q)P != oo , then N is prime.

Using the above theorem, the data in the vector C can be used to certify the primality of N assuming q_ℓ is prime.

  ? primecert(10^35 + 69)
  %1 = [[100000000000000000000000000000000069, 5468679110354
  52074, 2963504668391148, 0, [60737979324046450274283740674
  208692, 24368673584839493121227731392450025]], [3374383076
  4501150277, -11610830419, 734208843, 0, [26740412374402652
  72 4, 6367191119818901665]], [45959444779, 299597, 2331, 0
  , [18022351516, 9326882 51]]]
  ? primecert(nextprime(2^64))
  %2 = [[18446744073709551629, -8423788454, 160388, 1, [1059
  8342506117936052, 2225259013356795550]]]
  ? primecert(6)
  %3 = 0
  ? primecert(41)
  %4 = 41

If flag = 1 (very slow), return an N-1 certificate (Pocklington Lehmer)

A PARI N-1 Primality Certificate for the prime N is either a prime integer N < 264 or a pair [N, C], where C is a vector with ℓ elements which are either a single integer pi < 264 or a triple [pi,ai,Ci] with pi > 264 satisfying the following properties:

* pi is a prime divisor of N - 1;

* ai is an integer such that aiN-1 = 1 (mod N) and ai(N-1)/pi - 1 is coprime with N;

* Ci is an N-1 Primality Certificate for pi

* The product F of the pivpi(N-1) is strictly larger than N1/3. Provided that all pi are indeed primes, this implies that any divisor of N is congruent to 1 modulo F.

* The Billhart, Lehmer, Selfridge criterion is satisfied: when we write N = 1 + c1 F + c2 F^2 in base F the polynomial 1 + c1 X + c2 X^2 is irreducible over ℤ, i.e. c1^2 - 4c2 is not a square. This implies that N is prime.

This algorithm requires factoring partially p-1 for various prime integers p with an unfactored parted ≤ p2/3 and this may be exceedingly slow compared to the default.

The algorithm fails if one of the pseudo-prime factors is not prime, which is exceedingly unlikely and well worth a bug report. Note that if you monitor the algorithm at a high enough debug level, you may see warnings about untested integers being declared primes. This is normal: we ask for partial factorizations (sufficient to prove primality if the unfactored part is not too large), and factor warns us that the cofactor hasn't been tested. It may or may not be tested later, and may or may not be prime. This does not affect the validity of the whole Primality Certificate.

The library syntax is GEN primecert(GEN N, long flag).


primecertexport(cert, {format = 0})

Returns a string suitable for print/write to display a primality certificate from primecert, the format of which depends on the value of format:

* 0 (default): Human-readable format.

* 1: Primo format 4.

* 2: MAGMA format.

Currently, only ECPP Primality Certificates are supported.

  ? cert = primecert(10^35+69);
  ? s = primecertexport(cert); \\ Human-readable
  ? print(s)
  [1]
   N = 100000000000000000000000000000000069
   t = 546867911035452074
   s = 2963504668391148
  a4 = 0
  D = -3
  m = 99999999999999999453132088964547996
  q = 33743830764501150277
  E = [0, 1]
  P = [21567861682493263464353543707814204,
  49167839501923147849639425291163552]
  [2]
   N = 33743830764501150277
   t = -11610830419
   s = 734208843
  a4 = 0
  D = -3
  m = 33743830776111980697
  q = 45959444779
  E = [0, 25895956964997806805]
  P = [29257172487394218479, 3678591960085668324]
  
  \\ Primo format
  ? s = primecertexport(cert,1); write("cert.out", s);
  
  \\ Magma format, write to file
  ? s = primecertexport(cert,2); write("cert.m", s);
  
  ? cert = primecert(10^35+69, 1); \\ N-1 certificate
   ***   at top-level: primecertexport(cert)
   ***                 ^ —  —  —  —  —  —  — 
   *** primecertexport: sorry, N-1 certificate is not yet implemented.

The library syntax is GEN primecertexport(GEN cert, long format).


primecertisvalid(cert)

Verifies if cert is a valid PARI ECPP Primality certificate, as described in ??primecert.

  ? cert = primecert(10^35 + 69)
  %1 = [[100000000000000000000000000000000069, 5468679110354
  52074, 2963504668391148, 0, [60737979324046450274283740674
  208692, 24368673584839493121227731392450025]], [3374383076
  4501150277, -11610830419, 734208843, 0, [26740412374402652
  72 4, 6367191119818901665]], [45959444779, 299597, 2331, 0
  , [18022351516, 9326882 51]]]
  ? primecertisvalid(cert)
  %2 = 1
  
  ? cert[1][1]++; \\ random perturbation
  ? primecertisvalid(cert)
  %4 = 0  \\ no longer valid
  ? primecertisvalid(primecert(6))
  %5 = 0

The library syntax is long primecertisvalid(GEN cert).


primepi(x)

The prime counting function. Returns the number of primes p, p ≤ x.

  ? primepi(10)
  %1 = 4;
  ? primes(5)
  %2 = [2, 3, 5, 7, 11]
  ? primepi(10^11)
  %3 = 4118054813

Uses checkpointing and a naive O(x) algorithm; make sure to start gp with primelimit at least sqrt{x}.

The library syntax is GEN primepi(GEN x).


primes(n)

Creates a row vector whose components are the first n prime numbers. (Returns the empty vector for n ≤ 0.) A t_VEC n = [a,b] is also allowed, in which case the primes in [a,b] are returned

  ? primes(10)     \\ the first 10 primes
  %1 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
  ? primes([0,29])  \\ the primes up to 29
  %2 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
  ? primes([15,30])
  %3 = [17, 19, 23, 29]

The library syntax is GEN primes0(GEN n).


qfbclassno(D, {flag = 0})

Ordinary class number of the quadratic order of discriminant D, for "small" values of D.

* if D > 0 or flag = 1, use a O(|D|1/2) algorithm (compute L(1,χD) with the approximate functional equation). This is slower than quadclassunit as soon as |D| ~ 10^2 or so and is not meant to be used for large D.

* if D < 0 and flag = 0 (or omitted), use a O(|D|1/4) algorithm (Shanks's baby-step/giant-step method). It should be faster than quadclassunit for small values of D, say |D| < 1018.

Important warning. In the latter case, this function only implements part of Shanks's method (which allows to speed it up considerably). It gives unconditionnally correct results for |D| < 2.1010, but may give incorrect results for larger values if the class group has many cyclic factors. We thus recommend to double-check results using the function quadclassunit, which is about 2 to 3 times slower in the range |D| ∈ [1010, 1018], assuming GRH. We currently have no counter-examples but they should exist: we would appreciate a bug report if you find one.

Warning. Contrary to what its name implies, this routine does not compute the number of classes of binary primitive forms of discriminant D, which is equal to the narrow class number. The two notions are the same when D < 0 or the fundamental unit ϵ has negative norm; when D > 0 and Nϵ > 0, the number of classes of forms is twice the ordinary class number. This is a problem which we cannot fix for backward compatibility reasons. Use the following routine if you are only interested in the number of classes of forms:

  QFBclassno(D) =
  qfbclassno(D) * if (D < 0 || norm(quadunit(D)) < 0, 1, 2)

Here are a few examples:

  ? qfbclassno(400000028) \\ D > 0: slow
  time = 3,140 ms.
  %1 = 1
  ? quadclassunit(400000028).no
  time = 20 ms. \\ { much faster, assume GRH}
  %2 = 1
  ? qfbclassno(-400000028) \\ D < 0: fast enough
  time = 0 ms.
  %3 = 7253
  ? quadclassunit(-400000028).no
  time = 0 ms.
  %4 = 7253

See also qfbhclassno.

The library syntax is GEN qfbclassno0(GEN D, long flag).


qfbcompraw(x, y)

composition of the binary quadratic forms x and y, without reduction of the result. This is useful e.g. to compute a generating element of an ideal. The result is undefined if x and y do not have the same discriminant.

The library syntax is GEN qfbcompraw(GEN x, GEN y).


qfbhclassno(x)

Hurwitz class number of x, when x is non-negative and congruent to 0 or 3 modulo 4, and 0 for other values. For x > 5.10^5, we assume the GRH, and use quadclassunit with default parameters.

  ? qfbhclassno(1) \\ not 0 or 3 mod 4
  %1 = 0
  ? qfbhclassno(3)
  %2 = 1/3
  ? qfbhclassno(4)
  %3 = 1/2
  ? qfbhclassno(23)
  %4 = 3

The library syntax is GEN hclassno(GEN x).


qfbnucomp(x, y, L)

composition of the primitive positive definite binary quadratic forms x and y (type t_QFI) using the NUCOMP and NUDUPL algorithms of Shanks, à la Atkin. L is any positive constant, but for optimal speed, one should take L = |D/4|1/4, i.e. sqrtnint(abs(D) >> 2,4), where D is the common discriminant of x and y. When x and y do not have the same discriminant, the result is undefined.

The current implementation is slower than the generic routine for small D, and becomes faster when D has about 45 bits.

The library syntax is GEN nucomp(GEN x, GEN y, GEN L). Also available is GEN nudupl(GEN x, GEN L) when x = y.


qfbnupow(x, n, {L})

n-th power of the primitive positive definite binary quadratic form x using Shanks's NUCOMP and NUDUPL algorithms; if set, L should be equal to sqrtnint(abs(D) >> 2,4), where D < 0 is the discriminant of x.

The current implementation is slower than the generic routine for small discriminant D, and becomes faster for D ~ 245.

The library syntax is GEN nupow(GEN x, GEN n, GEN L = NULL).


qfbpowraw(x, n)

n-th power of the binary quadratic form x, computed without doing any reduction (i.e. using qfbcompraw). Here n must be non-negative and n < 231.

The library syntax is GEN qfbpowraw(GEN x, long n).


qfbprimeform(x, p)

Prime binary quadratic form of discriminant x whose first coefficient is p, where |p| is a prime number. By abuse of notation, p = ± 1 is also valid and returns the unit form. Returns an error if x is not a quadratic residue mod p, or if x < 0 and p < 0. (Negative definite t_QFI are not implemented.) In the case where x > 0, the "distance" component of the form is set equal to zero according to the current precision.

The library syntax is GEN primeform(GEN x, GEN p, long prec).


qfbred(x, {flag = 0}, {d}, {isd}, {sd})

Reduces the binary quadratic form x (updating Shanks's distance function if x is indefinite). The binary digits of flag are toggles meaning

  1: perform a single reduction step

  2: don't update Shanks's distance

The arguments d, isd, sd, if present, supply the values of the discriminant, floor{sqrt{d}}, and sqrt{d} respectively (no checking is done of these facts). If d < 0 these values are useless, and all references to Shanks's distance are irrelevant.

The library syntax is GEN qfbred0(GEN x, long flag, GEN d = NULL, GEN isd = NULL, GEN sd = NULL). Also available are

GEN redimag(GEN x) (for definite x),

and for indefinite forms:

GEN redreal(GEN x)

GEN rhoreal(GEN x) ( = qfbred(x,1)),

GEN redrealnod(GEN x, GEN isd) ( = qfbred(x,2,,isd)),

GEN rhorealnod(GEN x, GEN isd) ( = qfbred(x,3,,isd)).


qfbredsl2(x, {data})

Reduction of the (real or imaginary) binary quadratic form x, return [y,g] where y is reduced and g in SL(2,ℤ) is such that g.x = y; data, if present, must be equal to [D, sqrtint(D)], where D > 0 is the discriminant of x. In case x is a t_QFR, the distance component is unaffected.

The library syntax is GEN qfbredsl2(GEN x, GEN data = NULL).


qfbsolve(Q, n)

Return the list of solutions [x,y] of the equation Q(x,y) = n in coprime integers x and y, where Q is a binary quadratic form and n an integer, up to the action of the special orthogonal group S = SO(Q,ℤ), which is isomorphic to the group of units of positive norm of the quadratic order of discriminant D = disc Q.

The result is a vector of two-components vectors [x,y].

If D > 0, S is infinite. If D < -4, S is of order 2, if D = -3, S is of order 6 and if D = -4, S is of order 4.

The integer n can be given by its factorization matrix.

  ? qfbsolve(Qfb(1,0,2),201)
  %1 = [[-13,-4],[-1,10],[-1,-10],[-13,4]]
  ? N=2^128+1;F=factor(N);
  ? qfbsolve(Qfb(1,0,1),[N,F])
  %3 = [[-16382350221535464479,8479443857936402504],
        [18446744073709551616,-1],[-18446744073709551616,-1],
        [16382350221535464479,8479443857936402504]]

Assuming the factorisation of n is given, for each prime factor p of n, the algorithm used runs in probabilistic polynomial time in p (through the computation of square roots of D modulo 4 p); it is polynomial time in D if Q is imaginary, but exponential time if Q is real (through the computation of a full cycle of reduced forms). In the latter case, note that bnfisprincipal provides a solution in heuristic subexponential time in D assuming the GRH.

The library syntax is GEN qfbsolve(GEN Q, GEN n).


quadclassunit(D, {flag = 0}, {tech = []})

Buchmann-McCurley's sub-exponential algorithm for computing the class group of a quadratic order of discriminant D.

This function should be used instead of qfbclassno or quadregula when D < -1025, D > 1010, or when the structure is wanted. It is a special case of bnfinit, which is slower, but more robust.

The result is a vector v whose components should be accessed using member functions:

* v.no: the class number

* v.cyc: a vector giving the structure of the class group as a product of cyclic groups;

* v.gen: a vector giving generators of those cyclic groups (as binary quadratic forms).

* v.reg: the regulator, computed to an accuracy which is the maximum of an internal accuracy determined by the program and the current default (note that once the regulator is known to a small accuracy it is trivial to compute it to very high accuracy, see the tutorial).

The flag is obsolete and should be left alone. In older versions, it supposedly computed the narrow class group when D > 0, but this did not work at all; use the general function bnfnarrow.

Optional parameter tech is a row vector of the form [c1, c2], where c1 ≤ c2 are non-negative real numbers which control the execution time and the stack size, see se:GRHbnf. The parameter is used as a threshold to balance the relation finding phase against the final linear algebra. Increasing the default c1 means that relations are easier to find, but more relations are needed and the linear algebra will be harder. The default value for c1 is 0 and means that it is taken equal to c2. The parameter c2 is mostly obsolete and should not be changed, but we still document it for completeness: we compute a tentative class group by generators and relations using a factorbase of prime ideals ≤ c1 (log |D|)^2, then prove that ideals of norm ≤ c2 (log |D|)^2 do not generate a larger group. By default an optimal c2 is chosen, so that the result is provably correct under the GRH — a famous result of Bach states that c2 = 6 is fine, but it is possible to improve on this algorithmically. You may provide a smaller c2, it will be ignored (we use the provably correct one); you may provide a larger c2 than the default value, which results in longer computing times for equally correct outputs (under GRH).

The library syntax is GEN quadclassunit0(GEN D, long flag, GEN tech = NULL, long prec). If you really need to experiment with the tech parameter, it is usually more convenient to use GEN Buchquad(GEN D, double c1, double c2, long prec). If only the class number is needed, GEN quadclassno(GEN D) will be faster (still assuming the GRH), but will not provide the group structure. For negative D, |D| < 1020, qfbclassno should be faster but may return a wrong result.


quaddisc(x)

Discriminant of the étale algebra ℚ(sqrt{x}), where x ∈ ℚ*. This is the same as coredisc(d) where d is the integer square-free part of x, so x = d f^2 with f ∈ ℚ* and d ∈ ℤ. This returns 0 for x = 0, 1 for x square and the discriminant of the quadratic field ℚ(sqrt{x}) otherwise.

  ? quaddisc(7)
  %1 = 28
  ? quaddisc(-7)
  %2 = -7

The library syntax is GEN quaddisc(GEN x).


quadgen(D, {v = 'w})

Creates the quadratic number ω = (a+sqrt{D})/2 where a = 0 if D = 0 mod 4, a = 1 if D = 1 mod 4, so that (1,ω) is an integral basis for the quadratic order of discriminant D. D must be an integer congruent to 0 or 1 modulo 4, which is not a square. If v is given, the variable name is used to display g else 'w' is used.

  ? w = quadgen(5, 'w); w^2 - w - 1
  %1 = 0
  ? w = quadgen(0, 'w)
   ***   at top-level: w=quadgen(0)
   ***                   ^ —  —  — -
   *** quadgen: domain error in quadpoly: issquare(disc) = 1

The library syntax is GEN quadgen0(GEN D, long v = -1) where v is a variable number.

When v does not matter, the function GEN quadgen(GEN D) is also available.


quadhilbert(D)

Relative equation defining the Hilbert class field of the quadratic field of discriminant D.

If D < 0, uses complex multiplication (Schertz's variant).

If D > 0 Stark units are used and (in rare cases) a vector of extensions may be returned whose compositum is the requested class field. See bnrstark for details.

The library syntax is GEN quadhilbert(GEN D, long prec).


quadpoly(D, {v = 'x})

Creates the "canonical" quadratic polynomial (in the variable v) corresponding to the discriminant D, i.e. the minimal polynomial of quadgen(D). D must be an integer congruent to 0 or 1 modulo 4, which is not a square.

  ? quadpoly(5,'y)
  %1 = y^2 - y - 1
  ? quadpoly(0,'y)
   ***   at top-level: quadpoly(0,'y)
   ***                 ^ —  —  —  — --
   *** quadpoly: domain error in quadpoly: issquare(disc) = 1

The library syntax is GEN quadpoly0(GEN D, long v = -1) where v is a variable number.


quadray(D, f)

Relative equation for the ray class field of conductor f for the quadratic field of discriminant D using analytic methods. A bnf for x^2 - D is also accepted in place of D.

For D < 0, uses the σ function and Schertz's method.

For D > 0, uses Stark's conjecture, and a vector of relative equations may be returned. See bnrstark for more details.

The library syntax is GEN quadray(GEN D, GEN f, long prec).


quadregulator(x)

Regulator of the quadratic field of positive discriminant x. Returns an error if x is not a discriminant (fundamental or not) or if x is a square. See also quadclassunit if x is large.

The library syntax is GEN quadregulator(GEN x, long prec).


quadunit(D, {v = 'w})

Fundamental unit u of the real quadratic field ℚ(sqrt D) where D is the positive discriminant of the field. If D is not a fundamental discriminant, this probably gives the fundamental unit of the corresponding order. D must be an integer congruent to 0 or 1 modulo 4, which is not a square; the result is a quadratic number (see Section se:quadgen). If v is given, the variable name is used to display u else 'w' is used.

The library syntax is GEN quadunit0(GEN D, long v = -1) where v is a variable number.

When v does not matter, the function GEN quadunit(GEN D) is also available.


ramanujantau(n)

Compute the value of Ramanujan's tau function at an individual n, assuming the truth of the GRH (to compute quickly class numbers of imaginary quadratic fields using quadclassunit). Algorithm in Õ(n1/2) using O(log n) space. If all values up to N are required, then ∑ τ(n)q^n = q ∏n ≥ 1 (1-q^n)24 will produce them in time Õ(N), against Õ(N3/2) for individual calls to ramanujantau; of course the space complexity then becomes Õ(N).

  ? tauvec(N) = Vec(q*eta(q + O(q^N))^24);
  ? N = 10^4; v = tauvec(N);
  time = 26 ms.
  ? ramanujantau(N)
  %3 = -482606811957501440000
  ? w = vector(N, n, ramanujantau(n)); \\ much slower !
  time = 13,190 ms.
  ? v == w
  %4 = 1

The library syntax is GEN ramanujantau(GEN n).


randomprime({N = 231})

Returns a strong pseudo prime (see ispseudoprime) in [2,N-1]. A t_VEC N = [a,b] is also allowed, with a ≤ b in which case a pseudo prime a ≤ p ≤ b is returned; if no prime exists in the interval, the function will run into an infinite loop. If the upper bound is less than 264 the pseudo prime returned is a proven prime.

The library syntax is GEN randomprime(GEN N = NULL).


removeprimes({x = []})

Removes the primes listed in x from the prime number table. In particular removeprimes(addprimes()) empties the extra prime table. x can also be a single integer. List the current extra primes if x is omitted.

The library syntax is GEN removeprimes(GEN x = NULL).


sigma(x, {k = 1})

Sum of the k-th powers of the positive divisors of |x|. x and k must be of type integer.

The library syntax is GEN sumdivk(GEN x, long k). Also available is GEN sumdiv(GEN n), for k = 1.


sqrtint(x)

Returns the integer square root of x, i.e. the largest integer y such that y^2 ≤ x, where x a non-negative integer.

  ? N = 120938191237; sqrtint(N)
  %1 = 347761
  ? sqrt(N)
  %2 = 347761.68741970412747602130964414095216

The library syntax is GEN sqrtint(GEN x).


sqrtnint(x, n)

Returns the integer n-th root of x, i.e. the largest integer y such that y^n ≤ x, where x is a non-negative integer.

  ? N = 120938191237; sqrtnint(N, 5)
  %1 = 164
  ? N^(1/5)
  %2 = 164.63140849829660842958614676939677391

The special case n = 2 is sqrtint

The library syntax is GEN sqrtnint(GEN x, long n).


sumdedekind(h, k)

Returns the Dedekind sum attached to the integers h and k, corresponding to a fast implementation of

    s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))

The library syntax is GEN sumdedekind(GEN h, GEN k).


sumdigits(n, {B = 10})

Sum of digits in the integer n, when written in base B > 1.

  ? sumdigits(123456789)
  %1 = 45
  ? sumdigits(123456789, 2)
  %1 = 16

Note that the sum of bits in n is also returned by hammingweight. This function is much faster than vecsum(digits(n,B)) when B is 10 or a power of 2, and only slightly faster in other cases.

The library syntax is GEN sumdigits0(GEN n, GEN B = NULL). Also available is GEN sumdigits(GEN n), for B = 10.


znchar(D)

Given a datum D describing a group (ℤ/Nℤ)* and a Dirichlet character χ, return the pair [G, chi], where G is znstar(N, 1)) and chi is a GP character.

The following possibilities for D are supported

* a non-zero t_INT congruent to 0,1 modulo 4, return the real character modulo D given by the Kronecker symbol (D/.);

* a t_INTMOD Mod(m, N), return the Conrey character modulo N of index m (see znconreylog).

* a modular form space as per mfinit([N,k,χ]) or a modular form for such a space, return the underlying Dirichlet character χ (which may be defined modulo a divisor of N but need not be primitive).

In the remaining cases, G is initialized by znstar(N, 1).

* a pair [G, chi], where chi is a standard GP Dirichlet character c = (cj) on G (generic character t_VEC or Conrey characters t_COL or t_INT); given generators G = ⨁ (ℤ/djℤ) gj, χ(gj) = e(cj/dj).

* a pair [G, chin], where chin is a normalized representation [n, ~{c}] of the Dirichlet character c; χ(gj) = e(~{c}j / n) where n is minimal (order of χ).

  ? [G,chi] = znchar(-3);
  ? G.cyc
  %2 = [2]
  ? chareval(G, chi, 2)
  %3 = 1/2
  ?  kronecker(-3,2)
  %4 = -1
  ? znchartokronecker(G,chi)
  %5 = -3
  ? mf = mfinit([28, 5/2, Mod(2,7)]); [f] = mfbasis(mf);
  ? [G,chi] = znchar(mf); [G.mod, chi]
  %7 = [7, [2]~]
  ? [G,chi] = znchar(f); chi
  %8 = [28, [0, 2]~]

The library syntax is GEN znchar(GEN D).


zncharconductor(G, chi)

Let G be attached to (ℤ/qℤ)* (as per G = znstar(q, 1)) and chi be a Dirichlet character on (ℤ/qℤ)* (see Section se:dirichletchar or ??character). Return the conductor of chi:

  ? G = znstar(126000, 1);
  ? zncharconductor(G,11)   \\ primitive
  %2 = 126000
  ? zncharconductor(G,1)    \\ trivial character, not primitive!
  %3 = 1
  ? zncharconductor(G,1009) \\ character mod 5^3
  %4 = 125

The library syntax is GEN zncharconductor(GEN G, GEN chi).


znchardecompose(G, chi, Q)

Let N = ∏p pep and a Dirichlet character χ, we have a decomposition χ = ∏p χp into character modulo N where the conductor of χp divides pep; it equals pep for all p if and only if χ is primitive.

Given a znstar G describing a group (ℤ/Nℤ)*, a Dirichlet character chi and an integer Q, return ∏p | (Q,N) χp. For instance, if Q = p is a prime divisor of N, the function returns χp (as a character modulo N), given as a Conrey character (t_COL).

  ? G = znstar(40, 1);
  ? G.cyc
  %2 = [4, 2, 2]
  ? chi = [2, 1, 1];
  ? chi2 = znchardecompose(G, chi, 2)
  %4 = [1, 1, 0]~
  ? chi5 = znchardecompose(G, chi, 5)
  %5 = [0, 0, 2]~
  ? znchardecompose(G, chi, 3)
  %6 = [0, 0, 0]~
  ? c = charmul(G, chi2, chi5)
  %7 = [1, 1, 2]~  \\ t_COL: in terms of Conrey generators !
  ? znconreychar(G,c)
  %8 = [2, 1, 1]   \\ t_VEC: in terms of SNF generators

The library syntax is GEN znchardecompose(GEN G, GEN chi, GEN Q).


znchargauss(G, chi, {a = 1})

Given a Dirichlet character χ on G = (ℤ/Nℤ)* (see znchar), return the complex Gauss sum g(χ,a) = ∑n = 1^N χ(n) e(a n/N)

  ? [G,chi] = znchar(-3); \\ quadratic Gauss sum: I*sqrt(3)
  ? znchargauss(G,chi)
  %2 = 1.7320508075688772935274463415058723670*I
  ? [G,chi] = znchar(5);
  ? znchargauss(G,chi)  \\ sqrt(5)
  %2 = 2.2360679774997896964091736687312762354
  ? G = znstar(300,1); chi = [1,1,12]~;
  ? znchargauss(G,chi) / sqrt(300) - exp(2*I*Pi*11/25)  \\ = 0
  %4 = 2.350988701644575016 E-38 + 1.4693679385278593850 E-39*I
  ? lfuntheta([G,chi], 1)  \\ = 0
  %5 = -5.79[...] E-39 - 2.71[...] E-40*I

The library syntax is GEN znchargauss(GEN G, GEN chi, GEN a = NULL, long bitprec).


zncharinduce(G, chi, N)

Let G be attached to (ℤ/qℤ)* (as per G = znstar(q,1)) and let chi be a Dirichlet character on (ℤ/qℤ)*, given by

* a t_VEC: a standard character on bid.gen,

* a t_INT or a t_COL: a Conrey index in (ℤ/qℤ)* or its Conrey logarithm; see Section se:dirichletchar or ??character.

Let N be a multiple of q, return the character modulo N extending chi. As usual for arithmetic functions, the new modulus N can be given as a t_INT, via a factorization matrix or a pair [N, factor(N)], or by znstar(N,1).

  ? G = znstar(4, 1);
  ? chi = znconreylog(G,1); \\ trivial character mod 4
  ? zncharinduce(G, chi, 80)  \\ now mod 80
  %3 = [0, 0, 0]~
  ? zncharinduce(G, 1, 80)  \\ same using directly Conrey label
  %4 = [0, 0, 0]~
  ? G2 = znstar(80, 1);
  ? zncharinduce(G, 1, G2)  \\ same
  %4 = [0, 0, 0]~
  
  ? chi = zncharinduce(G, 3, G2)  \\ extend the non-trivial character mod 4
  %5 = [1, 0, 0]~
  ? [G0,chi0] = znchartoprimitive(G2, chi);
  ? G0.mod
  %7 = 4
  ? chi0
  %8 = [1]~

Here is a larger example:

  ? G = znstar(126000, 1);
  ? label = 1009;
  ? chi = znconreylog(G, label)
  %3 = [0, 0, 0, 14, 0]~
  ? [G0,chi0] = znchartoprimitive(G, label); \\ works also with 'chi'
  ? G0.mod
  %5 = 125
  ? chi0 \\ primitive character mod 5^3 attached to chi
  %6 = [14]~
  ? G0 = znstar(N0, 1);
  ? zncharinduce(G0, chi0, G) \\ induce back
  %8 = [0, 0, 0, 14, 0]~
  ? znconreyexp(G, %)
  %9 = 1009

The library syntax is GEN zncharinduce(GEN G, GEN chi, GEN N).


zncharisodd(G, chi)

Let G be attached to (ℤ/Nℤ)* (as per G = znstar(N,1)) and let chi be a Dirichlet character on (ℤ/Nℤ)*, given by

* a t_VEC: a standard character on G.gen,

* a t_INT or a t_COL: a Conrey index in (ℤ/qℤ)* or its Conrey logarithm; see Section se:dirichletchar or ??character.

Return 1 if and only if chi(-1) = -1 and 0 otherwise.

  ? G = znstar(8, 1);
  ? zncharisodd(G, 1)  \\ trivial character
  %2 = 0
  ? zncharisodd(G, 3)
  %3 = 1
  ? chareval(G, 3, -1)
  %4 = 1/2

The library syntax is long zncharisodd(GEN G, GEN chi).


znchartokronecker(G, chi, {flag = 0})

Let G be attached to (ℤ/Nℤ)* (as per G = znstar(N,1)) and let chi be a Dirichlet character on (ℤ/Nℤ)*, given by

* a t_VEC: a standard character on bid.gen,

* a t_INT or a t_COL: a Conrey index in (ℤ/qℤ)* or its Conrey logarithm; see Section se:dirichletchar or ??character.

If flag = 0, return the discriminant D if chi is real equal to the Kronecker symbol (D/.) and 0 otherwise. The discriminant D is fundamental if and only if chi is primitive.

If flag = 1, return the fundamental discriminant attached to the corresponding primitive character.

  ? G = znstar(8,1); CHARS = [1,3,5,7]; \\ Conrey labels
  ? apply(t->znchartokronecker(G,t), CHARS)
  %2 = [4, -8, 8, -4]
  ? apply(t->znchartokronecker(G,t,1), CHARS)
  %3 = [1, -8, 8, -4]

The library syntax is GEN znchartokronecker(GEN G, GEN chi, long flag).


znchartoprimitive(G, chi)

Let G be attached to (ℤ/qℤ)* (as per G = znstar(q, 1)) and chi be a Dirichlet character on (ℤ/qℤ)*, of conductor q0 | q.

  ? G = znstar(126000, 1);
  ? [G0,chi0] = znchartoprimitive(G,11)
  ? G0.mod
  %3 = 126000
  ? chi0
  %4 = 11
  ? [G0,chi0] = znchartoprimitive(G,1);\\ trivial character, not primitive!
  ? G0.mod
  %6 = 1
  ? chi0
  %7 = []~
  ? [G0,chi0] = znchartoprimitive(G,1009)
  ? G0.mod
  %4 = 125
  ? chi0
  %5 = [14]~

Note that znconreyconductor is more efficient since it can return χ0 and its conductor q0 without needing to initialize G0. (The price to pay is a more cryptic format and the need to initalize G0 later but the can be done once for all characters of conductor q0.)

The library syntax is GEN znchartoprimitive(GEN G, GEN chi).


znconreychar(G, m)

Given a znstar G attached to (ℤ/qℤ)* (as per G = znstar(q,1)), this function returns the Dirichlet character attached to m ∈ (ℤ/qℤ)* via Conrey's logarithm, which establishes a "canonical" bijection between (ℤ/qℤ)* and its dual.

Let q = ∏p pep be the factorization of q into distinct primes. For all odd p with ep > 0, let gp be the element in (ℤ/qℤ)* which is

* congruent to 1 mod q/pep,

* congruent mod pep to the smallest positive integer that generates (ℤ/p^2ℤ)*.

For p = 2, we let g4 (if 2e2 ≥ 4) and g8 (if furthermore (2e2 ≥ 8) be the elements in (ℤ/qℤ)* which are

* congruent to 1 mod q/2e2,

* g4 = -1 mod 2e2,

* g8 = 5 mod 2e2.

Then the gp (and the extra g4 and g8 if 2e2 ≥ 2) are independent generators of (ℤ/qℤ)*, i.e. every m in (ℤ/qℤ)* can be written uniquely as ∏p gpmp, where mp is defined modulo the order op of gp and p ∈ Sq, the set of prime divisors of q together with 4 if 4 | q and 8 if 8 | q. Note that the gp are in general not SNF generators as produced by znstar whenever ω(q) ≥ 2, although their number is the same. They however allow to handle the finite abelian group (ℤ/qℤ)* in a fast and elegant way. (Which unfortunately does not generalize to ray class groups or Hecke characters.)

The Conrey logarithm of m is the vector (mp)p ∈ Sq, obtained via znconreylog. The Conrey character χq(m,.) attached to m mod q maps each gp, p ∈ Sq to e(mp / op), where e(x) = exp(2iπ x). This function returns the Conrey character expressed in the standard PARI way in terms of the SNF generators G.gen.

  ? G = znstar(8,1);
  ? G.cyc
  %2 = [2, 2]  \\ Z/2 x Z/2
  ? G.gen
  %3 = [7, 3]
  ? znconreychar(G,1)  \\ 1 is always the trivial character
  %4 = [0, 0]
  ? znconreychar(G,2)  \\ 2 is not coprime to 8 !!!
    ***   at top-level: znconreychar(G,2)
    ***                 ^ —  —  —  —  — --
    *** znconreychar: elements not coprime in Zideallog:
      2
      8
    ***   Break loop: type 'break' to go back to GP prompt
  break>
  
  ? znconreychar(G,3)
  %5 = [0, 1]
  ? znconreychar(G,5)
  %6 = [1, 1]
  ? znconreychar(G,7)
  %7 = [1, 0]

We indeed get all 4 characters of (ℤ/8ℤ)*.

For convenience, we allow to input the Conrey logarithm of m instead of m:

  ? G = znstar(55, 1);
  ? znconreychar(G,7)
  %2 = [7, 0]
  ? znconreychar(G, znconreylog(G,7))
  %3 = [7, 0]

The library syntax is GEN znconreychar(GEN G, GEN m).


znconreyconductor(G, chi, {&chi0})

Let G be attached to (ℤ/qℤ)* (as per G = znstar(q, 1)) and chi be a Dirichlet character on (ℤ/qℤ)*, given by

* a t_VEC: a standard character on bid.gen,

* a t_INT or a t_COL: a Conrey index in (ℤ/qℤ)* or its Conrey logarithm; see Section se:dirichletchar or ??character.

Return the conductor of chi, as the t_INT bid.mod if chi is primitive, and as a pair [N, faN] (with faN the factorization of N) otherwise.

If chi0 is present, set it to the Conrey logarithm of the attached primitive character.

  ? G = znstar(126000, 1);
  ? znconreyconductor(G,11)   \\ primitive
  %2 = 126000
  ? znconreyconductor(G,1)    \\ trivial character, not primitive!
  %3 = [1, matrix(0,2)]
  ? N0 = znconreyconductor(G,1009, &chi0) \\ character mod 5^3
  %4 = [125, Mat([5, 3])]
  ? chi0
  %5 = [14]~
  ? G0 = znstar(N0, 1);      \\ format [N,factor(N)] accepted
  ? znconreyexp(G0, chi0)
  %7 = 9
  ? znconreyconductor(G0, chi0) \\ now primitive, as expected
  %8 = 125

The group G0 is not computed as part of znconreyconductor because it needs to be computed only once per conductor, not once per character.

The library syntax is GEN znconreyconductor(GEN G, GEN chi, GEN *chi0 = NULL).


znconreyexp(G, chi)

Given a znstar G attached to (ℤ/qℤ)* (as per G = znstar(q, 1)), this function returns the Conrey exponential of the character chi: it returns the integer m ∈ (ℤ/qℤ)* such that znconreylog(G, m) is chi.

The character chi is given either as a

* t_VEC: in terms of the generators G.gen;

* t_COL: a Conrey logarithm.

  ? G = znstar(126000, 1)
  ? znconreylog(G,1)
  %2 = [0, 0, 0, 0, 0]~
  ? znconreyexp(G,%)
  %3 = 1
  ? G.cyc \\ SNF generators
  %4 = [300, 12, 2, 2, 2]
  ? chi = [100, 1, 0, 1, 0]; \\ some random character on SNF generators
  ? znconreylog(G, chi)  \\ in terms of Conrey generators
  %6 = [0, 3, 3, 0, 2]~
  ? znconreyexp(G, %)  \\ apply to a Conrey log
  %7 = 18251
  ? znconreyexp(G, chi) \\ ... or a char on SNF generators
  %8 = 18251
  ? znconreychar(G,%)
  %9 = [100, 1, 0, 1, 0]

The library syntax is GEN znconreyexp(GEN G, GEN chi).


znconreylog(G, m)

Given a znstar attached to (ℤ/qℤ)* (as per G = znstar(q,1)), this function returns the Conrey logarithm of m ∈ (ℤ/qℤ)*.

Let q = ∏p pep be the factorization of q into distinct primes, where we assume e2 = 0 or e2 ≥ 2. (If e2 = 1, we can ignore 2 from the factorization, as if we replaced q by q/2, since (ℤ/qℤ)* ~ (ℤ/(q/2)ℤ)*.)

For all odd p with ep > 0, let gp be the element in (ℤ/qℤ)* which is

* congruent to 1 mod q/pep,

* congruent mod pep to the smallest positive integer that generates (ℤ/p^2ℤ)*.

For p = 2, we let g4 (if 2e2 ≥ 4) and g8 (if furthermore (2e2 ≥ 8) be the elements in (ℤ/qℤ)* which are

* congruent to 1 mod q/2e2,

* g4 = -1 mod 2e2,

* g8 = 5 mod 2e2.

Then the gp (and the extra g4 and g8 if 2e2 ≥ 2) are independent generators of ℤ/qℤ*, i.e. every m in (ℤ/qℤ)* can be written uniquely as ∏p gpmp, where mp is defined modulo the order op of gp and p ∈ Sq, the set of prime divisors of q together with 4 if 4 | q and 8 if 8 | q. Note that the gp are in general not SNF generators as produced by znstar whenever ω(q) ≥ 2, although their number is the same. They however allow to handle the finite abelian group (ℤ/qℤ)* in a fast and elegant way. (Which unfortunately does not generalize to ray class groups or Hecke characters.)

The Conrey logarithm of m is the vector (mp)p ∈ Sq. The inverse function znconreyexp recovers the Conrey label m from a character.

  ? G = znstar(126000, 1);
  ? znconreylog(G,1)
  %2 = [0, 0, 0, 0, 0]~
  ? znconreyexp(G, %)
  %3 = 1
  ? znconreylog(G,2)  \\ 2 is not coprime to modulus !!!
    ***   at top-level: znconreylog(G,2)
    ***                 ^ —  —  —  —  — --
    *** znconreylog: elements not coprime in Zideallog:
      2
      126000
    ***   Break loop: type 'break' to go back to GP prompt
  break>
  ? znconreylog(G,11) \\ wrt. Conrey generators
  %4 = [0, 3, 1, 76, 4]~
  ? log11 = ideallog(,11,G)   \\ wrt. SNF generators
  %5 = [178, 3, -75, 1, 0]~

For convenience, we allow to input the ordinary discrete log of m, ideallog(,m,bid), which allows to convert discrete logs from bid.gen generators to Conrey generators.

  ? znconreylog(G, log11)
  %7 = [0, 3, 1, 76, 4]~

We also allow a character (t_VEC) on bid.gen and return its representation on the Conrey generators.

  ? G.cyc
  %8 = [300, 12, 2, 2, 2]
  ? chi = [10,1,0,1,1];
  ? znconreylog(G, chi)
  %10 = [1, 3, 3, 10, 2]~
  ? n = znconreyexp(G, chi)
  %11 = 84149
  ? znconreychar(G, n)
  %12 = [10, 1, 0, 1, 1]

The library syntax is GEN znconreylog(GEN G, GEN m).


zncoppersmith(P, N, X, {B = N})

N being an integer and P ∈ ℤ[X], finds all integers x with |x| ≤ X such that gcd(N, P(x)) ≥ B, using Coppersmith's algorithm (a famous application of the LLL algorithm). X must be smaller than exp(log^2 B / (deg(P) log N)): for B = N, this means X < N1/deg(P). Some x larger than X may be returned if you are very lucky. The smaller B (or the larger X), the slower the routine will be. The strength of Coppersmith method is the ability to find roots modulo a general composite N: if N is a prime or a prime power, polrootsmod or polrootspadic will be much faster.

We shall now present two simple applications. The first one is finding non-trivial factors of N, given some partial information on the factors; in that case B must obviously be smaller than the largest non-trivial divisor of N.

  setrand(1); \\ to make the example reproducible
  interval = [10^30, 10^31];
  p = randomprime(interval);
  q = randomprime(interval); N = p*q;
  p0 = p % 10^20; \\ assume we know 1) p > 10^29, 2) the last 19 digits of p
  L = zncoppersmith(10^19*x + p0, N, 10^12, 10^29)
  
  \\ result in 10ms.
  %6 = [738281386540]
  ? gcd(L[1] * 10^19 + p0, N) == p
  %7 = 1

and we recovered p, faster than by trying all possibilities < 1012.

The second application is an attack on RSA with low exponent, when the message x is short and the padding P is known to the attacker. We use the same RSA modulus N as in the first example:

  setrand(1);
  P = random(N);    \\ known padding
  e = 3;            \\ small public encryption exponent
  X = floor(N^0.3); \\ N^(1/e - epsilon)
  x0 = random(X);   \\ unknown short message
  C = lift( (Mod(x0,N) + P)^e ); \\ known ciphertext, with padding P
  zncoppersmith((P + x)^3 - C, N, X)
  
  \\ result in 244ms.
  %14 = [2679982004001230401]
  
  ? %[1] == x0
  %15 = 1

We guessed an integer of the order of 1018, almost instantly.

The library syntax is GEN zncoppersmith(GEN P, GEN N, GEN X, GEN B = NULL).


znlog(x, g, {o})

This functions allows two distinct modes of operation depending on g:

* if g is the output of znstar (with initialization), we compute the discrete logarithm of x with respect to the generators contained in the structure. See ideallog for details.

* else g is an explicit element in (ℤ/Nℤ)*, we compute the discrete logarithm of x in (ℤ/Nℤ)* in base g. The rest of this entry describes the latter possibility.

The result is [] when x is not a power of g, though the function may also enter an infinite loop in this case.

If present, o represents the multiplicative order of g, see Section se:DLfun; the preferred format for this parameter is [ord, factor(ord)], where ord is the order of g. This provides a definite speedup when the discrete log problem is simple:

  ? p = nextprime(10^4); g = znprimroot(p); o = [p-1, factor(p-1)];
  ? for(i=1,10^4, znlog(i, g, o))
  time = 163 ms.
  ? for(i=1,10^4, znlog(i, g))
  time = 200 ms. \\ a little slower

The result is undefined if g is not invertible mod N or if the supplied order is incorrect.

This function uses

* a combination of generic discrete log algorithms (see below).

* in (ℤ/Nℤ)* when N is prime: a linear sieve index calculus method, suitable for N < 1050, say, is used for large prime divisors of the order.

The generic discrete log algorithms are:

* Pohlig-Hellman algorithm, to reduce to groups of prime order q, where q | p-1 and p is an odd prime divisor of N,

* Shanks baby-step/giant-step (q < 232 is small),

* Pollard rho method (q > 232).

The latter two algorithms require O(sqrt{q}) operations in the group on average, hence will not be able to treat cases where q > 1030, say. In addition, Pollard rho is not able to handle the case where there are no solutions: it will enter an infinite loop.

  ? g = znprimroot(101)
  %1 = Mod(2,101)
  ? znlog(5, g)
  %2 = 24
  ? g^24
  %3 = Mod(5, 101)
  
  ? G = znprimroot(2 * 101^10)
  %4 = Mod(110462212541120451003, 220924425082240902002)
  ? znlog(5, G)
  %5 = 76210072736547066624
  ? G^% == 5
  %6 = 1
  ? N = 2^4*3^2*5^3*7^4*11; g = Mod(13, N); znlog(g^110, g)
  %7 = 110
  ? znlog(6, Mod(2,3))  \\ no solution
  %8 = []

For convenience, g is also allowed to be a p-adic number:

  ? g = 3+O(5^10); znlog(2, g)
  %1 = 1015243
  ? g^%
  %2 = 2 + O(5^10)

The library syntax is GEN znlog0(GEN x, GEN g, GEN o = NULL). The function GEN znlog(GEN x, GEN g, GEN o) is also available


znorder(x, {o})

x must be an integer mod n, and the result is the order of x in the multiplicative group (ℤ/nℤ)*. Returns an error if x is not invertible. The parameter o, if present, represents a non-zero multiple of the order of x, see Section se:DLfun; the preferred format for this parameter is [ord, factor(ord)], where ord = eulerphi(n) is the cardinality of the group.

The library syntax is GEN znorder(GEN x, GEN o = NULL). Also available is GEN order(GEN x).


znprimroot(n)

Returns a primitive root (generator) of (ℤ/nℤ)*, whenever this latter group is cyclic (n = 4 or n = 2p^k or n = p^k, where p is an odd prime and k ≥ 0). If the group is not cyclic, the result is undefined. If n is a prime power, then the smallest positive primitive root is returned. This may not be true for n = 2p^k, p odd.

Note that this function requires factoring p-1 for p as above, in order to determine the exact order of elements in (ℤ/nℤ)*: this is likely to be costly if p is large.

The library syntax is GEN znprimroot(GEN n).


znstar(n, {flag = 0})

Gives the structure of the multiplicative group (ℤ/nℤ)*. The output G depends on the value of flag:

* flag = 0 (default), an abelian group structure [h,d,g], where h = φ(n) is the order (G.no), d (G.cyc) is a k-component row-vector d of integers di such that di > 1, di | di-1 for i ≥ 2 and (ℤ/nℤ)* ~ ∏i = 1^k (ℤ/diℤ), and g (G.gen) is a k-component row vector giving generators of the image of the cyclic groups ℤ/diℤ.

* flag = 1 the result is a bid structure; this allows computing discrete logarithms using znlog (also in the non-cyclic case!).

  ? G = znstar(40)
  %1 = [16, [4, 2, 2], [Mod(17, 40), Mod(21, 40), Mod(11, 40)]]
  ? G.no   \\ eulerphi(40)
  %2 = 16
  ? G.cyc  \\ cycle structure
  %3 = [4, 2, 2]
  ? G.gen  \\ generators for the cyclic components
  %4 = [Mod(17, 40), Mod(21, 40), Mod(11, 40)]
  ? apply(znorder, G.gen)
  %5 = [4, 2, 2]

For user convenience, we define znstar(0) as [2, [2], [-1]], corresponding to ℤ*, but flag = 1 is not implemented in this trivial case.

The library syntax is GEN znstar0(GEN n, long flag).