The Jacobi test by the binary algorithm

The binary algorithm of the C code (Ratpoints.Sieve.jacobi) computes the Jacobi symbol (lcf / b*) for b ≥ 1, where b* is the largest divisor of b prime to 2 lcf, b.coprimePart (2 * lcf.natAbs) (Nat.coprimePart): each round on two numbers f and b (Ratpoints.Sieve.jacobiRound) keeps the symbol with its sign, taking out the factors 2 of one number with the sign of (2 / m) for the other, odd, number m (Ratpoints.Sieve.twoSign) and exchanging the two numbers with the sign of the reciprocity, and gives the symbol (f / b) when f and b are coprime and otherwise their gcd, which the algorithm divides out of b. For a polynomial of even degree the symbol is 1 on the denominator b of every point, since c_d is a square modulo every prime factor of b prime to c_d (Ratpoints.Sieve.isSquare_getD_degree_of_sq_eq_binaryForm); so the binary algorithm lets the denominator of every point through, and so does the test by tables, which computes the same symbol (Ratpoints.Sieve.run_jacobiTest).

Main definitions

• Ratpoints.Sieve.IsJacobiReduction: what the subtraction loop of the binary algorithm returns.

Main statements

• Ratpoints.Sieve.twoSign_eq: the sign of (2 / m).
• Ratpoints.Sieve.jacobiReduce_spec, Ratpoints.Sieve.jacobiRound_spec: the contracts of the subtraction loop and of a round.
• Ratpoints.Sieve.jacobi_spec, Ratpoints.Sieve.run_jacobi: the binary algorithm computes (lcf / b*).
• Ratpoints.Sieve.jacobiSym_eq_one_of_sq_eq_binaryForm: the symbol of c_d over a divisor of the denominator of a point of a curve of even degree prime to c_d is 1, and so is (c_d / b*) (Ratpoints.Sieve.jacobiSym_coprimePart_eq_one_of_sq_eq).
• Ratpoints.Sieve.run_jacobi_eq_one_of_sq_eq: the binary algorithm gives 1 on the denominator of every point of a curve of even degree.

Implementation notes

The contracts are proved with vcgen, with loop invariants that are named definitions, in a file with Mathlib, since they speak of the Jacobi symbol; the early return of the main loop of Ratpoints.Sieve.jacobi is part of the invariant of that loop.

The binary algorithm

theorem Ratpoints.Sieve.twoSign_eq {m : ℕ} (hm : m % 2 = 1) (s : ℤ) : twoSign m s = if m % 8 = 3 ∨ m % 8 = 5 then -s else s

theorem Ratpoints.Sieve.twoSign_mul_jacobiSym_div_two {x y : ℕ} (hx : x % 2 = 0) (hy : y % 2 = 1) (s : ℤ) : twoSign y s * jacobiSym (↑(x / 2)) y = s * jacobiSym (↑x) y

Halving an even x against an odd y multiplies the symbol by (2 / y), the factor that Ratpoints.Sieve.twoSign takes.

theorem Ratpoints.Sieve.ite_neg_mul_jacobiSym {n m : ℕ} (hn : n % 2 = 1) (hm : m % 2 = 1) (s : ℤ) : (if (n &&& m &&& 2 != 0) = true then -s else s) * jacobiSym (↑m) n = s * jacobiSym (↑n) m

The exchange of two odd numbers in the binary algorithm, with the sign of the reciprocity, -1 exactly when both are 3 (mod 4), as Ratpoints.Sieve.jacobiRound tests it on their bits.

def Ratpoints.Sieve.IsJacobiReduction (x y : ℕ) (s : ℤ) (r : ℕ × ℤ) : Prop

What the subtraction loop Ratpoints.Sieve.jacobiReduce returns on the odd x and y with the sign s: an odd x' ≤ x, y and a sign s' with s' J(x' | y) = s J(x | y) and gcd(x', y) = gcd(x, y).

Equations

• Ratpoints.Sieve.IsJacobiReduction x y s r = (r.1 % 2 = 1 ∧ r.1 ≤ x ∧ r.1 ≤ y ∧ r.2 * jacobiSym (↑r.1) y = s * jacobiSym (↑x) y ∧ r.1.gcd y = x.gcd y)

def Ratpoints.Sieve.JacobiReduceInvariant (x y : ℕ) (s : ℤ) : Std.WP.WhileInvariant (ℕ × ℤ) Prop

The invariant of the outer loop of Ratpoints.Sieve.jacobiReduce on the odd x and y with the sign s: the current x' is odd, at most x, has the gcd of x with y, and s' J(x' | y) = s J(x | y); at the exit x' ≤ y.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sieve.JacobiReduceHalveInvariant (x y x₀ : ℕ) (s : ℤ) : Std.WP.WhileInvariant (ℕ × ℤ) Prop

The invariant of the inner loop of Ratpoints.Sieve.jacobiReduce, which halves x₀ - y, for the value x₀ of the outer loop: the current x' is positive and below x₀, keeps the symbol and the gcd of x and y, and is even while the loop goes on and odd at its exit.

Equations

• One or more equations did not get rendered due to their size.

theorem Ratpoints.Sieve.jacobiReduceInvariant_self {x y : ℕ} (s : ℤ) (h : x % 2 = 1 ∧ y % 2 = 1) : JacobiReduceInvariant x y s false (x, s)

theorem Ratpoints.Sieve.JacobiReduceInvariant.isJacobiReduction {x y x' : ℕ} {s s' : ℤ} (h : JacobiReduceInvariant x y s true (x', s')) : IsJacobiReduction x y s (x', s')

theorem Ratpoints.Sieve.JacobiReduceInvariant.exit {x y x' mb : ℕ} {s s' : ℤ} (hle : ¬y < x') (h : x' = mb ∧ JacobiReduceInvariant x y s false (x', s')) : JacobiReduceInvariant x y s true (x', s')

theorem Ratpoints.Sieve.JacobiReduceInvariant.jacobiReduceHalveInvariant {x y x₀ mb : ℕ} {s s₀ : ℤ} (hy : x % 2 = 1 ∧ y % 2 = 1) (hlt : y < x₀) (h : x₀ = mb ∧ JacobiReduceInvariant x y s false (x₀, s₀)) : JacobiReduceHalveInvariant x y x₀ s false (x₀ - y, s₀)

theorem Ratpoints.Sieve.JacobiReduceHalveInvariant.div_two_of_mod_two_eq {x y x₀ x' mb : ℕ} {s s' : ℤ} {e : Bool} (hy : x % 2 = 1 ∧ y % 2 = 1) (h : x' = mb ∧ JacobiReduceHalveInvariant x y x₀ s false (x', s')) (he : x' / 2 % 2 = if e = true then 1 else 0) : JacobiReduceHalveInvariant x y x₀ s e (x' / 2, twoSign y s')

theorem Ratpoints.Sieve.JacobiReduceHalveInvariant.exit {x y x₀ x' mb : ℕ} {s s' : ℤ} (hy : x % 2 = 1 ∧ y % 2 = 1) (h : x' = mb ∧ JacobiReduceHalveInvariant x y x₀ s false (x', s')) (hb : (x' / 2 % 2 == 1) = true) : JacobiReduceHalveInvariant x y x₀ s true (x' / 2, twoSign y s')

theorem Ratpoints.Sieve.JacobiReduceHalveInvariant.div_two {x y x₀ x' mb : ℕ} {s s' : ℤ} (hy : x % 2 = 1 ∧ y % 2 = 1) (h : x' = mb ∧ JacobiReduceHalveInvariant x y x₀ s false (x', s')) (hb : ¬(x' / 2 % 2 == 1) = true) : x' / 2 < mb ∧ JacobiReduceHalveInvariant x y x₀ s false (x' / 2, twoSign y s')

theorem Ratpoints.Sieve.JacobiReduceHalveInvariant.jacobiReduceInvariant {x y x₀ x' mb : ℕ} {s s' s₀ : ℤ} (h : JacobiReduceHalveInvariant x y x₀ s true (x', s')) (h₀ : x₀ = mb ∧ JacobiReduceInvariant x y s false (x₀, s₀)) : x' < mb ∧ JacobiReduceInvariant x y s false (x', s')

theorem Ratpoints.Sieve.jacobiReduce_spec (x y : ℕ) (s : ℤ) : ⦃ x % 2 = 1 ∧ y % 2 = 1 ⦄ jacobiReduce x y s ⦃ IsJacobiReduction x y s ⦄

The subtraction loop of the binary algorithm keeps the symbol with its sign and the gcd, and ends with an odd x' ≤ y; its loops end, since x' drops with every pass.

def Ratpoints.Sieve.JacobiRoundHalveInvariant (f b : ℕ) : Std.WP.WhileInvariant (ℤ × ℕ) Prop

The invariant of the loop of Ratpoints.Sieve.jacobiRound on f > 0 and the odd b that takes the factors 2 out of n = f: n is positive, s J(n | b) = J(f | b), the gcd with b is that of f, and n is odd at the exit.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sieve.JacobiRoundInvariant (f b : ℕ) : Std.WP.WhileInvariant (ℤ × ℕ × ℕ) Prop

The invariant of the main loop of Ratpoints.Sieve.jacobiRound on f and b, on the sign s and the odd n and m: s J(n | m) = J(f | b), the gcd of n and m is that of f and b, and n = m at the exit.

Equations

• One or more equations did not get rendered due to their size.

theorem Ratpoints.Sieve.jacobiRoundHalveInvariant_one {f b : ℕ} (h : 0 < f ∧ b % 2 = 1) : JacobiRoundHalveInvariant f b false (1, f)

theorem Ratpoints.Sieve.JacobiRoundHalveInvariant.div_two {f b n mb : ℕ} {s : ℤ} (h : 0 < f ∧ b % 2 = 1) (hn : (n % 2 == 0) = true) (hi : n = mb ∧ JacobiRoundHalveInvariant f b false (s, n)) : n / 2 < mb ∧ JacobiRoundHalveInvariant f b false (twoSign b s, n / 2)

theorem Ratpoints.Sieve.JacobiRoundHalveInvariant.exit {f b n mb : ℕ} {s : ℤ} (hn : ¬(n % 2 == 0) = true) (hi : n = mb ∧ JacobiRoundHalveInvariant f b false (s, n)) : JacobiRoundHalveInvariant f b true (s, n)

theorem Ratpoints.Sieve.JacobiRoundHalveInvariant.jacobiRoundInvariant {f b n : ℕ} {s : ℤ} (h : 0 < f ∧ b % 2 = 1) (hi : JacobiRoundHalveInvariant f b true (s, n)) : JacobiRoundInvariant f b false (s, n, b)

theorem Ratpoints.Sieve.JacobiRoundInvariant.eq_gcd {f b n m : ℕ} {s : ℤ} (h : JacobiRoundInvariant f b true (s, n, m)) : m = f.gcd b ∧ (m = 1 → s = jacobiSym (↑f) b)

theorem Ratpoints.Sieve.JacobiRoundInvariant.odd {f b n m mb : ℕ} {s : ℤ} (h : n + m = mb ∧ JacobiRoundInvariant f b false (s, n, m)) : m % 2 = 1 ∧ n % 2 = 1

theorem Ratpoints.Sieve.JacobiRoundInvariant.odd_of_mod_two_eq_one {f b n m m' mb : ℕ} {s : ℤ} (h : n + m = mb ∧ JacobiRoundInvariant f b false (s, n, m)) (hm' : m' % 2 = 1) : n % 2 = 1 ∧ m' % 2 = 1

theorem Ratpoints.Sieve.IsJacobiReduction.mul_jacobiSym_eq_and_gcd_eq {a c c' g : ℕ} {s s' t j : ℤ} (hr : IsJacobiReduction c a t (c', s')) (ha : a % 2 = 1) (hc : c % 2 = 1) (hj : s * jacobiSym (↑a) c = j) (hg : a.gcd c = g) (ht : t = if (a &&& c &&& 2 != 0) = true then -s else s) : s' * jacobiSym (↑c') a = j ∧ c'.gcd a = g

The exchange of the odd a and c with the sign of the reciprocity, followed by the subtraction loop Ratpoints.Sieve.jacobiReduce of c by a, keeps the symbol with its sign and the gcd; a pass of the main loop of Ratpoints.Sieve.jacobiRound takes two such steps.

theorem Ratpoints.Sieve.JacobiRoundInvariant.jacobiReduce_of_beq {f b n m m' mb : ℕ} {s s' t : ℤ} (h : n + m = mb ∧ JacobiRoundInvariant f b false (s, n, m)) (hr : IsJacobiReduction m n t (m', s')) (heq : (m' == n) = true) (ht : t = if (n &&& m &&& 2 != 0) = true then -s else s) : JacobiRoundInvariant f b true (s', n, m')

theorem Ratpoints.Sieve.JacobiRoundInvariant.jacobiReduce_jacobiReduce_of_beq {f b n m m' n' mb : ℕ} {s s' s'' t t' : ℤ} (h : n + m = mb ∧ JacobiRoundInvariant f b false (s, n, m)) (hr : IsJacobiReduction m n t (m', s')) (hr' : IsJacobiReduction n m' t' (n', s'')) (heq : (m' == n') = true) (ht : t = if (n &&& m &&& 2 != 0) = true then -s else s) (ht' : t' = if (n &&& m' &&& 2 != 0) = true then -s' else s') : JacobiRoundInvariant f b true (s'', n', m')

theorem Ratpoints.Sieve.JacobiRoundInvariant.jacobiReduce_jacobiReduce {f b n m m' n' mb : ℕ} {s s' s'' t t' : ℤ} (h : n + m = mb ∧ JacobiRoundInvariant f b false (s, n, m)) (hr : IsJacobiReduction m n t (m', s')) (hne : ¬(m' == n) = true) (hr' : IsJacobiReduction n m' t' (n', s'')) (ht : t = if (n &&& m &&& 2 != 0) = true then -s else s) (ht' : t' = if (n &&& m' &&& 2 != 0) = true then -s' else s') : n' + m' < mb ∧ JacobiRoundInvariant f b false (s'', n', m')

theorem Ratpoints.Sieve.jacobiRound_spec (f b : ℕ) : ⦃ 0 < f ∧ b % 2 = 1 ⦄ jacobiRound f b ⦃ fun (r : ℕ × ℤ) => r.1 = f.gcd b ∧ (r.1 = 1 → r.2 = jacobiSym (↑f) b) ⦄

A round of the binary algorithm on f > 0 and the odd b gives their gcd, with the symbol J(f | b) when it is 1; its main loop ends, since n + m drops with every pass.

def Ratpoints.Sieve.JacobiOddInvariant (b : ℕ) (lcf : ℤ) : Std.WP.WhileInvariant ℕ Prop

The invariant of the loop of Ratpoints.Sieve.jacobi that takes the factors 2 out of b: the current b' is positive with the part of b prime to 2 lcf, and odd at the exit.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sieve.JacobiInvariant (b : ℕ) (lcf : ℤ) : Std.WP.WhileInvariant (Option ℤ × ℕ × ℕ) Prop

The invariant of the main loop of Ratpoints.Sieve.jacobi on b and lcf, on the result of an early return, the current b' and f: without a result, b' is odd and positive, f < b', f ≡ lcf (mod b'), the part of b' prime to 2 lcf is that of b, and f = 0 at the exit; a result comes only at the exit and is the symbol (lcf / b*) of the part b* of b prime to 2 lcf.

Equations

• One or more equations did not get rendered due to their size.

theorem Ratpoints.Sieve.jacobiOddInvariant_self {b : ℕ} (lcf : ℤ) (hb : 0 < b) : JacobiOddInvariant b lcf false b

theorem Ratpoints.Sieve.JacobiOddInvariant.div_two {b b' mb : ℕ} {lcf : ℤ} (h : (b' % 2 == 0) = true) (hi : b' = mb ∧ JacobiOddInvariant b lcf false b') : b' / 2 < mb ∧ JacobiOddInvariant b lcf false (b' / 2)

theorem Ratpoints.Sieve.JacobiOddInvariant.exit {b b' mb : ℕ} {lcf : ℤ} (h : ¬(b' % 2 == 0) = true) (hi : b' = mb ∧ JacobiOddInvariant b lcf false b') : JacobiOddInvariant b lcf true b'

theorem Ratpoints.Sieve.JacobiOddInvariant.jacobiInvariant {b b' : ℕ} {lcf : ℤ} (hi : JacobiOddInvariant b lcf true b') : JacobiInvariant b lcf false (none, b', (lcf % ↑b').toNat)

theorem Ratpoints.Sieve.JacobiInvariant.eq_jacobiSym_of_eq_some {b b' f : ℕ} {lcf r : ℤ} {o : Option ℤ} (hi : JacobiInvariant b lcf true (o, b', f)) (h : o = some r) : r = jacobiSym lcf (b.coprimePart (2 * lcf.natAbs))

theorem Ratpoints.Sieve.JacobiInvariant.one_eq_jacobiSym_of_eq_none {b b' f : ℕ} {lcf : ℤ} {o : Option ℤ} (hi : JacobiInvariant b lcf true (o, b', f)) (h : o = none) : 1 = jacobiSym lcf (b.coprimePart (2 * lcf.natAbs))

At the exit of Ratpoints.Sieve.jacobi without a result, b' ∣ lcf, so that the part of b prime to 2 lcf, that of b', is 1.

theorem Ratpoints.Sieve.JacobiInvariant.eq_none {b b' f : ℕ} {lcf : ℤ} {o : Option ℤ} (hi : JacobiInvariant b lcf false (o, b', f)) : o = none

theorem Ratpoints.Sieve.JacobiInvariant.pos_and_odd {b b' f mb : ℕ} {lcf : ℤ} {o : Option ℤ} (hf : (f != 0) = true) (hi : b' = mb ∧ JacobiInvariant b lcf false (o, b', f)) : 0 < f ∧ b' % 2 = 1

theorem Ratpoints.Sieve.JacobiInvariant.exit {b b' f mb : ℕ} {lcf : ℤ} {o : Option ℤ} (hf : ¬(f != 0) = true) (hi : b' = mb ∧ JacobiInvariant b lcf false (o, b', f)) : JacobiInvariant b lcf true (none, b', f)

theorem Ratpoints.Sieve.JacobiInvariant.exit_some {b b' f g mb : ℕ} {lcf s : ℤ} {o : Option ℤ} (hr : g = f.gcd b' ∧ (g = 1 → s = jacobiSym (↑f) b')) (hg : (g == 1) = true) (hi : b' = mb ∧ JacobiInvariant b lcf false (o, b', f)) : JacobiInvariant b lcf true (some s, b', f)

When gcd(f, b') = 1, the symbol of the round is J(lcf | b'), and b' is the part of b prime to 2 lcf, being odd and prime to lcf.

theorem Ratpoints.Sieve.JacobiInvariant.div_gcd {b b' f g f' mb : ℕ} {lcf s : ℤ} {o : Option ℤ} (hf : (f != 0) = true) (hr : g = f.gcd b' ∧ (g = 1 → s = jacobiSym (↑f) b')) (hg : ¬(g == 1) = true) (hi : b' = mb ∧ JacobiInvariant b lcf false (o, b', f)) (hf' : f' = if b' / g ≤ f then f % (b' / g) else f) : b' / g < mb ∧ JacobiInvariant b lcf false (none, b' / g, f')

Dividing b' by g = gcd(f, b') > 1 keeps the part of b' prime to 2 lcf, since the primes of g divide lcf, and f ≡ lcf modulo the quotient, which stays odd.

theorem Ratpoints.Sieve.jacobi_spec (b : ℕ) (lcf : ℤ) : ⦃ 0 < b ⦄ jacobi b lcf ⦃ fun (r : ℤ) => r = jacobiSym lcf (b.coprimePart (2 * lcf.natAbs)) ⦄

The binary algorithm. Ratpoints.Sieve.jacobi b lcf is the Jacobi symbol (lcf / b*) for b ≥ 1, b* the largest divisor of b prime to 2 lcf; its main loop ends, since every round that does not return divides b by a gcd above 1.

theorem Ratpoints.Sieve.run_jacobi {b : ℕ} (hb : b ≠ 0) (lcf : ℤ) : (jacobi b lcf).run = jacobiSym lcf (b.coprimePart (2 * lcf.natAbs))

The denominators of the points

theorem Ratpoints.Sieve.jacobiSym_eq_one_of_sq_eq_binaryForm {c : List ℤ} {a y : ℤ} {b b' : ℕ} (hab : a.gcd ↑b = 1) (hd : degree c % 2 = 0) (hy : y ^ 2 = binaryForm c a ↑b) (hb' : b' ∣ b) (hc : (c.getD (degree c) 0).gcd ↑b' = 1) : jacobiSym (c.getD (degree c) 0) b' = 1

Soundness of every form of the Jacobi test. For a point (a, y, b) of a curve of even degree, every divisor b' of b prime to the leading coefficient c_d has J(c_d | b') = 1: a prime q ∣ b' divides b, so c_d is a square modulo q (Ratpoints.Sieve.isSquare_getD_degree_of_sq_eq_binaryForm), and a Jacobi symbol -1 would have a prime factor q with J(c_d | q) = -1 (jacobiSym.eq_neg_one_at_prime_divisor_of_eq_neg_one).

theorem Ratpoints.Sieve.jacobiSym_coprimePart_eq_one_of_sq_eq {c : List ℤ} (hd : degree c % 2 = 0) {a y : ℤ} {b : ℕ} (hab : a.gcd ↑b = 1) (hy : y ^ 2 = binaryForm c a ↑b) : jacobiSym (c.getD (degree c) 0) (b.coprimePart (2 * (c.getD (degree c) 0).natAbs)) = 1

The symbol (c_d / b*) of the largest divisor b* of the denominator b of a point prime to 2 c_d is 1, for even degree.

theorem Ratpoints.Sieve.run_jacobi_eq_one_of_sq_eq {c : List ℤ} (hd : degree c % 2 = 0) {a y : ℤ} {b : ℕ} (hb : b ≠ 0) (hab : a.gcd ↑b = 1) (hy : y ^ 2 = binaryForm c a ↑b) : (jacobi b (c.getD (degree c) 0)).run = 1

The binary algorithm gives the symbol 1 on the denominator b of every point of a curve of even degree: it computes the Jacobi symbol (c_d / b*), and modulo a prime q ∣ b with q ∤ c_d, c_d is the square of y / a ^ (d / 2).