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

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

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

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

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

  ? QFBclassno(D) = qfbclassno(D) * if (D > 0 && quadunitnorm(D) > 0, 2, 1)
  ? QFBclassno(136)
  %1 = 4
  ? qfbclassno(136)
  %2 = 2
  ? quadunitnorm(136)
  %3 = 1
  ? bnfnarrow(bnfinit(x^2 - 136)).cyc
  %4 = [4]  \\ narrow class group is cyclic ~ Z/4Z

Note that the use of bnfnarrow above is only valid because 136 is a fundamental discriminant: that function is asymptotically faster (and returns the group structure, not only its order) but only supports maximal orders. Here are a few more examples:

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

See also qfbhclassno.

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

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

  ? x=Qfb(2,3,-10);y=Qfb(5,3,-4);
  ? qfbcomp(x,y)
  %2 = Qfb(-2, 9, 1)
  ? qfbcomp(x,y)==qfbred(qfbcompraw(x,y))
  %3 = 1

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

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

  ? x=Qfb(2,3,-10);y=Qfb(5,3,-4);
  ? qfbcompraw(x,y)
  %2 = Qfb(10, 3, -2)
  ? x=Qfb(2,3,-10);y=Qfb(1,-1,1);
  ? qfbcompraw(x,y)
    ***   at top-level: qfbcompraw(x,y)
    ***                 ^ —  —  —  —  — 
    *** qfbcompraw: inconsistent qfbcompraw t_QFB , t_QFB.

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

Solves the equation x² + dy² = n in integers x and y, where d > 0 and n is prime. Returns the empty vector [] when no solution exists. It is also allowed to try n = 4 times a prime but the answer is then guaranteed only if d is 3 mod 4; more precisely if d != 3 mod 4, the algorithm may fail to find a non-primitive solution.

This function is a special case of qfbsolve applied to the principal form in the imaginary quadratic order of discriminant -4d (returning the solution with non-negative x and y). As its name implies, qfbcornacchia uses Cornacchia's algorithm and runs in time quasi-linear in log n (using halfgcd); in practical ranges, qfbcornacchia should be about twice faster than qfbsolve unless we indicate to the latter that its second argument is prime (see below).

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

We used Mat([p,1]) to indicate that p¹ was the integer factorization of p, i.e., that p is prime. Without it, qfbsolve attempts to factor p and wastes a little time.

The library syntax is GEN qfbcornacchia(GEN d, GEN n).

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

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

The library syntax is GEN hclassno(GEN x).

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

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

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

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

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

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

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

The library syntax is GEN qfbpow(GEN x, GEN n).

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

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

Prime binary quadratic form of discriminant x whose first coefficient is p, where |p| is a prime number. By abuse of notation, p = ± 1 is also valid and returns the unit form. Returns an error if x is not a quadratic residue mod p, or if x < 0 and p < 0. (Negative definite t_QFB are not implemented.)

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

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

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

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

Reduction of the (real or imaginary) binary quadratic form x, returns [y,g] where y is reduced and g in SL(2,ℤ) is such that g.x = y; isD, if present, must be equal to sqrtint(D), where D > 0 is the discriminant of x.

The action of g on x can be computed using qfeval(x,g)

  ? q1 = Qfb(33947,-39899,11650);
  ? [q2,U] = qfbredsl2(q1)
  %2 = [Qfb(749,2207,-1712),[-1,3;-2,5]]
  ? qfeval(q1,U)
  %3 = Qfb(749,2207,-1712)

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

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

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

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

The integer n can also be given by its factorization matrix fa = factor(n) or by the pair [n, fa].

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

For fixed Q, assuming the factorisation of n is given, the algorithm runs in probabilistic polynomial time in log p, where p is the largest prime divisor of n, through the computation of square roots of D modulo 4 p). The dependency on Q is more complicated: polynomial time in log |D| if Q is imaginary, but exponential time if Q is real (through the computation of a full cycle of reduced forms). In the latter case, note that bnfisprincipal provides a solution in heuristic subexponential time assuming the GRH.

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

Buchmann-McCurley's sub-exponential algorithm for computing the class group of a quadratic order of discriminant D. By default, the results are conditional on the GRH.

This function should be used instead of qfbclassno or quadregulator when D < -10²⁵, D > 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).

* v.normfu (for positive D only) return the norm of the fundamental unit, either 1 or -1. Note that a result of -1 is unconditional and no longer depends on the GRH.

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₁, c₂], where c₁ ≤ c₂ are nonnegative real numbers which control the execution time and the stack size, see se:GRHbnf. The parameter is used as a threshold to balance the relation finding phase against the final linear algebra. Increasing the default c₁ means that relations are easier to find, but more relations are needed and the linear algebra will be harder. The default value for c₁ is 0 and means that it is taken equal to c₂. The parameter c₂ is mostly obsolete and should not be changed, but we still document it for completeness: we compute a tentative class group by generators and relations using a factorbase of prime ideals ≤ c₁ (log |D|)², then prove that ideals of norm ≤ c₂ (log |D|)² do not generate a larger group. By default an optimal c₂ is chosen, so that the result is provably correct under the GRH — a result of Grenié and Molteni states that c₂ = 23/6 ~ 3.83 is fine (and even c₂ = 15/4 ~ 3.75 for large |D| > 2.41 E8). But it is possible to improve on this algorithmically. You may provide a smaller c₂, it will be ignored (we use the provably correct one); you may provide a larger c₂ 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 will be more convenient to use GEN Buchquad(GEN D, double c1, double c2, long prec).

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

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

The library syntax is GEN quaddisc(GEN x).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Regulator of the quadratic order of positive discriminant D in time Õ(D^1/2) using the continued fraction algorithm. Raise an error if D is not a discriminant (fundamental or not) or if D is a square. The function quadclassunit is asymptotically faster (and also in practice for D > 10¹⁰ or so) but depends on the GRH.

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

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

If v is given, the variable name is used to display u else 'w' is used. The algorithm computes the continued fraction of (1 + sqrt{D}) / 2 or sqrt{D}/2 (see GTM 138, algorithm 5.7.2). Although the continued fraction length is only O(sqrt{D}), the function still runs in time Õ(D), in part because the output size is not polynomially bounded in terms of log D. See bnfinit and bnfunits for a better alternative for large D, running in time subexponential in log D and returning the fundamental units in compact form (as a short list of S-units of size O(log D)³ raised to possibly large exponents).

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

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

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

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

The conductor f can be given in factored form or as [f, factor(f)]:

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

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

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

The library syntax is GEN quadunitindex(GEN D, GEN f).

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

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

This function computes the parity of the continued fraction expansion and runs in time Õ(D^1/2). If D is fundamental, the function bnfinit is asymptotically faster but depends of the GRH. If D = df² is not fundamental, it will usually be faster to first compute quadunitindex(d, f). If it is even, the result is 1, else the result is quadunitnorm(d). The narrow class number of the order of discriminant D is equal to the class number if the unit norm is 1 and to twice the class number otherwise.

Important remark. Assuming GRH, using bnfinit is much faster, running in time subexponential in log D (instead of exponential for quadunitnorm). We give examples for the maximal order:

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

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

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

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

The library syntax is long quadunitnorm(GEN D).