Reduced Binary Quadratic Forms

Over the past few weeks I've been writing some code for computing reduced binary quadratic forms:

(a, b, c) : = ax^2 + bxy + cy^2

The discriminant of a form is d = b^2 - 4ac. The code I wrote works for d < 0.

Such a form is reduced if |b| <= a <= c and |b| >= 0 if a = c or a = |b|.

You can see the code here:

https://github.com/wbhart/flint2/tree/qfb/qfb

There are two versions of the function qfb_reduced_forms.

Version 1

The first version works by iterating over all possible values "a" (it's a theorem that |b| <= a <= sqrt(|d|/3)) and c <= (1 - d)/4) and searching for valid triples (a, b, c).

To find valid triples, we note that if d = b^2 - 4ac then d = b^2 mod 4a. So we simply solve this quadratic congruence and search for all roots "b" in the correct range.

To speed things up, we factorise all the different possible values "a" by sieving. This makes the square root code much faster, as it doesn't have to factorise 4a before computing square roots.

This approach requires softly O(|d|^{1/4}) primes. There are about O(|d|^{1/2}) forms on average, and indeed the run time is softly O(|d|^{1/2}) (subject to the Generalised Riemann Hypothesis -- needed to guarantee we can find quadratic nonresidues quickly enough for the Tonelli-Shanks modular square root algorithm).

Version 2

The second approach iterates over all possible values |b| (the same bound applies as for "a").

This time we factorise (d - b^2)/4 for each "b" into products "ac". Again, we do this by sieving (this time quadratic sieving, which first finds square roots of "d" mod various primes "p"). 

However, we must sieve with primes "p" dividing values (d - b^2)/4. In other words we need softly O(|d|^{1/2}) primes. On a 64 bit machine, this exhausts the precomputed primes at about |d| = 100000000, so this technique is a bit limited.

The advantage of this technique is that it is about 50% faster as implemented. It still takes time softly O(|d|^{1/2}) though (again subject to GRH).

Timings

Magma computes reduced binary quadratic forms, so I did two comparisons:

Comparison A) Compute all reduced forms for discriminants 0 > d > -1000000.

Comparison B) Compute all reduced forms for disciminants -1000000000 > d > -1000000100.

Here are the timings:

Comparison	flint/qfb	Magma
A	118s	947s
B	1s	120s