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

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.

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 nonzero 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 pass it to the relevant functions, as in

? p = nextprime(10^30); ? 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

The finite abelian group G = (ℤ/Nℤ)^{*} can be written G = ⨁ _{i ≤
n} (ℤ/d_{i}ℤ) g_{i}, with d_{n} | ... | d_{2} | d_{1} (SNF condition),
all d_{i} > 0, and ∏_{i} d_{i} = φ(N).

The SNF condition makes the d_{i} unique, but the generators g_{i}, of
respective order d_{i}, are definitely not unique. The ⨁ notation
means that all elements of G can be written uniquely as ∏_{i} g_{i}^{ni}
where n_{i} ∈ ℤ/d_{i}ℤ. The g_{i} are the so-called *SNF generators*
of G.

***** a *character* on the abelian group
⨁ (ℤ/d_{j}ℤ) g_{j}
is given by a row vector χ = [a_{1},...,a_{n}] of integers 0 ≤ a_{i} <
d_{i} such that χ(g_{j}) = e(a_{j} / d_{j}) for all j, with the standard
notation e(x) := exp(2iπ x).
In other words,
χ(∏ g_{j}^{nj}) = e(∑ a_{j} n_{j} / d_{j}).

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 remainders instead. If N = ∏ p^{ep} is
the factorization of N into primes, the so-called *Conrey generators*
of G are the generators of the (ℤ/p^{ep}ℤ)^{*} lifted to (ℤ/Nℤ)^{*} by
requesting that they be congruent to 1 modulo N/p^{ep} (for p odd we
take the smallest positive primitive root mod p^2, and for p = 2
we take -1 if
e_{2} > 1 and additionally 5 if e_{2} > 2). We can again write G =
⨁ _{i ≤ n} (ℤ/D_{i}ℤ) G_{i}, where again ∏_{i} D_{i} = φ(N). These
generators don't satisfy the SNF condition in general since their orders are
now (p-1)p^{ep-1} for p odd; for p = 2, the generator -1 has order
2 and 5 has order 2^{e2-2} (e_{2} > 2). Nevertheless, any m ∈
(ℤ/Nℤ)^{*} can be uniquely decomposed as m = ∏ G_{i}^{mi} for some m_{i}
modulo D_{i} and we can define a character by χ(G_{j}) = e(m_{j} / D_{j}) for
all j.

***** The *column vector* of the m_{j}, 0 ≤ m_{j} < D_{j} 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_COL`

s.

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

To sum up a Dirichlet character can be defined by a `t_INTMOD`

`Mod`

(m, N), a `t_INT`

lift (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`

, `t_INTMOD`

and `t_COL`

to a SNF
character.

`znconreylog`

transforms `t_INT`

, `t_INTMOD`

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.

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.

? addprimes(37975227936943673922808872755445627854565536638199) ? factor(15226050279225333605356183781326374297180681149613806\ 88657908494580122963258952897654000350692006139) %2 = [37975227936943673922808872755445627854565536638199 1] [40094690950920881030683735292761468389214899724061 1] ? ## *** last result computed in 0 ms.

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)

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, returns 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, returns [].

The library syntax is `GEN `

.**bestappr**(GEN x, GEN B = NULL)

Using variants of the extended Euclidean algorithm (Padé approximants), 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 nonnegative 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_SER`

, we implicitly
convert the input to a `t_POLMOD`

modulo N = t^k where k is the
series absolute precision.

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

***** If x is a `t_RFRAC`

, we implicitly convert the input to a
`t_POLMOD`

modulo N = t^k where k = 2B + 1. If B was omitted,
we return x:

? T = (4*t^2 + 2*t + 3)/(t+1)^10; ? bestapprPade(T,1) %2 = [] \\ impossible ? bestapprPade(T,2) %3 = 27/(337*t^2 + 84*t + 9) ? bestapprPade(T,3) %4 = (4253*t - 3345)/(-39007*t^3 - 28519*t^2 - 8989*t - 1115)

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)

Deprecated alias for `gcdext`

The library syntax is `GEN `

.**gcdext0**(GEN x, GEN y)

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)

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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_{1},...,a_{n}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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(16, 1); \\ (Z/16Z)^{*}? charconj(G, 3) \\ Conrey label %2 = [1, 1]~ ? znconreyexp(G, %) %3 = 11 \\ attached Conrey label; indeed 11 = 3^(-1) mod 16 ? chi = znconreylog(G, 3); ? charconj(G, chi) \\ Conrey logarithm %5 = [1, 1]~

The library syntax is `GEN `

.
Also available is
**charconj0**(GEN cyc, GEN chi)`GEN `

, when **charconj**(GEN cyc, GEN chi)`cyc`

is known to
be a vector of elementary divisors and `chi`

a compatible character
(no checks).

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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 = [a_{1},...,a_{n}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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 `

.
Also available is
**chardiv0**(GEN cyc, GEN a, GEN b)`GEN `

, when **chardiv**(GEN cyc, GEN a, GEN b)`cyc`

is known to
be a vector of elementary divisors and a, b are compatible characters
(no checks).

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 = (ℤ/o_{1}ℤ) g_{1} ⨁ ...⨁ (ℤ/o_{d}ℤ) g_{d}
for some generators (g_{i}) of respective order d_{i}, where the group has
exponent o := lcm_{i} o_{i}. Since ζ^o = 1, the vector (c_{i}) in
∏ (ℤ/o_{i}ℤ) defines a character χ on G via χ(g_{i}) =
ζ^{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)

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)

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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_{1},...,a_{n}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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 g_{j} 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 `

.
Also available is
**charker0**(GEN cyc, GEN chi)`GEN `

, when **charker**(GEN cyc, GEN chi)`cyc`

is known to
be a vector of elementary divisors and `chi`

a compatible character
(no checks).

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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 = [a_{1},...,a_{n}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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 `

.
Also available is
**charmul0**(GEN cyc, GEN a, GEN b)`GEN `

, when **charmul**(GEN cyc, GEN a, GEN b)`cyc`

is known to
be a vector of elementary divisors and a, b are compatible characters
(no checks).

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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_{1},...,a_{k}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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 `

.
Also available is
**charorder0**(GEN cyc, GEN chi)`GEN `

, when **charorder**(GEN cyc, GEN chi)`cyc`

is known to
be a vector of elementary divisors and `chi`

a compatible character
(no checks).

Let *cyc* represent a finite abelian group by its elementary
divisors, i.e. (d_{j}) represents ∑_{j ≤ k} ℤ/d_{j}ℤ with d_{k}
| ... | d_{1}; 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 = [a_{1},...,a_{n}] such that
χ(∏ g_{j}^{nj}) = exp(2π i∑ a_{j} n_{j} / d_{j}), where g_{j} denotes
the generator (of order d_{j}) 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 `

.
Also available is
**charpow0**(GEN cyc, GEN a, GEN n)`GEN `

, when **charpow**(GEN cyc, GEN a, GEN n)`cyc`

is known to
be a vector of elementary divisors (no check).

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 `

is also available.**chinese1**(GEN x)

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)

Returns the row vector whose components are the partial quotients of the
continued fraction expansion of x. In other words, a result
[a_{0},...,a_{n}] means that x ~ a_{0}+1/(a_{1}+...+1/a_{n}). The
output is normalized so that a_{n} != 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/b_{0})(a_{0}+b_{1}/(a_{1}+...+b_{n}/a_{n})); for a numerical continued fraction
(x real), the a_{i} are integers, as large as possible; if x is a
rational function, they are polynomials with deg a_{i} = deg b_{i} + 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 a_{i}; 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/a_{1} + 1/(a_1a_{2}) +... ),
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 `

.
Also available are **contfrac0**(GEN x, GEN b = NULL, long nmax)`GEN `

,
**gboundcf**(GEN x, long nmax)`GEN `

and **gcf**(GEN x)`GEN `

.**gcf2**(GEN b, GEN x)

When x is a vector or a one-row matrix, x
is considered as the list of partial quotients [a_{0},a_{1},...,a_{n}] of a
rational number, and the result is the 2 by 2 matrix
[p_{n},p_{n-1};q_{n},q_{n-1}] in the standard notation of continued fractions,
so p_{n}/q_{n} = a_{0}+1/(a_{1}+...+1/a_{n}). If x is a matrix with two rows
[b_{0},b_{1},...,b_{n}] and [a_{0},a_{1},...,a_{n}], this is then considered as a
generalized continued fraction and we have similarly
p_{n}/q_{n} = (1/b_{0})(a_{0}+b_{1}/(a_{1}+...+b_{n}/a_{n})). Note that in this case one
usually has b_{0} = 1.

If n ≥ 0 is present, returns all convergents from p_{0}/q_{0} up to
p_{n}/q_{n}. (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 `

.
also available is **contfracpnqn**(GEN x, long n)`GEN `

for n = -1.**pnqn**(GEN x)

If n is an integer written as
n = df^2 with d squarefree, returns d. If *flag* is nonzero,
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 `

.
Also available are **core0**(GEN n, long flag)`GEN `

(**core**(GEN n)*flag* = 0) and
`GEN `

(**core2**(GEN n)*flag* = 1)

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 nonzero integer
n, this routine returns the (unique) fundamental discriminant d
such that n = df^2, f a positive rational number. If *flag* is nonzero,
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 `

.
Also available are **coredisc0**(GEN n, long flag)`GEN `

(**coredisc**(GEN n)*flag* = 0) and
`GEN `

(**coredisc2**(GEN n)*flag* = 1)

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)

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)

x and y being vectors of perhaps different lengths representing
the Dirichlet series ∑_{n} x_{n} n^{-s} and ∑_{n} y_{n} 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 nonzero coefficient.

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

The library syntax is `GEN `

.**dirmul**(GEN x, GEN y)

For positive integer N and complex number x, return the sum
f(1)1^x + f(2)2^x +...+ f(N)N^x, where f is a completely
multiplicative function. If f is omitted, return
1^x +...+ N^x. When N ≤ 0, the function returns 0.
If `both`

is set, return the pair for arguments (x,f) and
(-1-x,f). If `both = 2`

, assume in addition that f is
real-valued.

**Caveat.** when both is set, the present implementation
assumes that |f(n)| is either 0 or 1, which is the case for
Dirichlet characters.

Unlike variants using `dirpowers(N,x)`

, this function uses O(sqrt{N})
memory instead of O(N). And it is faster for large N. The return value
is usually a floating point number, but it will be exact if the result
is an integer. On the other hand, rational numbers are converted to
floating point approximations, since they are likely to blow up for large N.

? dirpowers(5, 2) %1 = [1, 4, 9, 16, 25] ? vecsum(%) %2 = 55 ? dirpowerssum(5, 2) %3 = 55 ? dirpowerssum(5, -2) %4 = 1.4636111111111111111111111111111111111 ? \p200 ? s = 1/2 + I * sqrt(3); N = 10^7; ? dirpowerssum(N, s); time = 11,425 ms. ? vecsum(dirpowers(N, s)) time = 19,365 ms. ? dirpowerssum(N, s, n->kronecker(-23,n)) time = 10,981 ms.

The `dirpowerssum`

commands work with default stack size,
the `dirpowers`

one requires a stacksize of at least 5GB.

The library syntax is

. When f = **dirpowerssumfun**(ulong N, GEN x, void *E, GEN (*f)(void*, ulong, long), long prec)`NULL`

, one may use
`GEN `

.**dirpowerssum**(ulong N, GEN x, long prec)

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]

This function requires a lot of memory if x has many divisors. The following idiom runs through all divisors using very little memory, in no particular order this time:

F = factor(x); P = F[,1]; E = F[,2]; forvec(e = vectorv(#E,i,[0,E[i]]), d = factorback(P,e); ...)

If the factorization of d is also desired, then [P,e] almost provides it but not quite: e may contain 0 exponents, which are not allowed in factorizations. These must be sieved out as in:

? tofact(P,E) = matreduce(Mat([P,E])); ? tofact([2,3,5,7]~, [4,0,2,0]~) %4 = [2 4] [5 2]

We can then run the above loop with `tofact(P,e)`

instead of,
or together with, `factorback`

.

The library syntax is `GEN `

.
The functions **divisors0**(GEN x, long flag)`GEN `

(**divisors**(GEN N)*flag* = 0) and
`GEN `

(**divisors_factored**(GEN N)*flag* = 1) are also available.

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)

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 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" > 2^{64} 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) %2 = [59649589127497217 1] [5704689200685129054721 1] ? isprime( fa[,1] ) %3 = [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 only trial divide
by all primes p < D and the `addprimes`

entries, then skip all
expensive factorization methods. The limit D must be nonnegative.
In this case, one entry in the factorization may be a composite number: all
factors less than D^2 and primes from the `addprimes`

table
are actual primes. But (at most) one entry may not verify this criterion,
and it may be prime or composite: it is only known to be coprime to all
other entries and not a pure power..

? factor(2^2^7 +1, 10^5) %4 = [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. %5 = [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]

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

Unless `factor_proven`

is set, large factors are actually pseudoprimes,
not proven primes; a rational factor is prime if less than 2^{64} and an
irrational one if its norm is less than 2^{64}.

? factor(5*I) %9 = [ 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) %10 = [ I 1] [ 2 + I 1] [1 + 2*I 1] ? factorback(%) %11 = -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) %12 = [t + 2 1] [2*t + 1 1] ? content(P, 1) \\ Q-content %13 = 1/2 ? pollead(factorback(F)) / pollead(P) %14 = 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 composite 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 F_{3}? factor(T, ffgen(ffinit(3,2,'t))^0) \\ over F_{3^2}? factor(T, Mod(Mod(1,3), t^2+t+2)) \\ over F_{3^2}, again ? factor(T, O(3^6)) \\ over Q_{3}, precision 6 ? factor(T, 1.) \\ over R, current precision ? factor(T, I*1.) \\ over C ? factor(T, Mod(1, y^3-2)) \\ over Q(2^{1/3})

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

? factormod(T, 3) \\ over F_{3}? factormod(T, [t^2+t+2, 3]) \\ over F_{3^2}? factormod(T, ffgen(3^2, 't)) \\ over F_{3^2}? factorpadic(T, 3,6) \\ over Q_{3}, precision 6 ? nffactor(y^3-2, T) \\ over Q(2^{1/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_INTMOD`

s and `t_POLMOD`

s (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 F_{3}? factor(T*ffgen(ffinit(3,2,'t))^0) \\ over F_{3^2}? factor(T*Mod(Mod(1,3), t^2+t+2)) \\ over F_{3^2}, again ? factor(T*(1 + O(3^6)) \\ over Q_{3}, 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(2^{1/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 Q_{2}%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[t_{1},..., t_{d}] 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)

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^2 * 3^1 %3 = 12 ? factorback([5,2,3]) %4 = 30

The library syntax is `GEN `

.
Also available is **factorback2**(GEN f, GEN e = NULL)`GEN `

(case e = **factorback**(GEN f)`NULL`

).

This function is obsolete, use factormod.

The library syntax is `GEN `

.**factmod**(GEN x, GEN p)

Obsolete, kept for backward compatibility: use factormod.

The library syntax is `GEN `

.**factorff**(GEN x, GEN p = NULL, GEN a = NULL)

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 `

returns x! as a **mpfact**(long x)`t_INT`

.

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
2^{64}, 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)

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(y) an irreducible polynomial over
𝔽_{p}: factor over 𝔽_{p}[y]/(T) (as usual the main variable of T must have
lower priority than the main variable of f);

***** 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 nonzero 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 F_{3}%1 = [Mod(1, 3)*x^2 + Mod(1, 3) 1] ? liftall( factormod(x^2 + 1, [t^2+1, 3]) ) \\ over F_{9}%2 = [ x + t 1] [x + 2*t 1] \\ same, now letting GP choose a model ? T = ffinit(3,2,'t) %3 = Mod(1, 3)*t^2 + Mod(1, 3)*t + Mod(2, 3) ? liftall( factormod(x^2 + 1, [T, 3]) ) %4 = \\ t is a root of T ! [ x + (t + 2) 1] [x + (2*t + 1) 1] ? t = ffgen(t^2+Mod(1,3)); factormod(x^2 + t^0) \\ same using t_FFELT %5 = [ x + t 1] [x + 2*t 1] ? factormod(x^2+Mod(1,3)) %6 = [Mod(1, 3)*x^2 + Mod(1, 3) 1] ? liftall( factormod(x^2 + Mod(Mod(1,3), y^2+1)) ) %7 = [ x + y 1] [x + 2*y 1]

If *flag* is nonzero, 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)

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.

If f is not squarefree, the result is undefined. 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
the same degree d;

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

The factorization is ordered by increasing degree d of irreducible factors, and the result is obviously canonical. This function is somewhat faster than full factorization.

? f = (x^2 + 1) * (x^2-1); ? factormodSQF(f,3) \\ squarefree over F_{3}%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 a little slower time = 464 ms. ? liftall( factormodDDF(x^2 + 1, [3, t^2+1]) ) \\ over F_{9}%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)

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.

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;

This is somewhat faster than full factorization. The factors are ordered by increasing exponent and the result is obviously canonical.

? f = (x^2 + 1)^3 * (x^2-1)^2; ? factormodSQF(f, 3) \\ over F_{3}%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 F_{9}%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)

Factors n-th cyclotomic polynomial Φ_{n}(x) mod p,
where p is a prime number not dividing n.
Much faster than `factormod(polcyclo(n), p)`

; the irreducible
factors should be identical and given in the same order.
If *single* is set, return a single irreducible factor; else (default)
return all the irreducible factors. Note that repeated calls of this
function with the *single* flag set may return different results because
the algorithm is probabilistic. Algorithms used are as follows.

Let F = ℚ(ζ_{n}). Let K be the splitting field of p in F and e the
conductor of K. Then Φ_{n}(x) and Φ_{e}(x) have the same
number of irreducible factors mod p and there is a simple algorithm
constructing irreducible factors of Φ_{n}(x) from irreducible
factors of Φ_{e}(x). So we may assume n is equal to the
conductor of K.
Let d be the order of p in (ℤ/nℤ)^{ x } and ϕ(n) = df.
Then Φ_{n}(x) has f irreducible factors g_{i}(x) (1 ≤ i ≤ f)
of degree d over 𝔽_{p} or ℤ_{p}.

***** If d is small, then we factor g_{i}(x) into
d linear factors g_{ij}(x), 1 ≤ j ≤ d in 𝔽_{q}[x] (q = p^d) and
construct G_{i}(x) = ∏_{j = 1}^d g_{ij}(x) ∈ 𝔽_{q}[x].
Then G_{i}(x) ∈ 𝔽_{p}[x] and g_{i}(x) = G_{i}(x).

***** If f is small, then we work in K, which is a Galois extension of
degree f over ℚ. The Gaussian period
θ_{k} = Tr_{F/K}(ζ_{n}^{k}) is a sum of k-th power of roots of
g_{i}(x) and K = ℚ(θ_{1}).

Now, for each k, there is a polynomial T_{k}(x) ∈ ℚ[x] satisfying
θ_{k} = T_{k}(θ_{1}) because all θ_{k} are in K. Let T(x) ∈ ℤ[x]
be the minimal polynomial of θ_{1} over ℚ. We get θ_{1} mod p
from T(x) and construct θ_{1},...,θ_{d} mod p using T_{k}(x).
Finally we recover g_{i}(x) from θ_{1},...,θ_{d} by Newton's
formula.

? lift(factormodcyclo(15, 11)) %1 = [x^2 + 9*x + 4, x^2 + 4*x + 5, x^2 + 3*x + 9, x^2 + 5*x + 3] ? factormodcyclo(15, 11, 1) \\ single %2 = Mod(1, 11)*x^2 + Mod(5, 11)*x + Mod(3, 11) ? z1 = lift(factormod(polcyclo(12345),11311)[,1]); time = 32,498 ms. ? z2 = factormodcyclo(12345,11311); time = 47 ms. ? z1 == z2 %4 = 1

The library syntax is `GEN `

where **factormodcyclo**(long n, GEN p, long single, long v = -1)`v`

is a variable number.

Let k, l, m be three finite fields and f a (partial) map from l to m and g a (partial) map from k to l, 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)

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)

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 `

where **ffextend**(GEN a, GEN P, long v = -1)`v`

is a variable number.

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)

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.
For efficiency, the characteristic is not checked to be prime; similarly
if a polynomial is given, we do not check whether it is irreducible.

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.

To obtain a multiplicative generator, call `ffprimroot`

on the result
(which is randomized). Its minimal polynomial then gives a *primitive*
polynomial, which can be used to redefine the finite field so that all
subsequent computations use the new primitive polynomial:

? g = ffgen(16, 't); ? g.mod \\ recover the underlying polynomial. %2 = t^4 + t^3 + t^2 + t + 1 ? g.pol \\ lift g as a t_POL %3 = t ? g.p \\ recover the characteristic %4 = 2 ? fforder(g) \\ g is not a multiplicative generator %5 = 5 ? a = ffprimroot(g) \\ recover a multiplicative generator %6 = t^3 + t^2 + t ? fforder(a) %7 = 15 ? T = minpoly(a) \\ primitive polynomial %8 = Mod(1, 2)*x^4 + Mod(1, 2)*x^3 + Mod(1, 2) ? G = ffgen(T); \\ is now a multiplicative generator ? fforder(G) %10 = 15

The library syntax is `GEN `

where **ffgen**(GEN k, long v = -1)`v`

is a variable number.

To create a generator for a prime finite field, the function
`GEN `

returns **p_to_GEN**(GEN p, long v)`ffgen(p,v)^0`

.

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 `

where **ffinit**(GEN p, long n, long v = -1)`v`

is a variable number.

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 ? ffmap(p, b) %8 = []

The library syntax is `GEN `

.**ffinvmap**(GEN m)

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 \\nota 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)

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)

Given a (partial) map m between two finite fields, express x as an algebraic element over the codomain of m in a way which is compatible with m. The function is applied recursively to the component of vectors, matrices and polynomials.

? a = ffgen([3,5],'a); ? b = ffgen([3,10],'b); ? m = ffembed(a, b); ? mi= ffinvmap(m); ? R = ffmaprel(mi,b) %5 = Mod(b,b^2+(a+1)*b+(a^2+2*a+2))

In particular, this function can be used to compute the relative minimal polynomial, norm and trace:

? minpoly(R) %6 = x^2+(a+1)*x+(a^2+2*a+2) ? trace(R) %7 = 2*a+2 ? norm(R) %8 = a^2+2*a+2

The library syntax is `GEN `

.**ffmaprel**(GEN m, GEN x)

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 `

.
Also available are
**ffnbirred0**(GEN q, long n, long fl)`GEN `

(for **ffnbirred**(GEN q, long n)*flag* = 0)
and `GEN `

(for **ffsumnbirred**(GEN q, long n)*flag* = 1).

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)

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)

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" nontrivial, is to use `polresultant`

: a value
close to 0 means that a small deformation of the inputs has nontrivial 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 `

.
Also available are **ggcd0**(GEN x, GEN y = NULL)`GEN `

, if **ggcd**(GEN x, GEN y)`y`

is not
`NULL`

, and `GEN `

, if **content**(GEN x)`y`

= `NULL`

.

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)

Let inputs x and y be both integers, or both polynomials in the same
variable. Return a vector `[M, [a,b]~]`

, where M is an invertible
2 x 2 matrix such that `M*[x,y]~ = [a,b]~`

, where b is
small. More precisely,

***** polynomial case: det M has degree 0 and we
have deg a ≥ ceil{max(deg x,deg y))/2} > deg b.

***** integer case: det M = ± 1 and we have
a ≥ ceil{sqrt{max(|x|,|y|)}} > b.
Assuming x and y are nonnegative, then M^{-1} has nonnegative
coefficients, and det M is equal to the sign of both main diagonal terms
M[1,1] and M[2,2].

The library syntax is `GEN `

.**ghalfgcd**(GEN x, GEN y)

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)

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)

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)

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 `

.
Also available is
**ispower**(GEN x, GEN k = NULL, GEN *n = NULL)`long `

(k omitted).**gisanypower**(GEN x, GEN *pty)

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)

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)

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)

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 V_{n} is ± 2, without computing U_{n}.)

There are no known composite numbers passing the above test, although it is
expected that infinitely many such numbers exist. In particular, all
composites ≤ 2^{64} 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)

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 ≤ 2^{64} the function returns k > 0 if and
only if n is indeed prime. When n > 2^{64} 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)

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 nonnegative `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_POLMOD`

s of `t_INTMOD`

s
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 `

.
Also available is **issquareall**(GEN x, GEN *n = NULL)`long `

. Deprecated
GP-specific functions **issquare**(GEN x)`GEN `

and
**gissquare**(GEN x)`GEN `

return **gissquareall**(GEN x, GEN *pt)`gen`

and _{0}`gen`

instead of a boolean value._{1}

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

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

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

Return the largest non-negative integer e so that b^e ≤ x, where b > 1 is an integer and x ≥ 1 is a real number. 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 ? logint(Pi^2, 2, &z) %5 = 3 ? z %6 = 8

The number of digits used to write x in base b 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 μ-function of |x|; x must be a nonzero integer.

The library syntax is `long `

.**moebius**(GEN 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`

.

? nextprime(2) %1 = 2 ? precprime(Pi) %2 = 5 ? nextprime(-10) %3 = 2 \\ primes are positive

Despite the name, please note that the function is not guaranteed to return
a prime number, although no counter-example is known at present. The return
value *is* a guaranteed prime if x ≤ 2^{64}. To rigorously prove
that the result is prime in all cases, use `isprime`

.

The library syntax is `GEN `

.**nextprime**(GEN x)

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

The library syntax is `GEN `

.**numdiv**(GEN 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)

Finds the largest pseudoprime (see `ispseudoprime`

) less than or equal
to x; the input 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.

? precprime(2) %1 = 2 ? precprime(Pi) %2 = 3 ? precprime(-10) %3 = 0 \\ primes are positive

The function name comes from *prec*eding *prime*.
Despite the name, please note that the function is not guaranteed to return
a prime number (although no counter-example is known at present); the return
value *is* a guaranteed prime if x ≤ 2^{64}. To rigorously prove
that the result is prime in all cases, use `isprime`

.

The library syntax is `GEN `

.**precprime**(GEN x)

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 10^{11}; make sure to start gp with
`primelimit`

at least sqrt{p_{n}}, e.g. the value
sqrt{nlog (nlog n)} is guaranteed to be sufficient.

The library syntax is `GEN `

.**prime**(long n)

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)

If *flag* = 0 and *partial* > 0, return a (potentially) partial
ECPP certificate.

A PARI ECPP Primality Certificate for the prime N is either a prime
integer N < 2^{64} or a vector `C`

of length ℓ whose ith
component `C[i]`

is a vector [N_{i}, t_{i}, s_{i}, a_{i}, P_{i}] of length 5
where N_{1} = N. It is said to be *valid* if for each
i = 1,..., ℓ, all of the following conditions are satisfied

***** N_{i} is a positive integer

***** t_{i} is an integer such that t_{i}^2 < 4N_{i}

***** s_{i} is a positive integer which divides m_{i} where
m_{i} = N_{i} + 1 - t_{i}

***** If we set q_{i} = (m_{i})/(s_{i}), then

***** q_{i} > (N_{i}^{1/4}+1)^2

***** q_{i} = N_{i+1} if 1 ≤ i < l

***** q_ℓ ≤ 2^{64} is prime

***** a_{i} is an integer

***** `P[i]`

is a vector of length 2 representing the affine
point P_{i} = (x_{i}, y_{i}) on the elliptic curve E: y^2 = x^3 + a_ix + b_{i}
modulo N_{i} where b_{i} = y_{i}^2 - x_{i}^3 - a_ix_{i} satisfying the following:

***** m_{i} P_{i} = oo

***** s_{i} P_{i} != oo

Using the following theorem, the data in the vector `C`

allows to
recursively certify the primality of N (and all the q_{i}) under the single
assumption that q_ℓ be prime.

**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 > (N^{1/4} + 1)^2, q is a prime divisor of m,
mP = oo and (m)/(q)P != oo , then N is prime.

A partial certificate is identical except that the condition q_ℓ ≤
2^{64} is replaced by q_ℓ ≤ 2^{partial}.
Such partial certificate C can be extended to a full certificate by calling
C = primecert(C), or to a longer partial certificate by calling
C = primecert(C,,b) with b < partial.

? 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 ? N = 2^2000+841; ? Cp1 = primecert(N,,1500); \\ partial certificate time = 16,018 ms. ? Cp2 = primecert(Cp1,,1000); \\ (longer) partial certificate time = 5,890 ms. ? C = primecert(Cp2); \\ full certificate for N time = 1,777 ms. ? primecertisvalid(C) %9 = 1 ? primecert(N); time = 23,625 ms.

As the last command shows, attempting a succession of partial certificates should be about as fast as a direct computation.

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 < 2^{64} or a pair [N, C], where C is a vector with ℓ
elements which are either a single integer p_{i} < 2^{64} or a
triple [p_{i},a_{i},C_{i}] with p_{i} > 2^{64} satisfying the following
properties:

***** p_{i} is a prime divisor of N - 1;

***** a_{i} is an integer such that a_{i}^{N-1} = 1 (mod N) and
a_{i}^{(N-1)/pi} - 1 is coprime with N;

***** C_{i} is an N-1 Primality Certificate for p_{i}

***** The product F of the p_{i}^{vpi(N-1)} is strictly larger than
N^{1/3}. Provided that all p_{i} are indeed primes, this implies that any
divisor of N is congruent to 1 modulo F.

***** The Brillhart-Lehmer-Selfridge criterion is satisfied: when we write
N = 1 + c_{1} F + c_{2} F^2 in base F the polynomial 1 + c_{1} X + c_{2} X^2
is irreducible over ℤ, i.e. c_{1}^2 - 4c_{2} 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 ≤ p^{2/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 `

.
Also available is
**primecert0**(GEN N, long flag, long partial)`GEN `

(**ecpp0**(GEN N, long partial)*flag* = 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. See `??primecert`

for the
meaning of the successive N, t, s, a, m, q, E, P. The integer D is the
negative fundamental discriminant `coredisc`

(t^2 - 4N).

***** 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 a = 0 D = -3 m = 99999999999999999453132088964547996 q = 33743830764501150277 E = [0, 1] P = [21567861682493263464353543707814204, 49167839501923147849639425291163552] [2] N = 33743830764501150277 t = -11610830419 s = 734208843 a = 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 ? primecertexport(cert) *** at top-level: primecertexport(cert) *** ^ — — — — — — — *** primecertexport: sorry, N-1 certificate is not yet implemented.

The library syntax is `GEN `

.**primecertexport**(GEN cert, long format)

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)

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)

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)

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| < 10^{18}.

**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.10^{10}, 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| ∈ [10^{10}, 10^{18}], 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 && quadunitnorm(D) > 0, 2, 1) ? QFBclassno(136) %1 = 4 ? qfbclassno(136) %2 = 2 ? quadunitnorm(136) %3 = 1 ? bnfnarrow(bnfinit(x^2 - 136)).cyc %4 = [4] \\ narrow class group is cyclic ~ Z/4Z

Note that the use of `bnfnarrow`

above is only valid
because 136 is a fundamental discriminant: that function is asymptotically
faster (and returns the group structure, not only its order) but only supports
*maximal* orders.
Here are a few more 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)

composition of the binary quadratic forms x and y, with reduction of the result.

The library syntax is `GEN `

.**qfbcomp**(GEN x, GEN 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)

Solve the equation x^2 + dy^2 = n in integers x and y, where
d > 0 and n is prime. Returns the empty vector `[]`

when no solution
exists. It is also allowed to try n = 4 times a prime but the answer is
then guaranteed only if d is 3 mod 4; more precisely if d != 3 mod
4, the algorithm may fail to find a non-primitive solution.

This function is a special case of `qfbsolve`

applied to the principal
form in the imaginary quadratic order of discriminant -4d (returning the
solution with non-negative x and y). As its name implies,
`qfbcornacchia`

uses Cornacchia's algorithm and runs in time quasi-linear
in log n (using `halfgcd`

); in practical ranges, `qfbcornacchia`

should be about twice faster than `qfbsolve`

unless we indicate to the
latter that its second argument is prime (see below).

? qfbcornacchia(1, 113) %1 = [8, 7] ? qfbsolve(Qfb(1,0,1), 113) %2 = [7, -8] ? qfbcornacchia(1, 4*113) \\ misses the non-primitive solution 2*[8,7] %3 = [] ? qfbcornacchia(1, 4*109) \\ finds a non-primitive solution %4 = [20, 6] ? p = 122838793181521; isprime(p) %5 = 1 ? qfbcornacchia(24, p) %6 = [10547339, 694995] ? Q = Qfb(1,0,24); qfbsolve(Q,p) %7 = [10547339, 694995] ? for (i=1, 10^5, qfbsolve(Q, p)) time = 345 ms. ? for (i=1, 10^5, qfbcornacchia(24,p)) \\ faster time = 251 ms. ? for (i=1, 10^5, qfbsolve(Q, Mat([p,1]))) \\ just as fast time = 251 ms.

We used `Mat([p,1])`

to indicate that p^1
was the integer factorization of p, i.e., that p is prime. Without it,
`qfbsolve`

further attempts to factor p and waste a little time.

The library syntax is `GEN `

.**qfbcornacchia**(GEN d, GEN n)

Hurwitz class number of x, when
x is nonnegative 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)

composition of the primitive positive
definite binary quadratic forms x and y (type `t_QFB`

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

.
Also available is **nucomp**(GEN x, GEN y, GEN L)`GEN `

when x = y.**nudupl**(GEN x, GEN 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 ~ 2^{45}.

The library syntax is `GEN `

.**nupow**(GEN x, GEN n, GEN L = NULL)

n-th power of the binary quadratic form
x, computed with reduction (i.e. using `qfbcomp`

).

The library syntax is `GEN `

.**qfbpow**(GEN x, GEN n)

n-th power of the binary quadratic form
x, computed without doing any reduction (i.e. using `qfbcompraw`

).
Here n must be nonnegative and n < 2^{31}.

The library syntax is `GEN `

.**qfbpowraw**(GEN x, long n)

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_QFB`

are not implemented.)

The library syntax is `GEN `

.**primeform**(GEN x, GEN p)

Reduces the binary quadratic form x (updating Shanks's distance function
d if x = [q,d] is and extended indefinite form).
If *flag* is 1, the function performs a single reduction step, and
a complete reduction otherwise.

The arguments *isd*, *sd*, if present, supply the values of
floor{sqrt{D}}, and sqrt{D} respectively, where D
is the discriminant (this is not checked).
If d < 0 these values are useless.

The library syntax is `GEN `

.
Also available is **qfbred0**(GEN x, long flag, GEN isd = NULL, GEN sd = NULL)`GEN `

(**qfbred**(GEN x)*flag* is 0, `isd`

and `sd`

are `NULL`

)

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; *isD*, if
present, must be equal to `sqrtint`

(D), where D > 0 is the
discriminant of x.

The library syntax is `GEN `

.**qfbredsl2**(GEN x, GEN isD = NULL)

Solve the equation Q(x,y) = n in coprime integers x and y (primitive solutions), where Q is a binary quadratic form and n an integer, up to the action of the special orthogonal group G = SO(Q,ℤ), which is isomorphic to the group of units of positive norm of the quadratic order of discriminant D = disc Q. If D > 0, G is infinite. If D < -4, G is of order 2, if D = -3, G is of order 6 and if D = -4, G is of order 4.

Binary digits of *flag* mean:
1: return all solutions if set, else a single solution; return [] if
a single solution is wanted (bit unset) but none exist.
2: also include imprimitive solutions.

When *flag* = 2 (return a single solution, possibly imprimitive), the
algorithm returns a solution with minimal content; in particular, a
primitive solution exists if and only if one is returned.

The integer n can also be given by its factorization matrix

or by the pair [n, *fa* = factor(n)*fa*].

? qfbsolve(Qfb(1,0,2), 603) \\ a single primitive solution %1 = [5, 17] ? qfbsolve(Qfb(1,0,2), 603, 1) \\ all primitive solutions %2 = [[5, 17], [-19, -11], [19, -11], [5, -17]] ? qfbsolve(Qfb(1,0,2), 603, 2) \\ a single, possibly imprimitive solution %3 = [5, 17] \\ actually primitive ? qfbsolve(Qfb(1,0,2), 603, 3) \\ all solutions %4 = [[5, 17], [-19, -11], [19, -11], [5, -17], [-21, 9], [-21, -9]] ? N = 2^128+1; F = factor(N); ? qfbsolve(Qfb(1,0,1),[N,F],1) %3 = [[-16382350221535464479,8479443857936402504], [18446744073709551616,-1],[-18446744073709551616,-1], [16382350221535464479,8479443857936402504]]

For fixed Q, assuming the factorisation of n is given, the algorithm
runs in probabilistic polynomial time in log p, where p is the largest
prime divisor of n, through the computation of square roots of D modulo
4 p). The dependency on Q is more complicated: polynomial time in log
|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
assuming the GRH.

The library syntax is `GEN `

.**qfbsolve**(GEN Q, GEN n, long flag)

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
`quadregulator`

when D < -10^{25}, D > 10^{10}, 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 [c_{1}, c_{2}],
where c_{1} ≤ c_{2} are nonnegative 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 c_{1} means that relations are easier
to find, but more relations are needed and the linear algebra will be
harder. The default value for c_{1} is 0 and means that it is taken equal
to c_{2}. The parameter c_{2} 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
≤ c_{1} (log |D|)^2, then prove that ideals of norm
≤ c_{2} (log |D|)^2 do
not generate a larger group. By default an optimal c_{2} is chosen, so that
the result is provably correct under the GRH — a result of Grenié
and Molteni states that c_{2} = 23/6 ~ 3.83 is fine (and even
c_{2} = 15/4 ~ 3.75 for large |D| > 2.41 E8). But it is possible
to improve on this algorithmically. You may provide a smaller c_{2}, it will
be ignored (we use the provably correct one); you may provide a larger c_{2}
than the default value, which results in longer computing times for equally
correct outputs (under GRH).

The library syntax is `GEN `

.
If you really need to experiment with the **quadclassunit0**(GEN D, long flag, GEN tech = NULL, long prec)*tech* parameter,
it will be more convenient to use
`GEN `

.**Buchquad**(GEN D, double c1, double c2, long prec)

Discriminant of the étale algebra ℚ(sqrt{x}), where x ∈ ℚ^{*}.
This is the same as `coredisc`

(d) where d is the integer
squarefree 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)

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 `

where **quadgen0**(GEN D, long v = -1)`v`

is a variable number.

When *v* does not matter, the function
`GEN `

is also available.**quadgen**(GEN 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)

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 `

where **quadpoly0**(GEN D, long v = -1)`v`

is a variable number.

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)

Regulator of the quadratic order of positive discriminant D in time
Õ(D^{1/2}) using the continued fraction algorithm. Raise
an error if D is not a discriminant (fundamental or not) or if D is a
square. The function `quadclassunit`

is asymptotically faster (and also
in practice for D > 10^{10} or so) but depends on the GRH.

The library syntax is `GEN `

.**quadregulator**(GEN D, long prec)

A fundamental unit u of the real quadratic order
of discriminant D. The integer D must be congruent to 0 or 1 modulo 4
and not a square; the result is a quadratic number (see Section se:quadgen).
If D is not a fundamental discriminant, the algorithm is wasteful: if D =
df^2 with d fundamental, it will be faster to compute `quadunit`

(d)
then raise it to the power `quadunitindex`

(d,f); or keep it in
factored form.

If *v* is given, the variable name is used to display u
else 'w' is used. The algorithm computes the continued fraction
of (1 + sqrt{D}) / 2 or sqrt{D}/2 (see GTM 138, algorithm 5.7.2).
Although the continued fraction length is only O(sqrt{D}),
the function still runs in time Õ(D), in part because the
output size is not polynomially bounded in terms of log D.
See `bnfinit`

and `bnfunits`

for a better alternative for large
D, running in time subexponential in log D and returning the
fundamental units in compact form (as a short list of S-units of size
O(log D)^3 raised to possibly large exponents).

The library syntax is `GEN `

where **quadunit0**(GEN D, long v = -1)`v`

is a variable number.

When *v* does not matter, the function
`GEN `

is also available.**quadunit**(GEN D)

Given a fundamental discriminant D, return the index of the unit group
of the order of conductor f in the units of ℚ(sqrt{D}). This function
uses the continued fraction algorithm and has O(D^{1/2 + ϵ}
f^ϵ) complexity; `quadclassunit`

is asymptotically faster but
depends on the GRH.

? quadunitindex(-3, 2) %1 = 3 ? quadunitindex(5, 2^32) \\ instantaneous %2 = 3221225472 ? quadregulator(5 * 2^64) / quadregulator(5) time = 3min, 1,488 ms. %3 = 3221225472.0000000000000000000000000000

The conductor f can be given in factored form or as
[f, `factor`

(f)]:

? quadunitindex(5, [100, [2,2;5,2]]) %4 = 150 ? quadunitindex(5, 100) %5 = 150 ? quadunitindex(5, [2,2;5,2]) %6 = 150

If D is not fundamental, the result is undefined; you may use the following script instead:

index(d, f) = { my([D,F] = coredisc(d, 1)); quadunitindex(D, f * F) / quadunitindex(D, F) } ? index(5 * 10^2, 10) %7 = 10

The library syntax is `GEN `

.**quadunitindex**(GEN D, GEN f)

Return the norm (1 or -1) of the fundamental unit of the quadratic
order of discriminant D. The integer D must be congruent to 0 or 1
modulo 4 and not a square. This is of course equal to `norm(quadunit(D))`

but faster.

? quadunitnorm(-3) \\ the result is always 1 in the imaginary case %1 = 1 ? quadunitnorm(5) %2 = -1 ? quadunitnorm(17345) %3 = -1 ? u = quadunit(17345) %4 = 299685042291 + 4585831442*w ? norm(u) %5 = -1

This function computes the parity of the continued fraction
expansion and runs in time Õ(D^{1/2}). If D is fundamental,
the function `bnfinit`

is asymptotically faster but depends of the GRH.
If D = df^2 is not fundamental, it will usually be faster to first compute
`quadunitindex`

(d, f). If it is even, the result is 1, else the result
is `quadunitnorm`

(d). The narrow class number of the order of
discriminant D is equal to the class number if the unit norm is 1 and to
twice the class number otherwise.

**Important remark.** Assuming GRH, using `bnfinit`

is *much*
faster, running in time subexponential in log D (instead of exponential
for `quadunitnorm`

). We give examples for the maximal order:

? GRHunitnorm(bnf) = vecprod(bnfsignunit(bnf)[,1]) ? bnf = bnfinit(x^2 - 17345, 1); GRHunitnorm(bnf) %2 = -1 ? bnf = bnfinit(x^2 - nextprime(2^60), 1); GRHunitnorm(bnf) time = 119 ms. %3 = -1 ? quadunitnorm(nextprime(2^60)) time = 24,086 ms. %4 = -1

Note that if the result is -1, it is unconditional because
(if GRH is false) it could happen that our tentative fundamental unit in
*bnf* is actually a power u^k of the true fundamental unit, but we
would still have Norm(u) = -1 (and k odd). We can also remove the
GRH assumption when the result is 1 with a little more work:

? v = bnfunits(bnf)[1][1] \\ a unit in factored form ? v[,2] %= 2; ? nfeltissquare(bnf, nffactorback(bnf, v)) %7 = 0

Under GRH, we know that v is the fundamental unit, but as
above it can be a power u^k of the true fundamental unit u. But the
final two lines prove that v is not a square, hence k is odd and
Norm(u) must also be 1. We modified the factorization matrix
giving v by reducing all exponents modulo 2: this allows to computed
`nffactorback`

even when the factorization involves huge exponents.
And of course the new v is a square if and only if the original one was.

The library syntax is `long `

.**quadunitnorm**(GEN D)

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`

). If `ell`

is 16, 18, 20, 22,
or 26, same for the newform of level 1 and corresponding weight. Otherwise,
compute the coefficient of the trace form at n.
The complexity is in Õ(n^{1/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}
and more generally, setting u = ℓ - 13 and C = 2/ζ(-u) for ℓ
> 12,
∑τ_{ℓ}(n)q^n = q ∏_{n ≥ 1}
(1-q^n)^{24} ( 1 + C∑_{n ≥ 1}n^u q^n / (1-q^n))
produces them in time Õ(N), against Õ(N^{3/2}) for
individual calls to `ramanujantau`

; of course the space complexity then
becomes Õ(N). For other values of ℓ,
`mfcoefs(mftraceform([1,ell]),N)`

is much faster.

? 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, long ell)

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 2^{64} the pseudo prime returned is a proven prime.

? randomprime(100) %1 = 71 ? randomprime([3,100]) %2 = 61 ? randomprime([1,1]) *** at top-level: randomprime([1,1]) *** ^ — — — — — — *** randomprime: domain error in randomprime: *** floor(b) - max(ceil(a),2) < 0 ? randomprime([24,28]) \\ infinite loop

If the optional parameter q is an integer, return a prime congruent to 1 mod q; if q is an intmod, return a prime in the given congruence class. If the class contains no prime in the given interval, the function will raise an exception if the class is not invertible, else run into an infinite loop

? randomprime(100, 4) \\ 1 mod 4 %1 = 71 ? randomprime(100, 4) %2 = 13 ? randomprime([10,100], Mod(2,5)) %3 = 47 ? randomprime(100, Mod(0,2)) \\ silly but works %4 = 2 ? randomprime([3,100], Mod(0,2)) \\ not invertible *** at top-level: randomprime([3,100],Mod(0,2)) *** ^ — — — — — — — — — -- *** randomprime: elements not coprime in randomprime: 0 2 ? randomprime(100, 97) \\ infinite loop

The library syntax is `GEN `

.
Also available is **randomprime0**(GEN N = NULL, GEN q = NULL)`GEN `

.**randomprime**(GEN N = NULL)

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)

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 `

.
Also available is **sumdivk**(GEN x, long k)`GEN `

, for k = 1.**sumdiv**(GEN n)

Returns the integer square root of x, i.e. the largest integer y such that y^2 ≤ x, where x a nonnegative real number. If r is present, set it to the remainder r = x - y^2, which satisfies 0 ≤ r < 2y + 1. Further, when x is an integer, r is an integer satisfying 0 ≤ r ≤ 2y.

? x = 120938191237; sqrtint(x) %1 = 347761 ? sqrt(x) %2 = 347761.68741970412747602130964414095216 ? y = sqrtint(x, &r); r %3 = 478116 ? x - y^2 %4 = 478116 ? sqrtint(9/4, &r) \\ not 3/2 ! %5 = 1 ? r %6 = 5/4

The library syntax is `GEN `

.
Also available is **sqrtint0**(GEN x, GEN *r = NULL)`GEN `

.**sqrtint**(GEN a)

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

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

The special case n = 2 is `sqrtint`

The library syntax is `GEN `

.**sqrtnint**(GEN x, long n)

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)

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 `

.
Also available is **sumdigits0**(GEN n, GEN B = NULL)`GEN `

, for B = 10.**sumdigits**(GEN n)

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 nonzero `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 = (c_{j}) on `G`

(generic character `t_VEC`

or
Conrey characters `t_COL`

or `t_INT`

); given
generators G = ⨁ (ℤ/d_{j}ℤ) g_{j}, χ(g_{j}) = e(c_{j}/d_{j}).

***** a pair `[G, chin]`

, where `chin`

is a *normalized*
representation [n, ~{c}] of the Dirichlet character c; χ(g_{j})
= 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)

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)

Let N = ∏_{p} p^{ep} and a Dirichlet character χ,
we have a decomposition χ = ∏_{p} χ_{p} into character modulo N
where the conductor of χ_{p} divides p^{ep}; it equals p^{ep} 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)

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)

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

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)

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)

Let *G* be attached to (ℤ/qℤ)^{*} (as per
`G = znstar(q, 1)`

) and `chi`

be a Dirichlet character on
(ℤ/qℤ)^{*}, of conductor q_{0} | 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 q_{0} without needing to initialize
G_{0}. The price to pay is a more cryptic format and the need to
initalize G_{0} later, but that needs to be done only once for all characters
with conductor q_{0}.

The library syntax is `GEN `

.**znchartoprimitive**(GEN G, GEN chi)

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} p^{ep} be the factorization of q into distinct primes.
For all odd p with e_{p} > 0, let g_{p} be the element in (ℤ/qℤ)^{*}
which is

***** congruent to 1 mod q/p^{ep},

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

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

***** congruent to 1 mod q/2^{e2},

***** g_{4} = -1 mod 2^{e2},

***** g_{8} = 5 mod 2^{e2}.

Then the g_{p} (and the extra g_{4} and g_{8} if 2^{e2} ≥ 2) are
independent generators of (ℤ/qℤ)^{*}, i.e. every m in (ℤ/qℤ)^{*} can be
written uniquely as ∏_{p} g_{p}^{mp}, where m_{p} is defined modulo the
order o_{p} of g_{p} and p ∈ S_{q}, the set of prime divisors of q
together with 4 if 4 | q and 8 if 8 | q. Note that the g_{p}
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 (m_{p})_{p ∈ Sq}, obtained
via `znconreylog`

. The Conrey character χ_{q}(m,.) attached to
m mod q maps
each g_{p}, p ∈ S_{q} to e(m_{p} / o_{p}), 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)

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)

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)

Given a *znstar* attached to (ℤ/qℤ)^{*} (as per
`G = znstar(q,1)`

), this function returns the Conrey logarithm of
m ∈ (ℤ/qℤ)^{*}.

Let q = ∏_{p} p^{ep} be the factorization of q into distinct primes,
where we assume e_{2} = 0 or e_{2} ≥ 2. (If e_{2} = 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 e_{p} > 0, let g_{p} be the element in (ℤ/qℤ)^{*}
which is

***** congruent to 1 mod q/p^{ep},

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

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

***** congruent to 1 mod q/2^{e2},

***** g_{4} = -1 mod 2^{e2},

***** g_{8} = 5 mod 2^{e2}.

Then the g_{p} (and the extra g_{4} and g_{8} if 2^{e2} ≥ 2) are
independent generators of ℤ/qℤ^{*}, i.e. every m in (ℤ/qℤ)^{*} can be
written uniquely as ∏_{p} g_{p}^{mp}, where m_{p} is defined modulo the
order o_{p} of g_{p} and p ∈ S_{q}, the set of prime divisors of q
together with 4 if 4 | q and 8 if 8 | q. Note that the g_{p}
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 (m_{p})_{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)

Coppersmith's algorithm. N being an integer and P ∈ ℤ[t],
finds in polynomial time in log(N) and d = deg(P) all integers x
with |x| ≤ X such that
gcd(N, P(x)) ≥ B.
This is a famous application of the LLL algorithm meant to help in the
factorization of N. Notice that P may be reduced modulo Nℤ[t] without
affecting the situation. The parameter X must not be too large: assume for
now that the leading coefficient of P is coprime to N, then we must have
d log X log N < log^2 B, i.e., X < N^{1/d} when B = N. Let now
P_{0} be the gcd of the leading coefficient of P and N. In applications to
factorization, we should have P_{0} = 1; otherwise, either P_{0} = N and we can
reduce the degree of P, or P_{0} is a non trivial factor of N. For
completeness, we nevertheless document the exact conditions that X must
satisfy in this case: let p := log_{N} P_{0}, b := log_{N} B, x := log_{N}
X, then

***** either p ≥ d / (2d-1) is large and we must have x d < 2b - 1;

***** or p < d / (2d-1) and we must have both p < b < 1 - p + p/d
and x(d + p(1-2d)) < (b - p)^2. Note that this reduces to
x d < b^2 when p = 0, i.e., the condition described above.

Some x larger than X may be returned if you are
very lucky. The routine runs in polynomial time in log N and d
but the smaller B, or the larger X, the slower.
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 nontrivial factors of N, given some partial information on the factors; in that case B must obviously be smaller than the largest nontrivial divisor of N.

setrand(1); \\ to make the example reproducible [a,b] = [10^30, 10^31]; D = 20; p = randomprime([a,b]); q = randomprime([a,b]); N = p*q; \\ assume we know 0) p | N; 1) p in [a,b]; 2) the last D digits of p p0 = p % 10^D; ? L = zncoppersmith(10^D*x + p0, N, b \ 10^D, a) time = 1ms. %6 = [738281386540] ? gcd(L[1] * 10^D + p0, N) == p %7 = 1

and we recovered p, faster than by trying all
possibilities x < 10^{11}.

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 10^{18}, almost instantly.

The library syntax is `GEN `

.**zncoppersmith**(GEN P, GEN N, GEN X, GEN B = NULL)

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 < 10^{50}, 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 < 2^{32} is small),

***** Pollard rho method (q > 2^{32}).

The latter two algorithms require O(sqrt{q}) operations in the group on
average, hence will not be able to treat cases where q > 10^{30}, 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 `

.
The function
**znlog0**(GEN x, GEN g, GEN o = NULL)`GEN `

is also available**znlog**(GEN x, GEN g, GEN 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 nonzero
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)

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 function will raise an
exception. 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)

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 d_{i} such that d_{i} > 1,
d_{i} | d_{i-1} for i ≥ 2 and
(ℤ/nℤ)^{*} ~ ∏_{i = 1}^k (ℤ/d_{i}ℤ),
and g (`G.gen`

) is a k-component row vector giving generators of
the image of the cyclic groups ℤ/d_{i}ℤ.

***** *flag* = 1 the result is a `bid`

structure;
this allows computing discrete logarithms using `znlog`

(also in the
noncyclic 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)

Finds a minimal set of generators for the subgroup of (ℤ/fℤ)^{*}
given by a vector (or vectorsmall) H of length f:
for 1 ≤ a ≤ f, `H[a]`

is 1 or 0 according as a ∈ H_{F}
or a ∉ H_{F}. In most PARI functions, subgroups of an abelian group
are given as HNF left-divisors of a diagonal matrix, representing the
discrete logarithms of the subgroup generators in terms of a fixed
generators for the group cyclic components. The present function
allows to convert an enumeration of the subgroup elements to this
representation as follows:

? G = znstar(f, 1); ? v = znsubgroupgenerators(H); ? subHNF(G, v) = mathnfmodid(Mat([znlog(h, G) | h<-v]), G.cyc);

The function `subHNF`

can be applied to any
elements of (ℤ/fℤ)^{*}, yielding the subgroup they generate, but using
`znsubgroupgenerators`

first allows to reduce the number of discrete
logarithms to be computed.

For example, if H = {1,4,11,14} ⊂ (ℤ/15ℤ)^{ x },
then we have

? f = 15; H = vector(f); H[1]=H[4]=H[11]=H[14] = 1; ? v = znsubgroupgenerators(H) %2 = [4, 11] ? G = znstar(f, 1); G.cyc %3 = [4, 2] ? subHNF(G, v) %4 = [2 0] [0 1] ? subHNF(G, [1,4,11,14]) %5 = [2 0] [0 1]

This function is mostly useful when f is large
and H has small index: if H has few elements, one may just use
`subHNF`

directly on the elements of H. For instance, let
K = ℚ(ζ_{p}, sqrt{m}) ⊂ L = ℚ(ζ_{f}), where p is
a prime, sqrt{m} is a quadratic number and f is the conductor of the
abelian extension K/ℚ. The following GP script creates H as the Galois
group of L/K, as a subgroub of (ℤ/fZ)^{*}:

HK(m, p, flag = 0)= { my(d = quaddisc(m), f = lcm(d, p), H); H = vectorsmall(f, a, a % p == 1 && kronecker(d,a) > 0); [f, znsubgroupgenerators(H,flag)]; } ? [f, v] = HK(36322, 5) time = 193 ms. %1 = [726440, [41, 61, 111, 131]] ? G = znstar(f,1); G.cyc; %2 = [1260, 12, 2, 2, 2, 2] ? A = subHNF(G, v) %3 = [2 0 1 1 0 1] [0 4 0 0 0 2] [0 0 1 0 0 0] [0 0 0 1 0 0] [0 0 0 0 1 0] [0 0 0 0 0 1] \\ Double check ? p = 5; d = quaddisc(36322); ? w = select(a->a % p == 1 && kronecker(d,a) > 0, [1..f]); #w time = 133 ms. %5 = 30240 \\ w enumerates the elements of H ? subHNF(G, w) == A \\ same result, about twice slower time = 242 ms. %6 = 1

This shows that K = ℚ(sqrt{36322},ζ_{5}) is contained in
ℚ(ζ_{726440}) and H = `<`

41, 61, 111, 131 `>`

.
Note that H = `<`

41`>`

`<`

61`>`

`<`

111 `>`

`<`

131`>`

is not a direct product. If *flag* = 1, then the function
finds generators which decompose H to direct factors:

? HK(36322, 5, 1) %3 = [726440, [41, 31261, 324611, 506221]]

This time
H = `<`

41`>`

x `<`

31261`>`

x
`<`

324611 `>`

x `<`

506221 `>`

.

The library syntax is `GEN `

.**znsubgroupgenerators**(GEN H, long flag)