### Arithmetic functions

#### 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 divise by the elements in this table. If x is empty or omitted, just returns the current list of extra primes.

The entries in x are not checked for primality, and in fact they need only be positive integers. The algorithm makes sure that all elements in the table are pairwise coprime, so it may end up containing divisors of the input integers.

It is a useful trick to add known composite numbers, which the function `factor`(x,0) was not able to factor. In case the message "impossible inverse modulo `<`some INTMOD`>`" shows up afterwards, you have just stumbled over a non-trivial factor. Note that the arithmetic functions in the narrow sense, like eulerphi, do not use this extra table.

To remove primes from the list use `removeprimes`.

The library syntax is addprimes(x).

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

if B is omitted, finds the best rational approximation to x belongs to R (or R[X], or R^n,...) with denominator at most equal to A using continued fractions.

If B is present, x is assumed to be of type `t_INTMOD` modulo M (or a recursive combination of those), and the routine returns the unique fraction a/b in coprime integers a `<=` A and b `<=` B which is congruent to x modulo M. If M `<=` 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.

The library syntax is bestappr0(x,A,B). Also available is `bestappr(x,A)` corresponding to an omitted B.

#### bezout(x,y)

finds u and v minimal in a natural sense such that x*u+y*v = gcd(x,y). The arguments must be both integers or both polynomials, and the result is a row vector with three components u, v, and gcd(x,y).

The library syntax is vecbezout(x,y) to get the vector, or `gbezout(x,y, &u, &v)` which gives as result the address of the created gcd, and puts the addresses of the corresponding created objects into u and v.

#### bezoutres(x,y)

as `bezout`, with the resultant of x and y replacing the gcd. The algorithm uses (subresultant) assumes the base ring is a domain.

The library syntax is vecbezoutres(x,y) to get the vector, or `subresext(x,y, &u, &v)` which gives as result the address of the created gcd, and puts the addresses of the corresponding created objects into u and v.

#### bigomega(x)

number of prime divisors of |x| counted with multiplicity. x must be an integer.

The library syntax is bigomega(x), the result is a `long`.

#### 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 binomial(x,y), where y must be a `long`.

#### 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 chinese(x,y). Also available is `chinese1`(x), corresponding to an ommitted `y`.

#### 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 content(x).

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

creates the row vector whose components are the partial quotients of the continued fraction expansion of x. That is 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 is limited to nmax. If x is a real number, the expansion stops at the last significant partial quotient if nmax is omitted. x can also be a rational function or a power series.

If a vector b is supplied, the numerators will be equal to the coefficients of b (instead of all equal to 1 as above). The length of the result is then equal to the length of b, unless a partial remainder is encountered which is equal to zero, in which case the expansion stops. In the case of real numbers, the stopping criterion is thus different from the one mentioned above since, if b is too long, some partial quotients may not be significant.

If b is an integer, the command is understood as `contfrac(x,nmax)`.

The library syntax is contfrac0(x,b,nmax). Also available are `gboundcf(x,nmax)`, `gcf(x)`, or `gcf2(b,x)`, where nmax is a C integer.

#### 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 pnqn(x).

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

if n is a non-zero integer written as n = df^2 with d squarefree, returns d. If flag is non-zero, returns the two-element row vector [d,f].

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

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

if n is a non-zero integer written as n = df^2 with d fundamental discriminant (including 1), returns d. If flag is non-zero, returns the two-element row vector [d,f]. Note that if n is not congruent to 0 or 1 modulo 4, f will be a half integer and not an integer.

The library syntax is coredisc0(n,flag). Also available are `coredisc(n)` ( = `coredisc(n,0)`) and `coredisc2(n)` ( = `coredisc(n,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 dirdiv(x,y).

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

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)(GEN,void*), 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 dirmul(x,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 divisors(x).

#### eulerphi(x)

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

The library syntax is phi(x).

#### factor(x,{lim = -1})

general factorization function. If x is of type integer, rational, polynomial or rational function, the result is a two-column matrix, the first column being the irreducibles dividing x (prime numbers or polynomials), and the second the exponents. If x is a vector or a matrix, the factoring is done componentwise (hence the result is a vector or matrix of two-column matrices). By definition, 0 is factored as 0^1.

If x is of type integer or rational, the factors are pseudoprimes (see `ispseudoprime`), and in general not rigorously proven primes. In fact, any factor which is `<=` 10^{13} is a genuine prime number. Use `isprime` to prove primality of other factors, as in

    fa = factor(2^2^7 +1)
     isprime( fa[,1] )

An argument lim can be added, meaning that we look only for prime factors p < lim, or up to `primelimit`, whichever is lowest (except when lim = 0 where the effect is identical to setting lim = `primelimit`). In this case, the remaining part may actually be a proven composite! See `factorint` for more information about the algorithms used.

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.

Note that PARI tries to guess in a sensible way over which ring you want to factor. Note also 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 `factornf` and `nffactor`.

The library syntax is factor0(x,lim), where lim is a C integer. Also available are `factor(x)` ( = `factor0(x,-1)`), `smallfact(x)` ( = `factor0(x,0)`).

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

gives back the factored object corresponding to a factorization. The integer 1 corresponds to the empty factorization. If the last argument is of number field type (e.g.created by `nfinit`), assume we are dealing with an ideal factorization in the number field. The resulting ideal product is given in HNF form.

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 1 (the product of the f[i] is returned). Finally, f can be a regular factorization, as produced with any `factor` command. A few examples:

    ? factorback([2,2; 3,1])
     %1 = 12
     ? factorback([2,2], [3,1])
     %2 = 12
     ? factorback([5,2,3])
     %3 = 30
     ? factorback([2,2], [3,1], nfinit(x^3+2))
     %4 =
     [16 0 0]

     [0 16 0]

     [0 0 16]
     ? nf = nfinit(x^2+1); fa = idealfactor(nf, 10)
     %5 =
     [[2, [1, 1]~, 2, 1, [1, 1]~] 2]

     [[5, [-2, 1]~, 1, 1, [2, 1]~] 1]

     [[5, [2, 1]~, 1, 1, [-2, 1]~] 1]
     ? factorback(fa)
       ***   forbidden multiplication t_VEC * t_VEC.
     ? factorback(fa, nf)
     %6 =
     [10 0]

     [0 10]

In the fourth example, 2 and 3 are interpreted as principal ideals in a cubic field. In the fifth one, `factorback(fa)` is meaningless since we forgot to indicate the number field, and the entries in the first column of `fa` can't be multiplied.

The library syntax is factorback0(f,e,nf), where an omitted nf or e is entered as `NULL`. Also available is `factorback`(f,nf) (case e = `NULL`) where an omitted nf is entered as `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 factcantor(x,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.

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

#### factorial(x)

or 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 mpfact(x) for x! and `mpfactr(x,prec)` for `factorial`(x). x must be a `long` integer and not a PARI integer.

#### factorint(n,{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 with exponents `<=` 10. The output is a two-column matrix as for `factor`. Use `isprime` on the result if you want to guarantee primality.

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 (very rarely) not be detected.

You are invited to play with the flag settings and watch the internals at work by using `gp`'s `debuglevel` default parameter (level 3 shows just the outline, 4 turns on time keeping, 5 and above show an increasing amount of internal details). If you see anything funny happening, please let us know.

The library syntax is factorint(n,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 factormod(x,p,flag). Also available are `factmod(x,p)` (which is equivalent to `factormod(x,p,0)`) and `simplefactmod(x,p)` ( = `factormod(x,p,1)`).

#### fibonacci(x)

x^{{th}} Fibonacci number.

The library syntax is fibo(x). x must be a `long`.

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

computes a monic polynomial of degree n which is irreducible over F_p. For instance if `P = ffinit(3,2,y)`, you can represent elements in F_{3^2} as polmods modulo `P`. This function uses a fast variant of Adleman-Lenstra's algorithm.

The library syntax is ffinit(p,n,v), where v is a variable number.

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

creates the greatest common divisor of x and y. x and y can be of 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 coeffients in the same number field (in particular rational): use modular gcd algorithm.

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

The library syntax is ggcd(x,y). For general polynomial inputs, `srgcd(x,y)` is also available. For univariate rational polynomials, one also has `modulargcd(x,y)`.

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

Hilbert symbol of x and y modulo p. If x and y are of type integer or fraction, an explicit third parameter p must be supplied, p = 0 meaning the place at infinity. Otherwise, p needs not be given, and x and y can be of compatible types integer, fraction, real, intmod a prime (result is undefined if the modulus is not prime), or p-adic.

The library syntax is hil(x,y,p).

#### isfundamental(x)

true (1) if x is equal to 1 or to the discriminant of a quadratic field, false (0) otherwise.

The library syntax is gisfundamental(x), but the simpler function `isfundamental(x)` which returns a `long` should be used if x is known to be of type integer.

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

if k is given, returns true (1) if x is a k-th power, false (0) if not. In this case, x may be an integer or polynomial, a rational number or function, or an intmod a prime or p-adic.

If k is omitted, only integers and fractions are allowed 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 a k-th root was computed in the process, then n is set to that root.

The library syntax is ispower(x, k, &n), the result is a `long`. Omitted k or n are coded as `NULL`.

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

true (1) if x is a (proven) prime number, false (0) otherwise. This can be very slow when x is indeed prime and has more than 1000 digits, say. Use `ispseudoprime` to quickly check for pseudo primality. See also `factor`.

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 < 10^{13}), 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 the prime factors p of x-1, 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).

If flag = 2, use APRCL.

The library syntax is gisprime(x,flag), but the simpler function `isprime(x)` which returns a `long` should be used if x is known to be of type integer.

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

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 (in particular, all composites `<=` 10^{13} are correctly detected), although it is expected that infinitely many such numbers exist.

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 gispseudoprime(x,flag), but the simpler function `ispseudoprime(x)` which returns a `long` should be used if x is known to be of type integer.

#### 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 and an exact square root had to be computed in the checking process, puts that square root in n. This is the case when x is a `t_INT`, `t_FRAC`, `t_POL` or `t_RFRAC` (or a vector of such objects):

        ? 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

This will not work for `t_INTMOD` (use quadratic reciprocity) or `t_SER` (only check the leading coefficient).

The library syntax is gissquarerem(x,&n). Also available is `gissquare(x)`.

#### issquarefree(x)

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

The library syntax is gissquarefree(x), but the simpler function `issquarefree(x)` which returns a `long` should be used if x is known to be of type integer. This issquarefree is just the square of the Moebius function, and is computed as a multiplicative arithmetic function much like the latter.

#### 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 kronecker(x,y), the result (0 or ± 1) is a `long`.

#### 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 glcm(x,y).

#### moebius(x)

Moebius mu-function of |x|. x must be of type integer.

The library syntax is mu(x), the result (0 or ± 1) is a `long`.

#### nextprime(x)

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

The library syntax is nextprime(x).

#### numdiv(x)

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

The library syntax is numbdiv(x).

#### numbpart(n)

gives the number of unrestricted partitions of n, usually called p(n) in the litterature; in other words the number of nonnegative integer solutions to a+2b+3c+.. .= n. n must be of type integer and 1 `<=` n < 10^{15}. The algorithm uses the Hardy-Ramanujan-Rademacher formula.

The library syntax is numbpart(n).

#### omega(x)

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

The library syntax is omega(x), the result is a `long`.

#### precprime(x)

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

The library syntax is precprime(x).

#### prime(x)

the x^{{th}} prime number, which must be among the precalculated primes.

The library syntax is prime(x). x must be a `long`.

#### 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 primepi(x).

#### primes(x)

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

The library syntax is primes(x). x must be a `long`.

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

ordinary class number of the quadratic order of discriminant D. In the present version 2.2.13, 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 a low exponent (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

The library syntax is qfbclassno0(D,flag). Also available: `classno(D)` ( = `qfbclassno(D)`), `classno2(D)` ( = `qfbclassno(D,1)`), and finally we have the function `hclassno(D)` which computes the class number of an imaginary quadratic field by counting reduced forms, an O(|D|) algorithm. See also `qfbhclassno`.

#### 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 compraw(x,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 hclassno(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 take advantadge of asymptotically fast operations and careful optimizations).

The library syntax is nucomp(x,y,l). The auxiliary function `nudupl(x,l)` can be used 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 nupow(x,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 powraw(x,n) where n must be a `long` integer.

#### qfbprimeform(x,p)

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

The library syntax is primeform(x,p,prec), where the third variable prec is a `long`, but is only taken into account when x > 0.

#### qfbred(x,{flag = 0},{D},{isqrtD},{sqrtD})

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

D, isqrtD, sqrtD, 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 qfbred0(x,flag,D,isqrtD,sqrtD). Use `NULL` to omit any of D, isqrtD, sqrtD.

Also available are

`redimag(x)` ( = `qfbred(x)` where x is definite),

and for indefinite forms:

`redreal(x)` ( = `qfbred(x)`),

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

`redrealnod(x,sq)` ( = `qfbred(x,2,,isqrtD)`),

`rhorealnod(x,sq)` ( = `qfbred(x,3,,isqrtD)`).

#### 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 qfbsolve(Q,n).

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

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

This function should be used instead of `qfbclassno` or `quadregula` when D < -10^{25}, D > 10^{10}, or when the structure is wanted. It is a special case of `bnfinit`, which is slower, but more robust.

If flag is non-zero and D > 0, computes the narrow class group and regulator, instead of the ordinary (or wide) ones. In the current version 2.2.13, this does 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. For a given c_1, set c_2 = c_1 to get maximum speed. To get a rigorous result under GRH, you must take c_2 `>=` 6. Reasonable values for c_1 are between 0.1 and 2. More precisely, the algorithm will assume that prime ideals of norm less than c_2 (log |D|)^2 generate the class group, but the bulk of the work is done with prime ideals of norm less than c_1 (log |D|)^2. A larger c_1 means that relations are easier to find, but more relations are needed and the linear algebra will be harder. The default is c_1 = c_2 = 0.2, so the result is not rigorously proven.

The result is a vector v with 3 components if D < 0, and 4 otherwise. The correspond respectively to

* v[1]: the class number

* v[2]: a vector giving the structure of the class group as a product of cyclic groups;

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

* v[4]: (omitted if D < 0) 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 library syntax is quadclassunit0(D,flag,tech). Also available are `buchimag(D,c_1,c_2)` and `buchreal(D,flag,c_1,c_2)`.

#### quaddisc(x)

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

The library syntax is quaddisc(x).

#### quadhilbert(D,{pq})

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

If D < 0, uses complex multiplication (Schertz's variant). The technical component pq, if supplied, is a vector [p,q] where p, q are the prime numbers needed for the Schertz's method. More precisely, prime ideals above p and q should be non-principal and coprime to all reduced representatives of the class group. In addition, if one of these ideals has order 2 in the class group, they should have the same class. Finally, for efficiency, gcd(24,(p-1)(q-1)) should be as large as possible. The routine returns 0 if [p,q] is not suitable.

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 quadhilbert(D,pq,prec).

#### quadgen(D)

creates the quadratic number omega = (a+sqrt{D})/2 where a = 0 if x = 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 quadgen(x).

#### 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 quadpoly0(x,v).

#### quadray(D,f,{lambda})

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. If supplied, lambda is is the technical element lambda of `bnf` necessary for Schertz's method. In that case, returns 0 if lambda is not suitable.

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 quadray(D,f,lambda,prec), where an omitted `lambda` is coded as `NULL`.

#### 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 regula(x,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 fundunit(x).

#### 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 removeprimes(x).

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

sum of the k^{{th}} powers of the positive divisors of |x|. x and k must be of type integer.

The library syntax is sumdiv(x) ( = `sigma(x)`) or `gsumdivk(x,k)` ( = `sigma(x,k)`), where k is a C long integer.

#### 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 sqrti(x). Also available is `sqrtremi`(x,&r) which returns s such that s^2 = x+r, with 0 `<=` r `<=` 2s.

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

finds all integers x_0 with |x_0| `<=` X such that gcd(N, P(x_0)) `>=` B. If N is prime or a prime power, `polrootsmod` or `polrootspadic` will be much faster. X must be smaller than exp(log^2 B / (deg(P) log N)).

The library syntax is zncoppersmith(P, N, X, B), where an omitted B is coded as `NULL`.

#### znlog(x,g)

g must be a primitive root mod a prime p, and the result is the discrete log of x in the multiplicative group (Z/pZ)^*. This function uses a simple-minded combination of Pohlig-Hellman algorithm and Shanks baby-step/giant-step which requires O(sqrt{q}) storage, where q is the largest prime factor of p-1. Hence it cannot be used when the largest prime divisor of p-1 is greater than about 10^{13}.

The library syntax is znlog(x,g).

#### 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. If optional parameter o is given it is assumed to be a multiple of the order (used to limit the search space).

The library syntax is znorder(x,o), where an omitted o is coded as `NULL`. Also available is `order(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).

The library syntax is gener(x).

#### 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 znstar(n).