The reductions by a reciprocal and of the row indices

The high half of the product of two 64-bit words from the products of their 32-bit halves (Ratpoints.Sieve.mulHi) is the high half of the product of natural numbers (Ratpoints.Sieve.toNat_mulHi), so that the reductions by a reciprocal on 64-bit words (Ratpoints.Sieve.mulMod, Ratpoints.Sieve.mulDiv) compute the remainder and the quotient where the C code relies on them: with the reciprocal Ratpoints.Sieve.magic p, which exceeds 2 ^ 64 / p by less than 1, for every u with u p < 2 ^ 64, and so for the values at most Ratpoints.mulModLimit and the moduli up to 2 ^ 32; the quotient for p ≥ 2, since the reciprocal 2 ^ 64 of 1 wraps to 0 on a word, as in the C code. So does the remainder of a value of either sign (Ratpoints.Sieve.modMul), and the helper mod of the C code (Ratpoints.Sieve.modLong) is the Euclidean % for every integer. The first two stages of the sieve reduce the sums of a word number and the offset of a row modulo the moduli of their sieve specs, by a multiplication with the reciprocal when the test small says that the word numbers of the piece allow it and by the helper mod otherwise; for a sieve spec whose offset lies in the range that the packing gives it (Ratpoints.Sieve.SieveSpec.Fits), both give the Euclidean remainder.

Main definitions

• Ratpoints.Sieve.SieveSpec.Wraps: the row of a sieve spec repeats its first bit arrays after the first p, where the chunked first stage reads it.
• Ratpoints.Sieve.SieveSpec.Fits: a sieve spec of a modulus 0 < p ≤ Ratpoints.maxPrimeEven with an offset from Ratpoints.rowBias to below Ratpoints.rowBias + 2 p, the end p and a row that repeats its first bit arrays.
• Ratpoints.Sieve.SmallFits: when the test small holds, the word numbers lie between -Ratpoints.rowBias and Ratpoints.rowBias - 2 Ratpoints.maxPrimeEven.
• Ratpoints.Sieve.StageFits: what the stages of the sieve need of their moduli, the reciprocals and the test small.

Main statements

• Ratpoints.maxPrimeEven_le_two_pow_10, Ratpoints.maxPrimeEven_le_two_pow_32: the moduli are below 2 ^ 10, and so below 2 ^ 32; Ratpoints.maxPrimeEven_le_sixty_four_mul_modWords: they fit the bit arrays of Ratpoints.modWords words.
• Ratpoints.rowBias_le_two_pow_31: Ratpoints.rowBias, the lower bound of the bias of a modulus (the multiple of the modulus in an offset), is at most 2 ^ 31.
• Ratpoints.mul_lt_two_pow_64_of_le_mulModLimit: u p < 2 ^ 64 for u ≤ Ratpoints.mulModLimit and p ≤ 2 ^ 32.
• Ratpoints.mul_le_mulModLimit_of_le_div, Ratpoints.mul_lt_two_pow_63: the products of fill_bp_list and of class_offsets are longs, those of the first branch of fill_bp_list at most Ratpoints.mulModLimit.
• Ratpoints.rbaLength_eq, Ratpoints.rbaLength_pos, Ratpoints.natCast_rbaLength_pos, Ratpoints.natCast_rbaLength_dvd_two_pow: the width Ratpoints.rbaLength of the bit arrays is 64 Ratpoints.rbaPack, positive, and divides 2 ^ n for n ≥ Ratpoints.rbaShift.
• Ratpoints.Sieve.toNat_mulHi: Ratpoints.Sieve.mulHi a b is ⌊a b / 2 ^ 64⌋; Ratpoints.Sieve.toNat_hi_mul_add_lo, Ratpoints.Sieve.toNat_ite_bne_zero: the 128-bit sum a m + r of exact_scale in its two words, and the quotient on two words capped at 2 ^ 64 - 1.
• Ratpoints.Sieve.mulMod_eq, Ratpoints.Sieve.mulDiv_eq: the reductions on 64-bit words are those of the macros on natural numbers reduced modulo 2 ^ 64.
• Ratpoints.Sieve.magic_mul_eq: magic p · p = 2 ^ 64 + δ with δ < p.
• Ratpoints.Sieve.mulMod_magic, Ratpoints.Sieve.mulDiv_magic: with the reciprocal Ratpoints.Sieve.magic p, the reductions give u % p and, for p ≥ 2, u / p for u p < 2 ^ 64; Ratpoints.Sieve.mulMod_magic_of_le_mulModLimit: the remainder of the values at most Ratpoints.mulModLimit modulo p ≤ 2 ^ 32.
• Ratpoints.Sieve.modMul_magic: mod_mul gives the Euclidean a % p for |a| p < 2 ^ 64.
• Ratpoints.Sieve.modLong_eq_emod: the helper mod of the C code is the Euclidean %.
• Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_eq_emod_of_natAbs_le: the reduction by mod_mul under the test small and by the helper mod otherwise is the Euclidean remainder of a value at most Ratpoints.mulModLimit in absolute value under small.
• Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_eq_emod, Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_add_eq_emod, Ratpoints.Sieve.SieveSpec.Fits.ite_mulMod_modLong_eq_toNat_emod: the first row index of the first stage one row at a time, the start of a row of the chunked first stage and the row index of the second stage are the Euclidean remainders of -w0 - offset, of w0 + offset and of i + offset.

Implementation notes

The file uses Lean core and Std only, so that the contracts of the sieve, in files without Mathlib, can use it. The proofs use only the facts on the constants that decide checks next to them, that Ratpoints.rbaPack > 0, that Ratpoints.mulModLimit < 2 ^ 32, that Ratpoints.maxBitsInPrime ≤ 10 and that 2 Ratpoints.rowBias - 1 ≤ Ratpoints.mulModLimit, never their values; Ratpoints.rowBias is generalized before lia, which otherwise fails on its value with a deep recursion of the kernel. Ratpoints.mul_lt_two_pow_63 and Ratpoints.rowBias_le_two_pow_31 (through Ratpoints.Sieve.SieveEntry.PackedFor.lt_two_pow_32) discharge range obligations of the longs that the Deviation paragraphs of the model cite, and no contract uses them.

The moduli

theorem Ratpoints.maxPrimeEven_le_two_pow_10 : maxPrimeEven ≤ 2 ^ 10

The moduli below Ratpoints.maxPrimeEven are below 2 ^ 10, since maxBitsInPrime ≤ 10 (Ratpoints.maxBitsInPrime_le_ten), as the comment on Horner's rule in examine_power assumes (find_points.c:2286-2287).

theorem Ratpoints.maxPrimeEven_le_two_pow_32 : maxPrimeEven ≤ 2 ^ 32

The moduli, below Ratpoints.maxPrimeEven, are below 2 ^ 32, the bound under which the reductions by a reciprocal of the values at most Ratpoints.mulModLimit are exact (Ratpoints.mul_lt_two_pow_64_of_le_mulModLimit) and the helper mod keeps its values longs (Ratpoints.fitsLong_sixteen_mul, Ratpoints.fitsLong_add_sixteen_mul).

theorem Ratpoints.maxPrimeEven_le_sixty_four_mul_modWords : maxPrimeEven ≤ 64 * modWords

A prime power below Ratpoints.maxPrimeEven, the moduli of the parameters, fits the Ratpoints.modWords words of the bit arrays of its tables and of the patterns of its rows.

theorem Ratpoints.rowBias_le_two_pow_31 : rowBias ≤ 2 ^ 31

Ratpoints.rowBias is at most 2 ^ 31, since 2 rowBias - 1 ≤ Ratpoints.mulModLimit < 2 ^ 32 (Ratpoints.two_mul_rowBias_sub_one_le_mulModLimit, Ratpoints.mulModLimit_lt_two_pow_32): the bound of the offsets of the rows (Ratpoints.Sieve.SieveEntry.PackedFor.lt_two_pow_32).

theorem Ratpoints.mul_lt_two_pow_64_of_le_mulModLimit {u p : Nat} (hu : u ≤ mulModLimit) (hp : p ≤ 2 ^ 32) : u * p < 2 ^ 64

A value at most Ratpoints.mulModLimit times a modulus p ≤ 2 ^ 32 is below 2 ^ 64, where the reductions by a reciprocal are exact (Ratpoints.Sieve.mulMod_magic, Ratpoints.Sieve.mulDiv_magic); this reads off only Ratpoints.mulModLimit_lt_two_pow_32.

theorem Ratpoints.mul_mul_lt_two_pow_64 {x d p : Nat} (hx : x ≤ 2 ^ 10) (hd : d < p) (hp : p ≤ maxPrimeEven) : x * d * p < 2 ^ 64

The products x d p of fill_bp_list and class_offsets that a reduction by the reciprocal of a modulus p ≤ Ratpoints.maxPrimeEven takes, for x ≤ 2 ^ 10 and d < p, are below 2 ^ 64, where it is exact.

theorem Ratpoints.mul_le_mulModLimit_of_le_div {b d p : Nat} (hb : b ≤ mulModLimit / maxPrimeEven) (hd : d < p) (hp : p ≤ maxPrimeEven) : b * d ≤ mulModLimit

The products b dinv of the first branch of fill_bp_list, for b ≤ Ratpoints.mulModLimit / Ratpoints.maxPrimeEven and dinv < p ≤ Ratpoints.maxPrimeEven, are at most Ratpoints.mulModLimit, so longs that the reduction by the reciprocal of p reduces exactly.

theorem Ratpoints.mul_lt_two_pow_63 {x d p : Nat} (hx : x ≤ 2 ^ 10) (hd : d < p) (hp : p ≤ maxPrimeEven) : x * d < 2 ^ 63

The products x d of the other branches of fill_bp_list, (b mod p) dinv, and of class_offsets, a₀ rinv, for x ≤ 2 ^ 10 and d < p ≤ Ratpoints.maxPrimeEven, are longs.

theorem Ratpoints.mod_le_two_pow_10 (b : Nat) {p : Nat} (hp : 0 < p) (hpm : p ≤ maxPrimeEven) : b % p ≤ 2 ^ 10

The width of the bit arrays

theorem Ratpoints.rbaLength_eq : rbaLength = 64 * rbaPack

theorem Ratpoints.rbaLength_pos : 0 < rbaLength

theorem Ratpoints.natCast_rbaLength_pos : 0 < ↑rbaLength

theorem Ratpoints.natCast_rbaLength_dvd_two_pow {n : Nat} (hn : rbaShift ≤ n) : ↑rbaLength ∣ 2 ^ n

The width Ratpoints.rbaLength = 2 ^ Ratpoints.rbaShift of the bit arrays divides 2 ^ n for n ≥ Ratpoints.rbaShift.

The high half of a product

theorem Ratpoints.Sieve.toNat_mulHi (a b : UInt64) : (mulHi a b).toNat = a.toNat * b.toNat / 2 ^ 64

The high half of the product of two 64-bit words from the products of their 32-bit halves is the high half of the product of natural numbers: no partial sum wraps.

theorem Ratpoints.Sieve.toNat_hi_mul_add_lo {a m rl rh : UInt64} (h : a.toNat * m.toNat + (rh.toNat * 2 ^ 64 + rl.toNat) < 2 ^ 128) : (mulHi a m + rh + if a * m + rl < rl then 1 else 0).toNat * 2 ^ 64 + (a * m + rl).toNat = a.toNat * m.toNat + (rh.toNat * 2 ^ 64 + rl.toNat)

The 128-bit sum b m + r of Ratpoints.Sieve.exactScale in its two words.

theorem Ratpoints.Sieve.toNat_ite_bne_zero {qh ql : UInt64} {N : Nat} (h : qh.toNat * 2 ^ 64 + ql.toNat = N) : (if (qh != 0) = true then ~~~0 else ql).toNat = min (2 ^ 64 - 1) N

The result of Ratpoints.Sieve.exactScale from the two words of its quotient: 2 ^ 64 - 1 when the high word is not 0.

The reductions by a reciprocal

theorem Ratpoints.Sieve.mulMod_eq (u p m : Nat) : mulMod u p m = m * u % 2 ^ 64 * (p % 2 ^ 64) / 2 ^ 64

RP_MULMOD on 64-bit words is the macro on natural numbers: the product m u reduced modulo 2 ^ 64, as the product of unsigned longs wraps, times p as an unsigned long, shifted right by 64 bits.

theorem Ratpoints.Sieve.mulDiv_eq (u m : Nat) : mulDiv u m = u % 2 ^ 64 * (m % 2 ^ 64) / 2 ^ 64

theorem Ratpoints.Sieve.magic_mul_eq {p : Nat} (hp : 0 < p) : magic p * p = 2 ^ 64 + (p - 1 - (2 ^ 64 - 1) % p)

The reciprocal of p > 0 exceeds 2 ^ 64 / p by less than 1: magic p · p = 2 ^ 64 + δ with δ < p.

theorem Ratpoints.Sieve.magic_lt_two_pow_64 {p : Nat} (hp : 1 < p) : magic p < 2 ^ 64

The reciprocal of a modulus p ≥ 2 is an unsigned long; that of 1, no modulus, is 2 ^ 64, which wraps to 0 in the C code.

theorem Ratpoints.Sieve.mulMod_magic {u p : Nat} (hp : 0 < p) (hu : u * p < 2 ^ 64) : mulMod u p (magic p) = u % p

RP_MULMOD with the reciprocal of p computes the remainder of every u with u p < 2 ^ 64, since magic p exceeds 2 ^ 64 / p by less than 1.

theorem Ratpoints.Sieve.mulDiv_magic {u p : Nat} (hp : 1 < p) (hu : u * p < 2 ^ 64) : mulDiv u (magic p) = u / p

RP_MULDIV with the reciprocal of p ≥ 2 computes the quotient of every u with u p < 2 ^ 64. For p = 1 the reciprocal 2 ^ 64 wraps to 0, as in the C code.

theorem Ratpoints.Sieve.mulMod_magic_of_le_mulModLimit {u p : Nat} (hp : 0 < p) (hp32 : p ≤ 2 ^ 32) (hu : u ≤ mulModLimit) : mulMod u p (magic p) = u % p

RP_MULMOD with the reciprocal of a modulus 0 < p ≤ 2 ^ 32 computes the remainder of every u ≤ Ratpoints.mulModLimit (Ratpoints.mul_lt_two_pow_64_of_le_mulModLimit).

theorem Ratpoints.Sieve.modMul_magic {a : Int} {p : Nat} (hp : 0 < p) (ha : a.natAbs * p < 2 ^ 64) : modMul a p (magic p) = a % ↑p

mod_mul with the reciprocal of p computes the Euclidean remainder of every a with |a| p < 2 ^ 64: the remainder of |a|, and its complement for a negative a.

The helper mod

theorem Ratpoints.Sieve.modLong_eq_emod (a : Int) {b : Int} (hb : 0 ≤ b) : modLong a b = a % b

The helper mod of the C code is the Euclidean %, for every integer a and every b ≥ 0; that its values stay longs is Ratpoints.fitsLong_sixteen_mul, Ratpoints.fitsLong_neg and Ratpoints.fitsLong_add_sixteen_mul.

The row indices of the first two stages

def Ratpoints.Sieve.SieveSpec.Wraps (s : SieveSpec) : Prop

The row of the sieve spec s repeats its first p bit arrays in the Ratpoints.chunk - 1 bit arrays after them: word κ < Ratpoints.rbaPack of bit array x < p + chunk - 1 is word κ of bit array x mod p, the copies that the chunked first stage reads past the first p bit arrays from a start below p (Ratpoints.Sieve.andRowsChunked).

Equations

• s.Wraps = ∀ (x : Nat), x < s.p + Ratpoints.chunk - 1 → ∀ (κ : Nat), κ < Ratpoints.rbaPack → s.row[Ratpoints.rbaPack * x + κ]! = s.row[Ratpoints.rbaPack * (x % s.p) + κ]!

structure Ratpoints.Sieve.SieveSpec.Fits (s : SieveSpec) : Prop

A sieve spec whose modulus p and offset lie in the range that the packing of the numerators gives them (Ratpoints.Sieve.classOffsets) and whose row has the copies of its first bit arrays with the end p, as the run sets them up (Ratpoints.Sieve.setRows): 0 < p ≤ Ratpoints.maxPrimeEven, the offset a residue below 2 p plus Ratpoints.rowBias, so that under the test small the sums of a word number and the offset that the first two stages reduce are at most Ratpoints.mulModLimit in absolute value (Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_eq_emod, Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_add_eq_emod, Ratpoints.Sieve.SieveSpec.Fits.ite_mulMod_modLong_eq_toNat_emod), and the row repeating its first bit arrays (Ratpoints.Sieve.SieveSpec.Wraps) where the chunked first stage reads it.

• pos : 0 < s.p

  The modulus is positive.

• le_maxPrimeEven : s.p ≤ maxPrimeEven

  The modulus is at most Ratpoints.maxPrimeEven.

• rowBias_le : rowBias ≤ s.offset

  The offset is at least Ratpoints.rowBias.

• lt_rowBias_add : s.offset < rowBias + 2 * s.p

  The offset is below Ratpoints.rowBias + 2 p.

• stop_eq : s.stop = s.p

  The end of the walk of the chunked first stage through the row is the modulus.

• wraps : s.Wraps

  The row repeats its first bit arrays after the first p.

def Ratpoints.Sieve.SmallFits (small : Bool) (lo hi : Int) : Prop

When the test small of Ratpoints.Sieve.siftPiece holds, the word numbers from lo to hi lie between -Ratpoints.rowBias and Ratpoints.rowBias - 2 Ratpoints.maxPrimeEven: the hypothesis under which the first two stages may reduce by a multiplication.

Equations

• Ratpoints.Sieve.SmallFits small lo hi = (small = true → -↑Ratpoints.rowBias ≤ lo ∧ hi ≤ ↑Ratpoints.rowBias - 2 * ↑Ratpoints.maxPrimeEven)

theorem Ratpoints.Sieve.smallFits_decide (w0 : Int) (n : Nat) : SmallFits (decide (-↑rowBias ≤ w0 ∧ w0 + ↑n ≤ ↑rowBias - 2 * ↑maxPrimeEven)) w0 (w0 + ↑n)

The test small of a piece, as Ratpoints.Sieve.siftPiece computes it, says what Ratpoints.Sieve.SmallFits asks of its word numbers.

theorem Ratpoints.Sieve.SmallFits.mono {small : Bool} {lo hi lo' hi' : Int} (h : SmallFits small lo hi) (hlo : lo ≤ lo') (hhi : hi' ≤ hi) : SmallFits small lo' hi'

structure Ratpoints.Sieve.StageFits (specs : List SieveSpec) (mgs : List Nat) (small : Bool) (lo hi : Int) : Prop

What the stages of the sieve need of the moduli specs with their reciprocals mgs, for the word numbers from lo to hi: the sieve specs in the range of the packing, the reciprocals of their moduli in the same order, and the test small of Ratpoints.Sieve.SmallFits.

• fits (s : SieveSpec) : s ∈ specs → s.Fits

  The sieve specs are in the range of the packing.

• magics : mgs = List.map (fun (x : SieveSpec) => magic x.p) specs

  The reciprocals are those of the moduli, in the same order.

• smallFits : SmallFits small lo hi

  The test small fits the word numbers.

theorem Ratpoints.Sieve.StageFits.mono {specs : List SieveSpec} {mgs : List Nat} {small : Bool} {lo hi lo' hi' : Int} (h : StageFits specs mgs small lo hi) (hlo : lo ≤ lo') (hhi : hi' ≤ hi) : StageFits specs mgs small lo' hi'

theorem Ratpoints.Sieve.stageFits_take {specs : List SieveSpec} {mgs : List Nat} {small : Bool} {lo hi : Int} (h : (∀ (s : SieveSpec), s ∈ specs → s.Fits) ∧ mgs = List.map (fun (x : SieveSpec) => magic x.p) specs) (hs : SmallFits small lo hi) (k : Nat) : StageFits (List.take k specs) (List.take k mgs) small lo hi

The first k of the moduli, with their reciprocals.

theorem Ratpoints.Sieve.stageFits_drop {specs : List SieveSpec} {mgs : List Nat} {small : Bool} {lo hi : Int} (h : (∀ (s : SieveSpec), s ∈ specs → s.Fits) ∧ mgs = List.map (fun (x : SieveSpec) => magic x.p) specs) (hs : SmallFits small lo hi) (k : Nat) : StageFits (List.drop k specs) (List.drop k mgs) small lo hi

The moduli after the first k, with their reciprocals.

theorem Ratpoints.Sieve.SieveSpec.Fits.natAbs_neg_sub_le {s : SieveSpec} (hs : s.Fits) {w0 : Int} (h0 : -↑rowBias ≤ w0) (h1 : w0 ≤ ↑rowBias - 2 * ↑maxPrimeEven) : (-w0 - ↑s.offset).natAbs ≤ mulModLimit

For a word number -Ratpoints.rowBias ≤ w0 ≤ Ratpoints.rowBias - 2 Ratpoints.maxPrimeEven, the value -w0 - offset that the first stage reduces is at most Ratpoints.mulModLimit in absolute value.

theorem Ratpoints.Sieve.SieveSpec.Fits.add_nonneg_and_toNat_le {s : SieveSpec} (hs : s.Fits) {i : Int} (h0 : -↑rowBias ≤ i) (h1 : i < ↑rowBias - 2 * ↑maxPrimeEven) : 0 ≤ i + ↑s.offset ∧ (i + ↑s.offset).toNat ≤ mulModLimit

For a word number -Ratpoints.rowBias ≤ i < Ratpoints.rowBias - 2 Ratpoints.maxPrimeEven, the value i + offset that the second stage reduces lies between 0 and Ratpoints.mulModLimit.

theorem Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_eq_emod_of_natAbs_le {s : SieveSpec} (hs : s.Fits) {small : Bool} {a : Int} (h : small = true → a.natAbs ≤ mulModLimit) : (if small = true then modMul a s.p (magic s.p) else modLong a ↑s.p) = a % ↑s.p

The reduction of a modulo the modulus of a sieve spec in the range of the packing, by mod_mul under the test small and by the helper mod otherwise, is the Euclidean remainder, for |a| ≤ Ratpoints.mulModLimit under small.

theorem Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_eq_emod {s : SieveSpec} (hs : s.Fits) {small : Bool} {w0 : Int} {n : Nat} (h : SmallFits small w0 (w0 + ↑n)) : (if small = true then modMul (-w0 - ↑s.offset) s.p (magic s.p) else modLong (-w0 - ↑s.offset) ↑s.p) = (-w0 - ↑s.offset) % ↑s.p

The index of the first row of the first stage over a piece from the word number w0 on, by mod_mul under the test small and by the helper mod otherwise, is the Euclidean remainder of -w0 - offset.

theorem Ratpoints.Sieve.SieveSpec.Fits.ite_modMul_modLong_add_eq_emod {s : SieveSpec} (hs : s.Fits) {small : Bool} {w0 : Int} (h : SmallFits small w0 w0) : (if small = true then modMul (w0 + ↑s.offset) s.p (magic s.p) else modLong (w0 + ↑s.offset) ↑s.p) = (w0 + ↑s.offset) % ↑s.p

The start of the row of the chunked first stage at the first word number w0 of a piece, by mod_mul under the test small and by the helper mod otherwise, is the Euclidean remainder of w0 + offset.

theorem Ratpoints.Sieve.SieveSpec.Fits.ite_mulMod_modLong_eq_toNat_emod {s : SieveSpec} (hs : s.Fits) {small : Bool} {i : Int} (h : SmallFits small i (i + 1)) : (if small = true then mulMod (i + ↑s.offset).toNat s.p (magic s.p) else (modLong (i + ↑s.offset) ↑s.p).toNat) = ((i + ↑s.offset) % ↑s.p).toNat

The row index of the second stage at the word number i, by a multiplication with the reciprocal under the test small and by the helper mod otherwise, is the Euclidean remainder of i + offset.