src/math/w64.bend checks

a + b >= 2^64

a - b mod 2^64

the full product of two U32

a * b mod 2^64

a * b >= 2^64

One digit of the long division by d: (r * base + g) / d and mod d, for a remainder r < d < 2^32 and a digit g < base = 2^16, so every partial dividend r * 2^16 + g is below 2^48 (the runtime's Nat bound). The low word is taken as two 16-bit digits; the base is written 256 * 256 so the proof checker never expands a 17-bit literal.

(a / d, a mod d) for d > 0: the high word natively, the low word by digits

a mod d for d > 0

2^k for k < 32

the number of bits of x (0 for 0), by a binary search on the top bit

a << k mod 2^64 for k < 64

a >> k (0 for k >= 64)

min(q, 2^32 - 1)

the product of a U32 and a U64 as (low 64 bits, top 32 bits)

the 96-bit product p exceeds x = xh * 2^64 + xl

q - 1 while q * b > x: at most q steps (fuel q + 1)

floor(x / b) from an estimate q >= floor(x / b): down while q * b > x (at most q steps, fuel q + 1)

Knuth's estimate for x = xh * 2^64 + xl < b * 2^32 and b >= 2^32, t = bits(hi(b)): with Y = x >> t (below 2^64) and B = b >> t in [2^31, 2^32), min(floor(Y / B), 2^32 - 1) is at least floor(x / b) (as B 2^t <= b) and at most 2 over it (Theorem B), so q_start takes a step or two

min(floor(y / d), 2^32 - 1) for d > 0: when hi(y) >= d the quotient needs more than 32 bits; otherwise the high word is its own remainder, and div32's two 16-bit Nat digit steps finish the division without dividing the high word.

floor(x / b) for x = xh * 2^64 + xl < b * 2^32 and b >= 2^32, t = bits(hi(b))

(a / b, a mod b) for b > 0

x mod m for x = xh * 2^32 + x0 with xh < m and m >= 2^32

(the zero test of hi(m) only guards the digit: red96 is called with hi(m) != 0; it also stops the proof checker from expanding the digit in every statement that names red96)

the full product as (low 64 bits, high 64 bits)

a * b mod m for a, b < m

r - 1 while r * r > x: from r <= 65535 at most r steps (fuel r + 1), and r * r never wraps

r + 1 while r < 65535 and (r + 1)^2 <= x: at most 65535 - r steps

floor(sqrt(x)) from any estimate r <= 65535: down while r^2 > x, then up while (r + 1)^2 <= x. The F32 square root lands within 1, so both loops take at most a step or two; the fuel covers any estimate, so the result never depends on F32's accuracy.

floor(sqrt(n)) from any estimate r1: down while r^2 > n, then up while (r + 1)^2 <= n; the fuel covers any estimate

one Newton step from the F32 estimate r0 > 0 lands within a few units of floor(sqrt(n)); the correction makes the result exact for any r0

a >> k with the bits shifted out OR-ed into bit 0 (SoftFloat's shiftRightJam64)

the number of leading zero bits of a (64 for 0)

a with bit 0 cleared when b

x (2 - m x): doubles the correct low bits of an inverse of m mod 2^64

-1 / m mod 2^64 for odd m (m is its own inverse mod 8; five Newton steps reach 96 >= 64 bits)

t / 2^64 mod m for t = (tl, th) < m^2

a b / 2^64 mod m for a, b < m

pow_mod_go of generic.bend on Montgomery forms x 2^64 mod m

x 2^64 mod m for x < m, m >= 2^32

m odd with 2^32 <= m < 2^63 and m minv(m) = -1 mod 2^64 (always so for odd m; checked, so the proof needs no Hensel lemma)

b^e mod m for b < m and mont_ok(m)

-1 / m mod 2^32 for odd m

a and b fit 48 bits (runtime Nats) and a has more than 16: the binary gcd then hands them to Euclid on Nat (native division); below 2^16 its few steps are cheaper

the words of g < 2^64 from h = g div 2^32 (two divisions by 2^16, so no closed 2^32 appears: the proof checker would expand it in unary)