### Arithmetic functions

#### Arithmetic functions and factorization

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

  default(factor_proven, 1)

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

#### Orders in finite groups and Discrete Logarithm functions

The following functions compute the order of an element in a finite group: ellorder (the rational points on an elliptic curve defined over a finite field), fforder (the multiplicative group of a finite field), znorder (the invertible elements in Z/nZ). 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.) The meaning of that optional argument is as follows and depends on its type

* 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 \kbd{[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


#### addprimes({x = []})

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

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

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

#### bestappr(x, {A},{B})

if B is omitted, finds the best rational approximation to x belongs to R using continued fractions. If A is omitted, return the best approximation affordable given the input accuracy; otherwise make sure that denominator is at most equal to A.

If B is present perform rational modular reconstruction (see below). In both cases, the function applies recursively to components of complex objects (polynomials, vectors,...).

? 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, n/d is the best rational approximation to x if |d x - n| < |v x - u| for all integers (u,v) with v <= A. (Which implies that n/d is a convergent of the continued fraction of x.)

If x is an t_INTMOD, (or a recursive combination of those), modulo N say, B must be present. The routine then returns the unique rational number a/b in coprime integers a <= A and b <= B which is congruent to x modulo N. If N <= 2AB, uniqueness is not guaranteed and the function fails with an error message. If rational reconstruction is not possible (no such a/b exists for at least one component of x), returns -1.

? bestappr(Mod(18526731858, 11^10), 10^10, 10^10)
***   at top-level: bestappr(Mod(1852673
***                 ^--------------------
*** bestappr: ratlift: must have 2*amax*bmax < m, found
amax=10000000000
bmax=10000000000
m=25937424601
? bestappr(Mod(18526731858, 11^10), 10^5, 10^5)
%1 = 1/7
? bestappr(Mod(18526731858, 11^20), 10^10, 10^10)
%2 = -1

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

If x is a t_POLMOD, modulo T say, B must be present. The routine then returns the unique rational function P/Q with deg P <= A and deg Q <= B which is congruent to x modulo T. If deg T <= A+B, uniqueness is not guaranteed and the function fails with an error message. If rational reconstruction is not possible, returns -1.

The library syntax is GEN bestappr0(GEN x, GEN A = NULL, GEN B = NULL). Also available is GEN bestappr(GEN x, GEN A).

#### bezout(x,y)

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

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

? bezout(a,a)
%2 = [0, 1, x + 0.E-28]

? bezout(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 vecbezout(GEN x, GEN y).

#### bezoutres(x,y)

finds u and v such that x*u + y*v = d, where d is the resultant of x and y. The result is the row vector [u,v,d]. The algorithm used (subresultant) assumes that the base ring is a domain.

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

#### bigomega(x)

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

? factor(392)
%1 =
[2 3]

[7 2]

? bigomega(392)
%2 = 5;  \\ = 3+2
? omega(392)
%3 = 2;  \\ without multiplicity


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

The library syntax is GEN gbigomega(GEN x). For a t_INT x, the variant long bigomega(GEN n) is generally easier to use.

#### binomial(x,y)

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

The library syntax is GEN binomial(GEN x, long y). The function GEN binomialuu(ulong n, ulong k) is also available, and so is GEN vecbinome(long n), which returns a vector v with n+1 components such that v[k+1] = binomial(n,k) for k from 0 up to n.

#### chinese(x,{y})

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

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

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

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

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

#### content(x)

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

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

returns the row vector whose components are the partial quotients of the continued fraction expansion of x. In other words, a result [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 contfrac0(GEN x, GEN b = NULL, long nmax). Also available are GEN gboundcf(GEN x, long nmax), GEN gcf(GEN x) and GEN gcf2(GEN b, GEN x).

#### contfracpnqn(x)

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.

The library syntax is GEN pnqn(GEN x).

#### core(n,{flag = 0})

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

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

#### coredisc(n,{flag = 0})

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

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

Note that quaddisc(n) returns the same value as coredisc(n), and also works with rational inputs n belongs to Q^*.

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

#### dirdiv(x,y)

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

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

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

computes the Dirichlet series associated to the \idx{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)

#### dirmul(x,y)

x and y being vectors of perhaps different lengths considered as Dirichlet series, computes the product of x by y, again as a vector.

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

#### divisors(x)

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

#### eulerphi(x)

Euler's phi (totient) function of |x|, in other words |(Z/xZ)^*|. Normally, x must be of type integer, but the function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN geulerphi(GEN x). For a t_INT x, the variant GEN eulerphi(GEN n) is also available.

#### factor(x,{lim})

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 pseudoprimes (see 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 and satisfy lim <= primelimit + 1; setting lim = 0 is the same as setting it to primelimit + 1. 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]


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]

? default(primelimit,10^9)
time = 4,360 ms.
%5 = 1000000000
? factor(F, 10^9)    \\ very slow
time = 6,120 ms.
%6 =
[1009 1]

[100003 1]

[340282366920938463463374607431768211457 1]

? factorint(F, 1+8)  \\ much faster, all small primes were found
time = 40 ms.
%7 =
[1009 1]

[100003 1]

[340282366920938463463374607431768211457 1]

? factorint(F)   \\ complete factorisation
time = 260 ms.
%8 =
[1009 1]

[100003 1]

[59649589127497217 1]

[5704689200685129054721 1]


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. 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 you want to factor: the smallest ring containing all coefficients, taking into account quotient structures induced by t_INTMODs and t_POLMODs (e.g.if a coefficient in Z/nZ is known, all rational numbers encountered are first mapped to Z/nZ; different moduli will produce an error). 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


See also nffactor.

The library syntax is GEN gp_factor0(GEN x, GEN lim = NULL). This function should only be used by the gp interface. Use directly GEN factor(GEN x) or GEN boundfact(GEN x, long lim). The obsolete function GEN factor0(GEN x, long lim) is kept for backward compatibility.

#### factorback(f,{e})

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

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

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

? factor(12)
%1 =
[2 2]

[3 1]

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


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

#### factorcantor(x,p)

factors the polynomial x modulo the prime p, using distinct degree plus Cantor-Zassenhaus. The coefficients of x must be operation-compatible with Z/pZ. 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).

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

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/pZ. 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(x)

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

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

#### factorint(x,{flag = 0})

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

By convention 0 is factored as 0^1, and 1 as the empty factorization; also the divisors are by default not proven primes 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).

#### factormod(x,p,{flag = 0})

factors the polynomial x modulo the prime integer p, using Berlekamp. The coefficients of x must be operation-compatible with Z/pZ. 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).

#### ffgen(P,{v})

return the generator g = X (mod P(X)) of the finite field defined by the polynomial P (which must have t_INTMOD coefficients). If v is given, the variable name is used to display g, else the variable of the polynomial P is used.

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

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

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-Lenstra's algorithm.

It is useful in conjunction with ffgen; for instance if \kbd{P = ffinit(3,2)}, you can represent elements in F_{3^2} in term of \kbd{g = ffgen(P,g)}.

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

#### fflog(x,g,{o})

discrete logarithm of the finite field element x in base 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. It may be set as a side effect of calling ffprimroot.

If no o is given, assume that g is a primitive root. See znlog for the limitations of the underlying discrete log algorithms.

? t = ffgen(ffinit(7,5));
? o = fforder(t)
%2 = 5602   \\ not a primitive root.
? fflog(t^10,t)
%3 = 11214  \\ Actually correct modulo o. We are lucky !
? 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  \\ no surprise there !
? fforder(g^10000, g, o)
? fflog(g^10000, g, o)
%9 = 10000


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

#### fforder(x,{o})

multiplicative order of the finite field element x. If o is present, it represents a multiple of the order of the element, see Section [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).

#### ffprimroot(x, {&o})

return a primitive root of the multiplicative group of the definition field of the finite field element x (not necessarily the same as the field generated by x). If present, o is set to a vector [ord, fa], where ord is the order of the group and fa its 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).

#### fibonacci(x)

x-th Fibonacci number.

The library syntax is GEN fibo(long x).

#### gcd(x,{y})

creates the greatest common divisor of x and y. If you also need the u and v such that x*u + y*v = gcd(x,y), use the 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 = 4

A good quantitative check to decide whether such a gcd "should be" non-trivial, is to use polresultant: a value close to 0 means that a small deformation of the inputs has non-trivial gcd. You may also use 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 ggcd0(GEN x, GEN y = NULL). Also available are GEN ggcd(GEN x, GEN y), if y is not NULL, and GEN content(GEN x), if y = NULL.

#### hilbert(x,y,{p})

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

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

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

#### isfundamental(x)

true (1) if x is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gisfundamental(GEN x).

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

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 ispower(GEN x, GEN k = NULL, GEN *n = NULL). Also available is long gisanypower(GEN x, GEN *pty) (k omitted).

#### isprime(x,{flag = 0})

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

This routine proves or disproves rigorously that a number is prime, which can be very slow when x is indeed 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 Proposition8.3.1 in GTM138 , 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).

#### ispseudoprime(x,{flag})

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

If flag = 0, checks whether x 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).

#### issquare(x,{&n})

true (1) if x is a square, false (0) if not. What "being a square" means depends on the type of x: all t_COMPLEX are squares, as well as all non-negative t_REAL; for exact types such as t_INT, t_FRAC and t_INTMOD, squares are numbers of the form s^2 with s in Z, Q and Z/NZ 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
? issquare([4, x^2], &n)
%3 = [1, 1]  \\ both are squares
? n
%4 = [2, x]  \\ the square roots


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 function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gissquareall(GEN x, GEN *n = NULL). Also available is GEN gissquare(GEN x).

#### issquarefree(x)

true (1) if x is squarefree, false (0) if not. Here x can be an integer or a polynomial. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gissquarefree(GEN x). For scalar arguments x ( t_INT or t_POL), the function long issquarefree(GEN x) is easier to use.

#### kronecker(x,y)

Kronecker symbol (x|y), where x and y must be of type integer. By definition, this is the extension of Legendre symbol to 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 GEN gkronecker(GEN x, GEN y).

#### lcm(x,{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).

#### moebius(x)

Moebius mu-function of |x|. x must be of type integer. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gmoebius(GEN x). For a t_INT x, the variant long moebius(GEN n) is generally easier to use.

#### nextprime(x)

finds the smallest pseudoprime (see ispseudoprime) greater than or equal to x. x can be of any real type. Note that if x is a pseudoprime, this function returns x and not the smallest pseudoprime strictly larger than x. To rigorously prove that the result is prime, use isprime. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gnextprime(GEN x). For a scalar x, long nextprime(GEN n) is also available.

#### numbpart(n)

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

#### numdiv(x)

number of divisors of |x|. x must be of type integer. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gnumbdiv(GEN x). If x is a t_INT, one may use GEN numbdiv(GEN n) directly.

#### omega(x)

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

? factor(392)
%1 =
[2 3]

[7 2]

? omega(392)
%2 = 2;  \\ without multiplicity
? bigomega(392)
%3 = 5;  \\ = 3+2, with multiplicity


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

The library syntax is GEN gomega(GEN x). For a t_INT x, the variant long omega(GEN n) is generally easier to use.

#### partitions(n,{restr = 0})

returns vector of partitions of the integer n (negative values return [], n = 0 returns the trivial partition of the empty set). The second optional argument may be set to a non-negative number smaller than n to restrict the value of each element in the partitions to that value. The default of 0 means that this maximum is n itself.

A partition is given by a t_VECSMALL:

? partitions(4, 2)
%1 = [Vecsmall([2, 2]), Vecsmall([1, 1, 2]), Vecsmall([1, 1, 1, 1])]


correspond to 2+2, 1+1+2, 1+1+1+1.

The library syntax is GEN partitions(long n, long restr).

#### polrootsff(x,{p},{a})

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

#### precprime(x)

finds the largest pseudoprime (see ispseudoprime) less than or equal to x. x can be of any real type. Returns 0 if x <= 1. Note that if x is a prime, this function returns x and not the largest prime strictly smaller than x. To rigorously prove that the result is prime, use isprime. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN gprecprime(GEN x). For a scalar x, long precprime(GEN n) is also available.

#### prime(n)

the x-th prime number, which must be among the precalculated primes.

The library syntax is GEN prime(long n).

#### primepi(x)

the prime counting function. Returns the number of primes p, p <= x. Uses a naive algorithm so that x must be less than primelimit.

The library syntax is GEN primepi(GEN x).

#### primes(x)

creates a row vector whose components are the first x prime numbers, which must be among the precalculated primes.

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


The library syntax is GEN primes(long x).

#### qfbclassno(D,{flag = 0})

ordinary class number of the quadratic order of discriminant D. In the present version 2.5.1, 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 qfbclassno0(GEN D, long flag). The following functions are also available:

GEN classno(GEN D) (flag = 0)

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

Finally

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

#### qfbcompraw(x,y)

composition of the binary quadratic forms x and y, without reduction of the result. This is useful e.g.to compute a generating element of an ideal.

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

#### qfbhclassno(x)

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

#### qfbnucomp(x,y,L)

composition of the primitive positive definite binary quadratic forms x and y (type t_QFI) using the NUCOMP and NUDUPL algorithms of Shanks, à la Atkin. L is any positive constant, but for optimal speed, one should take L = |D|^{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 nucomp(GEN x, GEN y, GEN L). Also available is GEN nudupl(GEN x, GEN L) when x = y.

#### qfbnupow(x,n)

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

#### qfbpowraw(x,n)

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

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

#### qfbprimeform(x,p)

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

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

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

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

1: perform a single reduction step

2: don't update Shanks's distance

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

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

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

and for indefinite forms:

GEN redreal(GEN x)

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

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

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

#### qfbsolve(Q,p)

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

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

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

This function should be used instead of qfbclassno or quadregula when D < -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 positive 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 = 0.2 means that relations are easier to find, but more relations are needed and the linear algebra will be harder. 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 \leq 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 quadclassunit0(GEN D, long flag, GEN tech = NULL, long prec). If you really need to experiment with the tech parameter, it is usually more convenient to use GEN Buchquad(GEN D, double c1, double c2, long prec)

#### quaddisc(x)

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

The library syntax is GEN quaddisc(GEN x).

#### quadgen(D)

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

#### quadhilbert(D)

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

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

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

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

#### quadpoly(D,{v = x})

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

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

#### quadray(D,f)

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

For D < 0, uses the 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).

#### quadregulator(x)

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

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

#### quadunit(D)

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

#### removeprimes({x = []})

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

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

#### sigma(x,{k = 1})

sum of the k-th powers of the positive divisors of |x|. x and k must be of type integer. The function accepts vector/matrices arguments for x, and is then applied componentwise.

The library syntax is GEN gsumdivk(GEN x, long k). Also available are GEN gsumdiv(GEN n) (k = 1), GEN sumdivk(GEN n,long k) (n a t_INT) and GEN sumdiv(GEN n) (k = 1, n a t_INT)

#### sqrtint(x)

integer square root of x, which must be a non-negative integer. The result is non-negative and rounded towards zero.

The library syntax is GEN sqrtint(GEN x).

#### stirling(n,k,{flag = 1})

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 stirling(long n, long k, long flag). Also available are GEN stirling1(ulong n, ulong k) (flag = 1) and GEN stirling2(ulong n, ulong k) (flag = 2).

#### sumdedekind(h,k)

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

#### zncoppersmith(P, N, X, {B = 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 = 1

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


We guessed an integer of the order of 10^{18} in a couple of seconds.

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

#### znlog(x,g,{o})

discrete logarithm of x in (Z/NZ)^* in base 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. If no o is given, assume that g generate (Z/NZ)^*.

This function uses a simple-minded combination of generic discrete log algorithms (index calculus methods are not yet implemented).

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

* Pollard rho method (q large).

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.

? 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

The result is undefined when x is not a power of g or when x is not invertible mod N:

? znlog(6, Mod(2,3))
***   at top-level: znlog(6,Mod(2,3))
***                 ^-----------------
*** znlog: impossible inverse modulo: Mod(0, 3).

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

#### znorder(x,{o})

x must be an integer mod n, and the result is the order of x in the multiplicative group (Z/nZ)^*. 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 znorder(GEN x, GEN o = NULL). Also available is GEN order(GEN x).

#### znprimroot(n)

returns a primitive root (generator) of (Z/nZ)^*, 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, then the smallest positive primitive root is returned. This is no longer true for composites.

Note that this function requires factoring p-1 for p as above, in order to determine the exact order of elements in (Z/nZ)^*: this is likely to be very costly if p is large. The function accepts vector/matrices arguments, and is then applied componentwise.

The library syntax is GEN znprimroot0(GEN n). For a t_INT x, the special case GEN znprimroot(GEN n) is also available.

#### znstar(n)

gives the structure of the multiplicative group (Z/nZ)^* 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/nZ)^* ~ 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.

The library syntax is GEN znstar(GEN n).