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 **Z**/n**Z**). The following functions
compute discrete logarithms in the same groups (whenever this is meaningful)
`elllog`

, `fflog`

, `znlog`

.

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

***** `t_INT`

: the integer N,

***** `t_MAT`

: the factorization `fa = factor(N)`

,

***** `t_VEC`

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

.

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

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

Adds the integers contained in the
vector x (or the single integer x) to a special table of
"user-defined primes", and returns that table. Whenever `factor`

is
subsequently called, it will trial divide by the elements in this table.
If x is empty or omitted, just returns the current list of extra
primes.

The entries in x must be primes: there is no internal check, even if
the `factor_proven`

default is set. To remove primes from the list use
`removeprimes`

.

The library syntax is `GEN `

.**addprimes**(GEN x = NULL)

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

B).

***** If x is a `t_REAL`

or a `t_FRAC`

, this function uses continued
fractions.

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

By definition, a/b is the best rational approximation to x if
|b x - a| < |v x - u| for all integers (u,v) with 0 < v `<=`

B.
(Which implies that n/d is a convergent of the continued fraction of x.)

***** If x is a `t_INTMOD`

modulo N or a `t_PADIC`

of precision N =
p^k, this function performs rational modular reconstruction modulo N. The
routine then returns the unique rational number a/b in coprime integers
|a| < N/2B and b `<=`

B which is congruent to x modulo N. Omitting
B amounts to choosing it of the order of sqrt{N/2}. If rational
reconstruction is not possible (no suitable a/b exists), returns [].

? bestappr(Mod(18526731858, 11^10)) %1 = 1/7 ? bestappr(Mod(18526731858, 11^20)) %2 = [] ? bestappr(3 + 5 + 3*5^2 + 5^3 + 3*5^4 + 5^5 + 3*5^6 + O(5^7)) %2 = -1/3In most concrete uses, B is a prime power and we performed Hensel lifting to obtain x.

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

The library syntax is `GEN `

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

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

{degree}(b) `<=`

B).

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

or `t_SER`

, this function uses continued
fractions.

? bestapprPade((1-x^11)/(1-x)+O(x^11)) %1 = 1/(-x + 1) ? bestapprPade([1/(1+x+O(x^10)), (x^3-2)/(x^3+1)], 1) %2 = [1/(x + 1), -2]

***** If x is a `t_POLMOD`

modulo N or a `t_SER`

of precision N =
t^k, 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 which is congruent to x modulo
N. Omitting B amounts to choosing it of the order of N/2. If rational
reconstruction is not possible (no suitable a/b exists), returns [].

? bestapprPade(Mod(1+x+x^2+x^3+x^4, x^4-2)) %1 = (2*x - 1)/(x - 1) ? % * Mod(1,x^4-2) %2 = Mod(x^3 + x^2 + x + 3, x^4 - 2) ? bestapprPade(Mod(1+x+x^2+x^3+x^5, x^9)) %2 = [] ? bestapprPade(Mod(1+x+x^2+x^3+x^5, x^10)) %3 = (2*x^4 + x^3 - x - 1)/(-x^5 + x^3 + x^2 - 1)

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

The library syntax is `GEN `

.**bestapprPade**(GEN x, long B)

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)

binomial coefficient binom{x}{y}. Here y must be an integer, but x can be any PARI object.

The library syntax is `GEN `

.
The function
**binomial**(GEN x, long y)`GEN `

is also available, and so is
**binomialuu**(ulong n, ulong k)`GEN `

, which returns a vector v
with n+1 components such that v[k+1] = **vecbinome**(long n)`binomial`

(n,k) for k from
0 up to n.

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 library syntax is `GEN `

.**content**(GEN x)

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,20) %1 = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2] ? contfracpnqn(a,3) %2 = [3 22 333 355] [1 7 106 113] ? contfracpnqn(a,7) %3 = [3 22 333 355 103993 104348 208341 312689] [1 7 106 113 33102 33215 66317 99532]

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 non-zero,
returns the two-element row vector [d,f]. By convention, we write 0 = 0
x 1^2, so `core(0, 1)`

returns [0,1].

The library syntax is `GEN `

.
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 non-zero integer
n, this routine returns the (unique) fundamental discriminant d
such that n = df^2, f a positive rational number. If *flag* is non-zero,
returns the two-element row vector [d,f]. If n is congruent to
0 or 1 modulo 4, f is an integer, and a half-integer otherwise.

By convention, `coredisc(0, 1))`

returns [0,1].

Note that `quaddisc`

(n) returns the same value as `coredisc`

(n),
and also works with rational inputs n belongs to **Q**^*.

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 associated 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 present, output
only the first c coefficients in the series. The following command computes
the **sigma** function, associated to zeta(s)zeta(s-1):

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

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 sum_n x_n n^{-s} and sum_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 non-zero coefficient.

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

The library syntax is `GEN `

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

Creates a row vector whose components are the
divisors of x. The factorization of x (as output by `factor`

) can
be used instead.

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.

The library syntax is `GEN `

.**divisors**(GEN x)

Euler's phi (totient) function of the
integer |x|, in other words |(**Z**/x**Z**)^*|.

? eulerphi(40) %1 = 16

According to this definition we let phi(0) := 2, since **Z**^ *= {-1,1};
this is consistant with `znstar(0)`

: we have \kbd{znstar(n).no =
eulerphi(n)} for all n belongs to **Z**.

The library syntax is `GEN `

.**eulerphi**(GEN x)

General factorization function, where x is a rational (including integers), a complex number with rational real and imaginary parts, or a rational function (including polynomials). 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.

** Q and Q(i).**
See

`factorint`

for more information about the algorithms used.
The rational or Gaussian primes are in fact `ispseudoprime`

), a priori not rigorously proven primes. In fact,
any factor which is `<=`

10^{15} (whose norm is `<=`

10^{15} for an
irrational Gaussian prime) is a genuine prime. Use
`isprime`

to prove primality of other factors, as in

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

Another possibility is to set the global default `factor_proven`

, which
will perform a rigorous primality proof for each pseudoprime factor.

A `t_INT`

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

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

**Deprecated feature.** Setting *lim* = 0 is the same
as setting it to `primelimit`

+ 1. Don't use this: it is unwise to
rely on global variables when you can specify an explicit argument.

This routine uses trial division and perfect power tests, and should not be
used for huge values of *lim* (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 much more
efficient way.

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

**Rational functions.**
The polynomials or rational functions to be factored must have scalar
coefficients. In particular PARI does not know how to factor
*multivariate* polynomials. The following domains are currently
supported: **Q**, **R**, **C**, **Q**_p, finite fields and number fields.
See `factormod`

and `factorff`

for
the algorithms used over finite fields, `factornf`

for the algorithms
over number fields. Over **Q**, van Hoeij's method is used, which is
able to cope with hundreds of modular factors.

The routine guesses a sensible ring over which to factor: the
smallest ring containing all coefficients, taking into account quotient
structures induced by `t_INTMOD`

s and `t_POLMOD`

s (e.g. if a coefficient
in **Z**/n**Z** is known, all rational numbers encountered are first mapped to
**Z**/n**Z**; different moduli will produce an error). Factoring modulo a
non-prime number is not supported; to factor in **Q**_p, use `t_PADIC`

coefficients not `t_INTMOD`

modulo p^n.

? 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})In most cases, it is clearer and simpler to call an explicit variant than to rely on the generic

`factor`

function and
the above detection mechanism:

? factormod(T, 3) \\ over F_3 ? factorff(T, 3, t^2+t+2)) \\ 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

Note that factorization of polynomials is done up to multiplication by a constant. In particular, the factors of rational polynomials will have integer coefficients, and the content of a polynomial or rational function is discarded and not included in the factorization. If needed, you can always ask for the content explicitly:

? factor(t^2 + 5/2*t + 1) %1 = [2*t + 1 1] [t + 2 1] ? content(t^2 + 5/2*t + 1) %2 = 1/2

The irreducible factors are sorted by increasing degree.
See also `nffactor`

.

The library syntax is `GEN `

.
This function should only be used by the **gp_factor0**(GEN x, GEN lim = NULL)`gp`

interface. Use
directly `GEN `

or **factor**(GEN x)`GEN `

.
The obsolete function **boundfact**(GEN x, ulong lim)`GEN `

is kept for
backward compatibility.**factor0**(GEN x, long 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^3 * 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`

).

Factors the polynomial x modulo the
prime p, using distinct degree plus
Cantor-Zassenhaus. The coefficients of x must be
operation-compatible with **Z**/p**Z**. The result is a two-column matrix, the
first column being the irreducible polynomials dividing x, and the second
the exponents. If you want only the *degrees* of the irreducible
polynomials (for example for computing an L-function), use
`factormod`

(x,p,1). Note that the `factormod`

algorithm is
usually faster than `factorcantor`

.

The library syntax is `GEN `

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

Factors the polynomial x in the field
**F**_q defined by the irreducible polynomial a over **F**_p. The
coefficients of x must be operation-compatible with **Z**/p**Z**. The result
is a two-column matrix: the first column contains the irreducible factors of
x, and the second their exponents. If all the coefficients of x are in
**F**_p, a much faster algorithm is applied, using the computation of
isomorphisms between finite fields.

Either a or p can omitted (in which case both are ignored) if x has
`t_FFELT`

coefficients; the function then becomes identical to `factor`

:

? factorff(x^2 + 1, 5, y^2+3) \\ over F_5[y]/(y^2+3) ~ F_25 %1 = [Mod(Mod(1, 5), Mod(1, 5)*y^2 + Mod(3, 5))*x + Mod(Mod(2, 5), Mod(1, 5)*y^2 + Mod(3, 5)) 1] [Mod(Mod(1, 5), Mod(1, 5)*y^2 + Mod(3, 5))*x + Mod(Mod(3, 5), Mod(1, 5)*y^2 + Mod(3, 5)) 1] ? t = ffgen(y^2 + Mod(3,5), 't); \\ a generator for F_25 as a t_FFELT ? factorff(x^2 + 1) \\ not enough information to determine the base field *** at top-level: factorff(x^2+1) *** ^--------------- *** factorff: incorrect type in factorff. ? factorff(x^2 + t^0) \\ make sure a coeff. is a t_FFELT %3 = [x + 2 1] [x + 3 1] ? factorff(x^2 + t + 1) %11 = [x + (2*t + 1) 1] [x + (3*t + 4) 1]

Notice that the second syntax is easier to use and much more readable.

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 is 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 x modulo the prime integer p, using
Berlekamp. The coefficients of x must be operation-compatible with
**Z**/p**Z**. The result is a two-column matrix, the first column being the
irreducible polynomials dividing x, and the second the exponents. If *flag*
is non-zero, outputs only the *degrees* of the irreducible polynomials
(for example, for computing an L-function). A different algorithm for
computing the mod p factorization is `factorcantor`

which is sometimes
faster.

The library syntax is `GEN `

.**factormod0**(GEN x, GEN p, long flag)

Return a `t_FFELT`

generator for the finite field with q elements;
q = p^f must be a prime power. This functions computes an irreducible
monic polynomial P belongs to **F**_p[X] of degree f (via `ffinit`

) and
returns g = X (mod P(X)). If `v`

is given, the variable name is used
to display g, else the variable x is used.

? g = ffgen(8, 't); ? g.mod %2 = t^3 + t^2 + 1 ? g.p %3 = 2 ? g.f %4 = 3 ? ffgen(6) *** at top-level: ffgen(6) *** ^-------- *** ffgen: not a prime number in ffgen: 6.Alternative syntax: instead of a prime power q, one may input directly the polynomial P (monic, irreducible, with

`t_INTMOD`

coefficients), and the function returns the generator g = X (mod P(X)),
inferring p from the coefficients of P. If `v`

is given, the
variable name is used to display g, else the variable of the polynomial
P is used. If P is not irreducible, we create an invalid object and
behaviour of functions dealing with the resulting `t_FFELT`

is undefined; in fact, it is much more costly to test P for
irreducibility than it would be to produce it via `ffinit`

.
The library syntax is `GEN `

, where **ffgen**(GEN q, 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)`1+ffgen(x*Mod(1,p),v)`

.

Computes a monic polynomial of degree n which is irreducible over
**F**_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 **F**_{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.

Discrete logarithm of the finite field element x in base g, i.e.
an e in **Z** such that g^e = o. If
present, o represents the multiplicative order of g, see
Section [**Label: 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`

.

If no o is given, assume that g is a primitive root. 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)

Computes the number of monic irreducible polynomials over **F**_q of degree exactly n,
(*flag* = 0 or omited) 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 [**Label: 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 factorisation `factor(ord)`

. This last parameter is
useful in `fflog`

and `fforder`

, see Section [**Label: 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)

x-th Fibonacci number.

The library syntax is `GEN `

.**fibo**(long x)

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

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) %2 = 4A good quantitative check to decide whether such a gcd "should be" non-trivial, is to use

`polresultant`

: a value
close to 0 means that a small deformation of the inputs has non-trivial gcd.
You may also use `bezout`

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

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 x is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise.

The library syntax is `long `

.**isfundamental**(GEN x)

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.

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 ^ poldegree(x.pol) - 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 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 prime and has more than 1000 digits, say.
Use `ispseudoprime`

to quickly check for compositeness. See also
`factor`

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

If *flag* = 0, use a combination of Baillie-PSW pseudo primality test (see
`ispseudoprime`

), Selfridge "p-1" test if x-1 is smooth enough, and
Adleman-Pomerance-Rumely-Cohen-Lenstra (APRCL) for general x.

If *flag* = 1, use Selfridge-Pocklington-Lehmer "p-1" test and output a
primality certificate as follows: return

***** 0 if x is composite,

***** 1 if x is small enough that passing Baillie-PSW test guarantees
its primality (currently x < 2^{64}, as checked by Jan Feitsma),

***** 2 if x is a large prime whose primality could only sensibly be
proven (given the algorithms implemented in PARI) using the APRCL test.

***** Otherwise (x is large and x-1 is smooth) output a three column
matrix as a primality certificate. The first column contains prime
divisors p of x-1 (such that prod p^{v_p(x-1)} > x^{1/3}), the second
the corresponding elements a_p as in Proposition 8.3.1 in GTM 138 , and the
third the output of isprime(p,1).

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

at a high enough debug level, you may see warnings about
untested integers being declared primes. This is normal: we ask for partial
factorisations (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 `isprime`

procedure.

If *flag* = 2, use APRCL.

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 is a Baillie-Pomerance-Selfridge-Wagstaff
pseudo prime (strong Rabin-Miller pseudo prime for base 2, followed by
strong Lucas test for the sequence (P,-1), P smallest positive integer
such that P^2 - 4 is not a square mod x).

There are no known composite numbers passing this 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)

True (1) if x is a square, false (0)
if not. What "being a square" means depends on the type of x: all
`t_COMPLEX`

are squares, as well as all non-negative `t_REAL`

; for
exact types such as `t_INT`

, `t_FRAC`

and `t_INTMOD`

, squares are
numbers of the form s^2 with s in **Z**, **Q** and **Z**/N**Z** 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.

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

and `gen_1`

instead of a boolean value.

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

The library syntax is `long `

.**issquarefree**(GEN x)

True (1) if x = phi(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 **Z** x **Z**
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) = {abs}(x*y). If y is omitted and x is a vector, returns the {lcm} of all components of x.

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^4, i, random); ? lcm(v); time = 323 ms. ? l = v[1]; for (i = 1, #v, l = lcm(l, v[i])) time = 833 ms.

The library syntax is `GEN `

.**glcm0**(GEN x, GEN y = NULL)

Return the largest integer e so that b^e `<=`

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

? logint(1000, 2) %1 = 9 ? 2^9 %2 = 512 ? logint(1000, 2, &z) %3 = 9 ? z %4 = 512The number of digits used to write b in base x is

`1 + logint(x,b)`

:

? #digits(1000!, 10) %5 = 2568 ? logint(1000!, 10) %6 = 2567This 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 mu-function of |x|. x must be of type 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`

.

The library syntax is `GEN `

.**nextprime**(GEN x)

Gives the number of unrestricted partitions of
n, usually called p(n) in the literature; in other words the number of
nonnegative integer solutions to a+2b+3c+.. .= n. n must be of type
integer and n < 10^{15} (with trivial values p(n) = 0 for n < 0 and
p(0) = 1). The algorithm uses the Hardy-Ramanujan-Rademacher formula.
To explicitly enumerate them, see `partitions`

.

The library syntax is `GEN `

.**numbpart**(GEN n)

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)

Returns the vector of partitions of the integer k as a sum of positive
integers (parts); for k < 0, it returns the empty set `[]`

, and for k
= 0 the trivial partition (no parts). A partition is given by a
`t_VECSMALL`

, where parts are sorted in nondecreasing order:

? partitions(3) %1 = [Vecsmall([3]), Vecsmall([1, 2]), Vecsmall([1, 1, 1])]correspond to 3, 1+2 and 1+1+1. The number of (unrestricted) partitions of k is given by

`numbpart`

:

? #partitions(50) %1 = 204226 ? numbpart(50) %2 = 204226

Optional parameters n and a are as follows:

***** n = *nmax* (resp. n = [*nmin*,*nmax*]) restricts
partitions to length less than *nmax* (resp. length between
*nmin* and nmax), where the *length* is the number of nonzero
entries.

***** a = *amax* (resp. a = [*amin*,*amax*]) restricts the parts
to integers less than *amax* (resp. between *amin* and
*amax*).

? partitions(4, 2) \\ parts bounded by 2 %1 = [Vecsmall([2, 2]), Vecsmall([1, 1, 2]), Vecsmall([1, 1, 1, 1])] ? partitions(4,, 2) \\ at most 2 parts %2 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])] ? partitions(4,[0,3], 2) \\ at most 2 parts %3 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])]

By default, parts are positive and we remove zero entries unless
amin `<=`

0, in which case nmin is ignored and X is of constant length
*nmax*:

? partitions(4, [0,3]) \\ parts between 0 and 3 %1 = [Vecsmall([0, 0, 1, 3]), Vecsmall([0, 0, 2, 2]),\ Vecsmall([0, 1, 1, 2]), Vecsmall([1, 1, 1, 1])]

The library syntax is `GEN `

.**partitions**(long k, GEN a = NULL, GEN n) = NULL)

Returns the vector of distinct roots of the polynomial x in the field
**F**_q defined by the irreducible polynomial a over **F**_p. The
coefficients of x must be operation-compatible with **Z**/p**Z**.
Either a or p can omitted (in which case both are ignored) if x has
`t_FFELT`

coefficients:

? polrootsff(x^2 + 1, 5, y^2+3) \\ over F_5[y]/(y^2+3) ~ F_25 %1 = [Mod(Mod(3, 5), Mod(1, 5)*y^2 + Mod(3, 5)), Mod(Mod(2, 5), Mod(1, 5)*y^2 + Mod(3, 5))] ? t = ffgen(y^2 + Mod(3,5), 't); \\ a generator for F_25 as a t_FFELT ? polrootsff(x^2 + 1) \\ not enough information to determine the base field *** at top-level: polrootsff(x^2+1) *** ^----------------- *** polrootsff: incorrect type in factorff. ? polrootsff(x^2 + t^0) \\ make sure one coeff. is a t_FFELT %3 = [3, 2] ? polrootsff(x^2 + t + 1) %4 = [2*t + 1, 3*t + 4]

Notice that the second syntax is easier to use and much more readable.

The library syntax is `GEN `

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

Finds the largest pseudoprime (see
`ispseudoprime`

) less than or equal to x. x can be of any real type.
Returns 0 if x `<=`

1. Note that if x is a prime, this function returns x
and not the largest prime strictly smaller than x. To rigorously prove that
the result is prime, use `isprime`

.

The library syntax is `GEN `

.**precprime**(GEN x)

The n-th prime number

? prime(10^9) %1 = 22801763489Uses checkpointing and a naive O(n) algorithm.

The library syntax is `GEN `

.**prime**(long n)

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 = 4118054813Uses checkpointing and a naive O(x) algorithm.

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. In the present version **2.7.0**, a O(D^{1/2})
algorithm is used for D > 0 (using Euler product and the functional
equation) so D should not be too large, say D < 10^8, for the time to be
reasonable. On the other hand, for D < 0 one can reasonably compute
`qfbclassno(D)`

for |D| < 10^{25}, since the routine uses
Shanks's method which is in O(|D|^{1/4}). For larger values of |D|,
see `quadclassunit`

.

If *flag* = 1, compute the class number using Euler products and the
functional equation. However, it is in O(|D|^{1/2}).

**Important warning.** For D < 0, this function may give incorrect
results when the class group has many cyclic factors,
because implementing Shanks's method in full generality slows it down
immensely. It is therefore strongly recommended to double-check results using
either the version with *flag* = 1 or the function `quadclassunit`

.

**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 varepsilon has negative norm; when D
> 0 and Nvarepsilon > 0, the number of classes of forms is twice the
ordinary class number. This is a problem which we cannot fix for backward
compatibility reasons. Use the following routine if you are only interested
in the number of classes of forms:

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

Here are a few examples:

? qfbclassno(400000028) time = 3,140 ms. %1 = 1 ? quadclassunit(400000028).no time = 20 ms. \\ { much faster} %2 = 1 ? qfbclassno(-400000028) time = 0 ms. %3 = 7253 \\ { correct, and fast enough} ? quadclassunit(-400000028).no time = 0 ms. %4 = 7253

See also `qfbhclassno`

.

The library syntax is `GEN `

.
The following functions are also available:**qfbclassno0**(GEN D, long flag)

`GEN `

(**classno**(GEN D)*flag* = 0)

`GEN `

(**classno2**(GEN D)*flag* = 1).

Finally

`GEN `

computes the class number of an imaginary
quadratic field by counting reduced forms, an O(|D|) algorithm.**hclassno**(GEN D)

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)

Hurwitz class number of x, where
x is non-negative and congruent to 0 or 3 modulo 4. For x > 5.
10^5, we assume the GRH, and use `quadclassunit`

with default
parameters.

The library syntax is `GEN `

.**hclassno**(GEN x)

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

) using the NUCOMP
and NUDUPL algorithms of Shanks, à la Atkin. L is any positive
constant, but for optimal speed, one should take L = |D|^{1/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 straightforward and in general *slower*
than the generic routine (since the latter takes advantage of asymptotically
fast operations and careful optimizations).

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
(see `qfbnucomp`

, in particular the final warning).

The library syntax is `GEN `

.**nupow**(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 non-negative 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_QFI`

are not implemented.) In the case where x > 0,
the "distance" component of the form is set equal to zero according to the
current precision.

The library syntax is `GEN `

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

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

1: perform a single reduction step

2: don't update Shanks's distance

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

The library syntax is `GEN `

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

`GEN `

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

and for indefinite forms:

`GEN `

**redreal**(GEN x)

`GEN `

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

),

`GEN `

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

),

`GEN `

( = **rhorealnod**(GEN x, GEN isd)`qfbred(x,3,,isd)`

).

Solve the equation Q(x,y) = p over the integers, where Q is a binary quadratic form and p a prime number.

Return [x,y] as a two-components vector, or zero if there is no solution. Note that this function returns only one solution and not all the solutions.

Let D = disc Q. The algorithm used runs in probabilistic polynomial time
in p (through the computation of a square root of D modulo p); it is
polynomial time in D if Q is imaginary, but exponential time if Q is
real (through the computation of a full cycle of reduced forms). In the
latter case, note that `bnfisprincipal`

provides a solution in heuristic
subexponential time in D assuming the GRH.

The library syntax is `GEN `

.**qfbsolve**(GEN Q, GEN p)

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

This function should be used instead of `qfbclassno`

or `quadregula`

when D < -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 non-negative real numbers which control the execution
time and the stack size, see [**Label: 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 famous result of Bach
states that c_2 = 6 is fine, 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 is
usually more convenient to use
`GEN `

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

Discriminant of the quadratic field **Q**(sqrt{x}), where x belongs to **Q**.

The library syntax is `GEN `

.**quaddisc**(GEN x)

Creates the quadratic number omega = (a+sqrt{D})/2 where a = 0 if D = 0 mod 4, a = 1 if D = 1 mod 4, so that (1,omega) 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.

The library syntax is `GEN `

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

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 sigma 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 field of positive discriminant x. Returns
an error if x is not a discriminant (fundamental or not) or if x is a
square. See also `quadclassunit`

if x is large.

The library syntax is `GEN `

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

Fundamental unit of the
real quadratic field **Q**(sqrt D) where D is the positive discriminant
of the field. If D is not a fundamental discriminant, this probably gives
the fundamental unit of the corresponding order. D must be an integer
congruent to 0 or 1 modulo 4, which is not a square; the result is a
quadratic number (see Section [**Label: se:quadgen**]).

The library syntax is `GEN `

.**quadunit**(GEN D)

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.

The library syntax is `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 non-negative integer.

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

The library syntax is `GEN `

.**sqrtint**(GEN x)

Returns the integer n-th root of x, i.e. the largest integer y such
that y^n `<=`

x, where x is a non-negative integer.

? N = 120938191237; sqrtnint(N, 5) %1 = 164 ? N^(1/5) %2 = 164.63140849829660842958614676939677391The special case n = 2 is

`sqrtint`

The library syntax is `GEN `

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

Stirling number of the first kind s(n,k) (*flag* = 1, default) or
of the second kind S(n,k) (*flag* = 2), where n, k are non-negative
integers. The former is (-1)^{n-k} times the
number of permutations of n symbols with exactly k cycles; the latter is
the number of ways of partitioning a set of n elements into k non-empty
subsets. Note that if all s(n,k) are needed, it is much faster to compute
sum_k s(n,k) x^k = x(x-1)...(x-n+1).
Similarly, if a large number of S(n,k) are needed for the same k,
one should use
sum_n S(n,k) x^n = (x^k)/((1-x)...(1-kx)).
(Should be implemented using a divide and conquer product.) Here are
simple variants for n fixed:

/* list of s(n,k), k = 1..n */ vecstirling(n) = Vec( factorback(vector(n-1,i,1-i*'x)) ) /* list of S(n,k), k = 1..n */ vecstirling2(n) = { my(Q = x^(n-1), t); vector(n, i, t = divrem(Q, x-i); Q=t[1]; t[2]); }

The library syntax is `GEN `

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

(**stirling1**(ulong n, ulong k)*flag* = 1) and `GEN `

(**stirling2**(ulong n, ulong k)*flag* = 2).

Returns the Dedekind sum associated 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 (decimal) digits in the integer n.

? sumdigits(123456789) %1 = 45Other bases that 10 are not supported. Note that the sum of bits in n is returned by

`hammingweight`

.
The library syntax is `GEN `

.**sumdigits**(GEN n)

N being an integer and P belongs to **Z**[X], finds all integers x with
|x| `<=`

X such that
gcd(N, P(x)) `>=`

B,
using Coppersmith's algorithm (a famous application of the LLL
algorithm). X must be smaller than exp(log^2 B / (deg(P) log N)):
for B = N, this means X < N^{1/deg(P)}. Some x larger than X may
be returned if you are very lucky. The smaller B (or the larger X), the
slower the routine will be. The strength of Coppersmith method is the
ability to find roots modulo a general *composite* N: if N is a prime
or a prime power, `polrootsmod`

or `polrootspadic`

will be much
faster.

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

setrand(1); \\ to make the example reproducible p = nextprime(random(10^30)); q = nextprime(random(10^30)); N = p*q; p0 = p % 10^20; \\ assume we know 1) p > 10^29, 2) the last 19 digits of p p1 = zncoppersmith(10^19*x + p0, N, 10^12, 10^29) \\ result in 10ms. %1 = [35023733690] ? gcd(p1[1] * 10^19 + p0, N) == p %2 = 1and we recovered p, faster than by trying all possibilities < 10^{12}.

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. %3 = [265174753892462432] ? %[1] == x0 %4 = 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)

Discrete logarithm of x in (**Z**/N**Z**)^* in base g.
The result is [] when x is not a power of g.
If present, o represents the multiplicative order of g, see
Section [**Label: 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 = 205 ms. ? for(i=1,10^4, znlog(i, g)) time = 244 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 (**Z**/N**Z**)^* 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 `

.**znlog**(GEN x, GEN g, GEN o = NULL)

x must be an integer mod n, and the
result is the order of x in the multiplicative group (**Z**/n**Z**)^*. Returns
an error if x is not invertible.
The parameter o, if present, represents a non-zero
multiple of the order of x, see Section [**Label: 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 `

.
Also available is **znorder**(GEN x, GEN o = NULL)`GEN `

.**order**(GEN x)

Returns a primitive root (generator) of (**Z**/n**Z**)^*, whenever this
latter group is cyclic (n = 4 or n = 2p^k or n = p^k, where p is an
odd prime and k `>=`

0). If the group is not cyclic, the result is
undefined. If n is a prime power, then the smallest positive primitive
root is returned. This may not be true for n = 2p^k, p odd.

Note that this function requires factoring p-1 for p as above,
in order to determine the exact order of elements in
(**Z**/n**Z**)^*: 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
(**Z**/n**Z**)^* as a 3-component row vector v, where v[1] = phi(n) is the
order of that group, v[2] is a k-component row-vector d of integers
d[i] such that d[i] > 1 and d[i] | d[i-1] for i `>=`

2 and
(**Z**/n**Z**)^* ~ prod_{i = 1}^k(**Z**/d[i]**Z**), and v[3] is a k-component row
vector giving generators of the image of the cyclic groups **Z**/d[i]**Z**.

? 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]According to the above definitions,

`znstar(0)`

is
`[2, [2], [-1]]`

, corresponding to
The library syntax is `GEN `

.**znstar**(GEN n)