Lean proofs for exact quotient and remainder recovery using an upper magic number

README.ja.md

ExactDivRem.lean

上側 magic number を使って、商と剰余を正確に復元する条件を Lean で形式化したファイルです。

このファイルでは、次の上側近似

$$ M = \left\lceil \frac{\beta}{D} \right\rceil $$

を用いて乗算により商を近似し、積 $Mx$ から商と剰余の両方を復元するための正しさ条件を証明しています。主な用途は Barrett 型除算や magic number による除算で、特に $\beta$ が 2 の冪である場合を想定しています。

概要

次のようにおきます。

$$ q := \left\lfloor \frac{x}{D} \right\rfloor, \qquad r := x \bmod D, \qquad M := \left\lceil \frac{\beta}{D} \right\rceil, \qquad \delta := MD - \beta. $$

$\delta$ は、 $\beta / D$ を上側整数近似 $M$ に置き換えたことで生じる余剰を表します。

このファイルでは、補助残差

$$ Y := q\delta + Mr $$

を導入します。この量は次の 2 つの恒等式を満たすため有用です。

$$ Mx = q\beta + Y $$

および

$$ DY = x\delta + \beta r. $$

これらの恒等式により、商と剰余の正しさは $x\delta$ に関する不等式へ帰着されます。

approxQuot_eq_div_iff_mul_delta_lt:

$$ x\delta < \beta(D-r) \iff \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor. $$

approxDivRem_iff_mul_delta_lt:

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). $$

approxDivRem_pow2_of_delta_le:

$$ \begin{aligned} & \beta = 2^{N+L}, \quad D < 2^N, \quad x < 2^N, \quad \delta \le 2^L \\ \implies & \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). \end{aligned} $$

主な定義

magicUp

def magicUp (β D : Nat) : Nat :=
  β ⌈/⌉ D

上側 magic number

$$ M = \left\lceil \frac{\beta}{D} \right\rceil $$

を定義します。

delta

def delta (β D : Nat) : Nat :=
  magicUp β D * D - β

余剰

$$ \delta = MD - \beta $$

を定義します。

approxQuot

def approxQuot (β D x : Nat) : Nat :=
  (magicUp β D * x) / β

近似商

$$ \left\lfloor \frac{Mx}{\beta} \right\rfloor $$

を定義します。

Good

def Good (β D x : Nat) : Prop :=
  x / D = approxQuot β D x

近似商が真の商と一致することを表します。

upperY

def upperY (β D x : Nat) : Nat :=
  (x / D) * delta β D + magicUp β D * (x % D)

残差

$$ Y = q\delta + Mr $$

を定義します。これは証明全体で中心となる量です。

主定理

approxQuot_eq_div_iff_mul_delta_lt

theorem approxQuot_eq_div_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    x * δ < β * (D - r) ↔ (M * x) / β = q

この定理は、商を正確に復元できるための必要十分条件を与えます。

数学的には、

$$ x\delta < \beta(D-r) \iff \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor $$

です。つまり、上側 magic number による商が正確であることは、累積誤差 $x\delta$ が利用可能な余裕 $\beta(D-r)$ より小さいことと同値です。

approxDivRem_iff_mul_delta_lt

theorem approxDivRem_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    let y := (M * x) % β
    x * δ < β ↔ ((M * x) / β = q ∧ ((y * D) / β = r))

この定理は、商と剰余を同時に正確に復元できる条件を特徴付けます。

$$ y := Mx \bmod \beta $$

とおくと、定理は次の同値を主張します。

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). $$

同じことを $y := Mx \bmod \beta$ を使って書くと、

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = q \land \left\lfloor \frac{yD}{\beta} \right\rfloor = r \right) $$

です。言い換えると、この条件がちょうど成り立つとき、 $Mx$ の上位部分が商を、 $Mx$ の下位部分をスケールしたものが剰余を与えます。

approxDivRem_pow2_of_delta_le

theorem approxDivRem_pow2_of_delta_le
    {N L D x : Nat} (hD : 0 < D) (hDlt : D < 2 ^ N)
    (hδ : ((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L) ≤ 2 ^ L)
    (hx : x < 2 ^ N) :
    let β := 2 ^ (N + L)
    let M := β ⌈/⌉ D
    let y := (M * x) % β
    (M * x) / β = x / D ∧ (y * D) / β = x % D

これは機械整数実装で有用な、2 の冪基数に対する系です。

$$ \beta = 2^{N+L}, \qquad M = \left\lceil \frac{2^{N+L}}{D} \right\rceil, \qquad y = Mx \bmod \beta $$

とおきます。結論は、商と剰余を同時に復元する次の式です。

$$ \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{yD}{\beta} \right\rfloor = x \bmod D. $$

$\beta = 2^{N+L}$ なので、同値に次のようにも書けます。

$$ \left\lfloor \frac{Mx}{2^{N+L}} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod 2^{N+L})D}{2^{N+L}} \right\rfloor = x \bmod D. $$

もし

$$ D < 2^N, \qquad x < 2^N, \qquad \delta \le 2^L $$

ならば、

$$ x\delta \le (2^N - 1)2^L < 2^{N+L} = \beta $$

が成り立ちます。したがって approxDivRem_iff_mul_delta_lt を適用でき、上に表示した商と剰余の式がどちらも成立します。

$\beta = 2^{N+L}$ の場合、これは次の実装パターンに対応します。

t = M * x
q = t >> (N + L)
y = t mod 2^(N + L)
r = (y * D) >> (N + L)

このとき、上記の仮定のもとで

q = x / D
r = x % D

が成り立ちます。

証明構造

証明は残差項 upperY を中心に構成されています。

主要な内部恒等式は次の 2 つです。

magicUp β D * x = (x / D) * β + upperY β D x

および

D * upperY β D x = x * delta β D + β * (x % D)

数学的には、これらは

$$ Mx = q\beta + Y $$

および

$$ DY = x\delta + \beta r $$

です。

最初の恒等式から、近似商が正確であることは

$$ Y < \beta $$

と同値です。

2 つ目の恒等式から、

$$ Y < \beta \iff x\delta < \beta(D-r) $$

が得られます。

商と剰余を同時に復元する場合は、より強い条件

$$ x\delta < \beta $$

により

$$ Mx \bmod \beta = Y $$

が保証されます。すると

$$ \left\lfloor \frac{YD}{\beta} \right\rfloor = r $$

が成り立ちます。

逆向きは同じ恒等式を逆にたどることで得られます。

$$ \left\lfloor \frac{YD}{\beta} \right\rfloor = r + \left\lfloor \frac{x\delta}{\beta} \right\rfloor. $$

したがって、剰余が正確に復元されるなら

$$ \left\lfloor \frac{x\delta}{\beta} \right\rfloor = 0 $$

が強制され、これは

$$ x\delta < \beta $$

と同値です。

直接剰余計算との関係

2 の冪基数に対する系は、事前計算した逆数を用いた直接剰余計算と密接に関係します。 対応する文献は次の通りです。

Daniel Lemire, Owen Kaser, Nathan Kurz,
Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries, arXiv:1902.01961, 2019.
https://arxiv.org/abs/1902.01961

Lean の定理は、商の復元と剰余の復元の両方を同時に扱う形を証明しています。

要件

このファイルは次を import します。

import Mathlib.Algebra.Order.Floor.Div
import Mathlib.Tactic

Lean 4 と Mathlib を想定しています。

README.md

ExactDivRem.lean

Lean formalization of exact quotient and remainder recovery using an upper magic number.

This file proves correctness conditions for computing a quotient approximation by multiplication with

$$ M = \left\lceil \frac{\beta}{D} \right\rceil $$

and for recovering both the quotient and the remainder from the product $Mx$. The main use case is a Barrett-style or magic-number division scheme, especially when $\beta$ is a power of two.

Overview

Let

$$ q := \left\lfloor \frac{x}{D} \right\rfloor, \qquad r := x \bmod D, \qquad M := \left\lceil \frac{\beta}{D} \right\rceil, \qquad \delta := MD - \beta. $$

The number $\delta$ measures the excess introduced by replacing $\beta / D$ with the upper integer approximation $M$.

The file introduces the auxiliary residual

$$ Y := q\delta + Mr. $$

This residual is useful because it satisfies the two identities

$$ Mx = q\beta + Y $$

and

$$ DY = x\delta + \beta r. $$

These identities reduce quotient and remainder correctness to inequalities involving $x\delta$.

approxQuot_eq_div_iff_mul_delta_lt:

$$ x\delta < \beta(D-r) \iff \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor. $$

approxDivRem_iff_mul_delta_lt:

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). $$

approxDivRem_pow2_of_delta_le:

$$ \begin{aligned} & \beta = 2^{N+L}, \quad D < 2^N, \quad x < 2^N, \quad \delta \le 2^L \\ \implies & \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). \end{aligned} $$

Main definitions

magicUp

def magicUp (β D : Nat) : Nat :=
  β ⌈/⌉ D

Defines the upper magic number

$$ M = \left\lceil \frac{\beta}{D} \right\rceil. $$

delta

def delta (β D : Nat) : Nat :=
  magicUp β D * D - β

Defines the excess

$$ \delta = MD - \beta. $$

approxQuot

def approxQuot (β D x : Nat) : Nat :=
  (magicUp β D * x) / β

Defines the approximate quotient

$$ \left\lfloor \frac{Mx}{\beta} \right\rfloor. $$

Good

def Good (β D x : Nat) : Prop :=
  x / D = approxQuot β D x

States that the approximate quotient agrees with the true quotient.

upperY

def upperY (β D x : Nat) : Nat :=
  (x / D) * delta β D + magicUp β D * (x % D)

Defines the residual

$$ Y = q\delta + Mr. $$

This is the central quantity used in the proofs.

Main theorems

approxQuot_eq_div_iff_mul_delta_lt

theorem approxQuot_eq_div_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    x * δ < β * (D - r) ↔ (M * x) / β = q

This theorem gives an exact necessary and sufficient condition for quotient recovery.

In mathematical notation,

$$ x\delta < \beta(D-r) \iff \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor. $$

Thus the upper-magic quotient is exact precisely when the accumulated error $x\delta$ is smaller than the available margin $\beta(D-r)$.

approxDivRem_iff_mul_delta_lt

theorem approxDivRem_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    let y := (M * x) % β
    x * δ < β ↔ ((M * x) / β = q ∧ ((y * D) / β = r))

This theorem characterizes exact simultaneous quotient and remainder recovery.

Let

$$ y := Mx \bmod \beta. $$

The theorem states the following equivalence:

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod \beta)D}{\beta} \right\rfloor = x \bmod D \right). $$

Equivalently, using $y := Mx \bmod \beta$,

$$ x\delta < \beta \iff \left( \left\lfloor \frac{Mx}{\beta} \right\rfloor = q \land \left\lfloor \frac{yD}{\beta} \right\rfloor = r \right). $$

In other words, under exactly this condition, the high part of $Mx$ gives the quotient and the scaled low part of $Mx$ gives the remainder.

approxDivRem_pow2_of_delta_le

theorem approxDivRem_pow2_of_delta_le
    {N L D x : Nat} (hD : 0 < D) (hDlt : D < 2 ^ N)
    (hδ : ((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L) ≤ 2 ^ L)
    (hx : x < 2 ^ N) :
    let β := 2 ^ (N + L)
    let M := β ⌈/⌉ D
    let y := (M * x) % β
    (M * x) / β = x / D ∧ (y * D) / β = x % D

This is a power-of-two corollary useful for machine-integer implementations.

Set

$$ \beta = 2^{N+L}, \qquad M = \left\lceil \frac{2^{N+L}}{D} \right\rceil, \qquad y = Mx \bmod \beta. $$

The conclusion is the following simultaneous quotient and remainder recovery formula:

$$ \left\lfloor \frac{Mx}{\beta} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{yD}{\beta} \right\rfloor = x \bmod D. $$

Equivalently, since $\beta = 2^{N+L}$,

$$ \left\lfloor \frac{Mx}{2^{N+L}} \right\rfloor = \left\lfloor \frac{x}{D} \right\rfloor \land \left\lfloor \frac{(Mx \bmod 2^{N+L})D}{2^{N+L}} \right\rfloor = x \bmod D. $$

If

$$ D < 2^N, \qquad x < 2^N, \qquad \delta \le 2^L, $$

then

$$ x\delta \le (2^N - 1)2^L < 2^{N+L} = \beta. $$

Therefore approxDivRem_iff_mul_delta_lt applies, and the displayed quotient and remainder formulas both hold.

For $\beta = 2^{N+L}$, this corresponds to the implementation pattern

t = M * x
q = t >> (N + L)
y = t mod 2^(N + L)
r = (y * D) >> (N + L)

with

q = x / D
r = x % D

under the stated assumptions.

Proof structure

The proof is organized around the residual term upperY.

The key internal identities are:

magicUp β D * x = (x / D) * β + upperY β D x

and

D * upperY β D x = x * delta β D + β * (x % D)

In mathematical notation, these are

$$ Mx = q\beta + Y $$

and

$$ DY = x\delta + \beta r. $$

From the first identity, the approximate quotient is exact precisely when

$$ Y < \beta. $$

From the second identity,

$$ Y < \beta \iff x\delta < \beta(D-r). $$

For simultaneous quotient and remainder recovery, the stronger condition

$$ x\delta < \beta $$

ensures that

$$ Mx \bmod \beta = Y. $$

Then

$$ \left\lfloor \frac{YD}{\beta} \right\rfloor = r. $$

The converse direction follows by reversing the same identity:

$$ \left\lfloor \frac{YD}{\beta} \right\rfloor = r + \left\lfloor \frac{x\delta}{\beta} \right\rfloor. $$

Thus exact remainder recovery forces

$$ \left\lfloor \frac{x\delta}{\beta} \right\rfloor = 0, $$

which is equivalent to

$$ x\delta < \beta. $$

Relation to direct remainder computation

The power-of-two corollary is closely related to direct remainder computation using a precomputed reciprocal, as in:

Daniel Lemire, Owen Kaser, Nathan Kurz,
Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries, arXiv:1902.01961, 2019. https://arxiv.org/abs/1902.01961

The Lean theorem proves a form that simultaneously covers both quotient recovery and remainder recovery.

Requirements

This file imports:

import Mathlib.Algebra.Order.Floor.Div
import Mathlib.Tactic

It is intended for Lean 4 with Mathlib.

ExactDivRem.lean

import Mathlib.Algebra.Order.Floor.Div
import Mathlib.Tactic

namespace ExactDivRem

/--
目的: 上側 magic number `M⁺` を定義する。
定義: `M⁺ := β ⌈/⌉ D`。
入力/前提: `β D : Nat`。
出力: 近似商 `⌊M⁺ * x / β⌋` に使う倍率。
役割: 主定理の右辺に現れる近似商を構成する。
-/
def magicUp (β D : Nat) : Nat :=
  β ⌈/⌉ D

/--
目的: 上側近似の余剰 `δ` を定義する。
定義: `M⁺ := β ⌈/⌉ D` として `δ := M⁺ * D - β`。
入力/前提: `β D : Nat`。
出力: `M⁺ * D` が `β` をどれだけ上回るかを表す自然数。
役割: 主定理の左辺 `xδ < β(D-r)` に現れる誤差項である。
-/
def delta (β D : Nat) : Nat :=
  magicUp β D * D - β

/--
目的: 上側 magic number による近似商を定義する。
定義: `M⁺ := β ⌈/⌉ D` として `⌊M⁺ * x / β⌋`。
入力/前提: `β D x : Nat`。
出力: 近似商。
役割: 主定理の右辺に現れる商である。
-/
def approxQuot (β D x : Nat) : Nat :=
  (magicUp β D * x) / β

/--
目的: 真の商と近似商が一致することを表す。
定義: `x / D = approxQuot β D x`。
入力/前提: `β D x : Nat`。
出力: 命題。
役割: residual 条件と商一致条件を結ぶ中間述語として使う。
-/
def Good (β D x : Nat) : Prop :=
  x / D = approxQuot β D x

/--
目的: 上側近似整除法で現れる補助量 `y = qδ + M⁺r` を定義する。
定義: `q = x / D`, `r = x % D`, `M⁺ = β ⌈/⌉ D`,
  `δ = M⁺ * D - β` として `qδ + M⁺r`。
入力/前提: `β D x : Nat`。
出力: 商判定・剰余判定の両方で共通に現れる自然数。
役割: `Mx = qβ + y`, `Dy = xδ + βr` の共通中間項となる。
-/
def upperY (β D x : Nat) : Nat :=
  (x / D) * delta β D + magicUp β D * (x % D)

/--
入力/前提: `A D : Nat`, `hD : 0 < D`。
主張: `A ≤ (A ⌈/⌉ D) * D`。
内容: 切り上げ除算の基本的な上側評価を、後続で使う乗法順に整える。
証明: Mathlib の `le_smul_ceilDiv` を、乗法の可換性で書き換える。
役割: `magicUp β D * D = β + delta β D` と `delta_lt` の基礎補題である。
-/
private lemma le_mul_ceilDiv {A D : Nat} (hD : 0 < D) :
    A ≤ (A ⌈/⌉ D) * D := by
  simpa [Nat.mul_comm] using
    (le_smul_ceilDiv (a := D) (b := A) hD)

/--
入力/前提: `β D : Nat`, `hD : 0 < D`。
主張: `magicUp β D * D = β + delta β D`。
内容: `delta β D` の定義を、後続の代数変形で使う等式に直す。
証明: `β ≤ magicUp β D * D` から `Nat.add_sub_of_le` で差を戻す。
役割: `M * D = β + δ` という基本変形として各分解補題で使う。
-/
private lemma magicUp_mul_eq_beta_add_delta {β D : Nat} (hD : 0 < D) :
    magicUp β D * D = β + delta β D := by
  have hle : β ≤ magicUp β D * D := by
    simpa [magicUp] using le_mul_ceilDiv (A := β) (D := D) hD
  simpa [delta, Nat.add_comm] using (Nat.add_sub_of_le hle).symm

/--
入力/前提: `x D q r : Nat`, `hD : 0 < D`, `hqr : x = q * D + r`,
  `hr : r < D`。
主張: `x / D = q`。
内容: 標準的な商剰余分解から、自然数除算の商が `q` であることを示す。
証明: `Nat.add_mul_div_right` と `Nat.div_eq_of_lt hr` に帰着する。
役割: `Good` を residual 条件に変換する補題で使う。
-/
private lemma div_of_decompose {x D q r : Nat}
    (hD : 0 < D) (hqr : x = q * D + r) (hr : r < D) :
    x / D = q := by
  rw [hqr]
  calc
    (q * D + r) / D = (r + q * D) / D := by rw [Nat.add_comm]
    _ = r / D + q := Nat.add_mul_div_right r q hD
    _ = q := by rw [Nat.div_eq_of_lt hr, zero_add]

/--
入力/前提: `β D q r : Nat`, `hD : 0 < D`。
主張: `magicUp β D * (q * D + r)` を
  `(q * delta β D + magicUp β D * r) + q * β` に分解する。
内容: `M * D = β + δ` を用いて、近似商の分子を `qβ + y` 型に変形する。
証明: 分配法則、`magicUp_mul_eq_beta_add_delta`、`ring` で正規化する。
役割: 商一致条件を residual `qδ + Mr` の大小に落とす代数補題である。
-/
private lemma approxQuot_numerator_decompose {β D q r : Nat}
    (hD : 0 < D) :
    magicUp β D * (q * D + r) =
      (q * delta β D + magicUp β D * r) + q * β := by
  calc
    magicUp β D * (q * D + r) = q * (magicUp β D * D) + magicUp β D * r := by
      rw [Nat.mul_add]
      ring
    _ = q * (β + delta β D) + magicUp β D * r := by
      rw [magicUp_mul_eq_beta_add_delta (β := β) (D := D) hD]
    _ = (q * delta β D + magicUp β D * r) + q * β := by
      ring

/--
入力/前提: `β D x q r : Nat`, `hβ : 0 < β`, `hD : 0 < D`,
  `hqr : x = q * D + r`。
主張: `approxQuot β D x = q + ((q * delta β D + magicUp β D * r) / β)`。
内容: 近似商を、真の商候補 `q` と residual の商に分ける。
証明: 分子分解後、`Nat.add_mul_div_right` で `q * β` を除去する。
役割: residual が `β` 未満であることと商一致を結ぶ。
-/
private lemma approxQuot_of_decompose {β D x q r : Nat}
    (hβ : 0 < β) (hD : 0 < D) (hqr : x = q * D + r) :
    approxQuot β D x = q + ((q * delta β D + magicUp β D * r) / β) := by
  dsimp [approxQuot]
  rw [hqr]
  rw [approxQuot_numerator_decompose (β := β) (D := D) (q := q) (r := r) hD]
  rw [Nat.add_mul_div_right _ _ hβ]
  simp [Nat.add_comm]

/--
入力/前提: `β D x q r : Nat`, `hβ : 0 < β`, `hD : 0 < D`,
  `hqr : x = q * D + r`, `hr : r < D`。
主張: `Good β D x ↔ q * delta β D + magicUp β D * r < β`。
内容: 商一致条件を、分子分解で出る residual の一つの不等式に置き換える。
証明: 真の商と近似商を展開し、`Nat.div_eq_zero_iff` で割り算を消す。
役割: 商一致の基本判定を与える中間補題である。
-/
private lemma good_iff_residual_lt {β D x q r : Nat}
    (hβ : 0 < β) (hD : 0 < D) (hqr : x = q * D + r) (hr : r < D) :
    Good β D x ↔ q * delta β D + magicUp β D * r < β := by
  rw [Good, div_of_decompose hD hqr hr, approxQuot_of_decompose hβ hD hqr]
  constructor
  · intro h
    have hq0 : q + ((q * delta β D + magicUp β D * r) / β) = q + 0 := by
      simpa using h.symm
    have hdiv : (q * delta β D + magicUp β D * r) / β = 0 := Nat.add_left_cancel hq0
    rcases (Nat.div_eq_zero_iff).1 hdiv with hβ0 | hlt
    · exact False.elim (Nat.ne_of_gt hβ hβ0)
    · exact hlt
  · intro hlt
    rw [Nat.div_eq_of_lt hlt, Nat.add_zero]

/--
入力/前提: `β D : Nat`, `hD : 0 < D`。
主張: `delta β D < D`。
内容: 切り上げ倍率の余剰 `δ = M * D - β` が 1 除数未満であることを示す。
証明: 下側評価 `β ≤ M * D` と上側評価 `M * D < β + D` を差に移す。
役割: `xδ < β` から商条件 `xδ < β(D-r)` を導く際の補助境界である。
-/
private lemma delta_lt {β D : Nat} (hD : 0 < D) :
    delta β D < D := by
  have hle : β ≤ magicUp β D * D := by
    simpa [magicUp] using le_mul_ceilDiv (A := β) (D := D) hD
  have hlt : magicUp β D * D < β + D := by
    have hmul : magicUp β D * D ≤ β + D - 1 := by
      simpa [magicUp, Nat.mul_comm] using
        Nat.div_mul_le_self (β + D - 1) D
    have hpred : β + D - 1 < β + D := by omega
    exact lt_of_le_of_lt hmul hpred
  have hlt' : magicUp β D * D < D + β := by
    simpa [Nat.add_comm] using hlt
  simpa [delta] using (Nat.sub_lt_iff_lt_add hle).2 hlt'

/--
入力/前提: `β D x : Nat`, `hD : 0 < D`。
主張: `magicUp β D * x = (x / D) * β + upperY β D x`。
内容: `x = qD + r` と `M * D = β + δ` から、`Mx = qβ + y` を得る。
証明: `Nat.div_add_mod` による商剰余分解と分子分解補題を使う。
役割: 近似商と剰余の解析で、`upperY` を `Mx` の β 未満部分として扱う。
-/
private lemma upperY_eq {β D x : Nat} (hD : 0 < D) :
    magicUp β D * x = (x / D) * β + upperY β D x := by
  let q := x / D
  let r := x % D
  have hqr : x = q * D + r := by
    dsimp [q, r]
    simpa [Nat.mul_comm, Nat.add_comm, Nat.add_left_comm, Nat.add_assoc] using
      (Nat.div_add_mod x D).symm
  calc
    magicUp β D * x = magicUp β D * (q * D + r) := by rw [hqr]
    _ = (q * delta β D + magicUp β D * r) + q * β := by
      simpa [q, r] using
        approxQuot_numerator_decompose (β := β) (D := D) (q := q) (r := r) hD
    _ = q * β + upperY β D x := by
      dsimp [upperY, q, r]
      ring

/--
入力/前提: `β D x : Nat`, `hD : 0 < D`。
主張: `D * upperY β D x = x * delta β D + β * (x % D)`。
内容: `upperY = qδ + Mr` を `D` 倍し、`M * D = β + δ` を代入する。
証明: 商剰余分解、`magicUp_mul_eq_beta_add_delta`、`ring` による正規化で示す。
役割: `upperY < β` と `xδ < β(D-r)`、剰余復元条件を接続する。
-/
private lemma upperY_scaled_eq {β D x : Nat} (hD : 0 < D) :
    D * upperY β D x = x * delta β D + β * (x % D) := by
  let q := x / D
  let r := x % D
  have hqr : x = q * D + r := by
    dsimp [q, r]
    simpa [Nat.mul_comm, Nat.add_comm, Nat.add_left_comm, Nat.add_assoc] using
      (Nat.div_add_mod x D).symm
  have hMD : magicUp β D * D = β + delta β D := by
    exact magicUp_mul_eq_beta_add_delta (β := β) (D := D) hD
  have hr : r < D := by
    dsimp [r]
    exact Nat.mod_lt _ hD
  calc
    D * upperY β D x = D * (q * delta β D + magicUp β D * r) := by
      dsimp [upperY, q, r]
    _ = q * D * delta β D + (magicUp β D * D) * r := by ring
    _ = q * D * delta β D + (β + delta β D) * r := by rw [hMD]
    _ = (q * D + r) * delta β D + β * r := by ring
    _ = x * delta β D + β * (x % D) := by
      rw [hqr]
      have hmod : (q * D + r) % D = r := by
        rw [Nat.add_comm, Nat.mul_comm q D, Nat.add_mul_mod_self_left]
        exact Nat.mod_eq_of_lt hr
      rw [hmod]

/--
入力/前提: `β D x : Nat`, `hD : 0 < D`, `hβ : 0 < β`。
主張: `approxQuot β D x = x / D ↔ upperY β D x < β`。
内容: `Mx = qβ + upperY` から、近似商の余分な項が 0 である条件に直す。
証明: `upperY_eq` で分子を展開し、`Nat.div_eq_zero_iff` を使う。
役割: 主定理の商一致側を共通補助量 `upperY` の不等式へ変換する。
-/
private lemma approxQuot_eq_div_iff_upperY_lt {β D x : Nat}
    (hD : 0 < D) (hβ : 0 < β) :
    approxQuot β D x = x / D ↔ upperY β D x < β := by
  have hdiv : approxQuot β D x = x / D + upperY β D x / β := by
    calc
      approxQuot β D x = (magicUp β D * x) / β := rfl
      _ = ((x / D) * β + upperY β D x) / β := by rw [upperY_eq (β := β) (D := D) (x := x) hD]
      _ = (upperY β D x + (x / D) * β) / β := by rw [Nat.add_comm]
      _ = upperY β D x / β + x / D := Nat.add_mul_div_right (upperY β D x) (x / D) hβ
      _ = x / D + upperY β D x / β := by rw [Nat.add_comm]
  rw [hdiv]
  constructor
  · intro hEq
    have h0 : x / D + upperY β D x / β = x / D + 0 := by simpa using hEq
    have hydiv0 : upperY β D x / β = 0 := Nat.add_left_cancel h0
    rcases (Nat.div_eq_zero_iff).1 hydiv0 with hβ0 | hylt
    · exact False.elim (Nat.ne_of_gt hβ hβ0)
    · exact hylt
  · intro hylt
    rw [Nat.div_eq_of_lt hylt, Nat.add_zero]

/--
入力/前提: `β D x : Nat`, `hD : 0 < D`。
主張: `upperY β D x < β ↔ x * delta β D < β * (D - x % D)`。
内容: `D * upperY = xδ + βr` を使い、商一致条件を元の不等式へ戻す。
証明: 正の `D` を両辺に掛け、`upperY_scaled_eq` と自然数線形算術で整理する。
役割: `approxQuot_eq_div_iff_mul_delta_lt` の左辺条件を導く最終変換である。
-/
private lemma upperY_lt_iff_mul_delta_lt {β D x : Nat} (hD : 0 < D) :
    upperY β D x < β ↔ x * delta β D < β * (D - x % D) := by
  calc
    upperY β D x < β ↔ D * upperY β D x < D * β := by
      exact (Nat.mul_lt_mul_left hD).symm
    _ ↔ x * delta β D + β * (x % D) < D * β := by
      rw [upperY_scaled_eq (β := β) (D := D) (x := x) hD]
    _ ↔ x * delta β D + β * (x % D) < β * D := by rw [Nat.mul_comm D β]
    _ ↔ x * delta β D < β * D - β * (x % D) := by
      omega
    _ ↔ x * delta β D < β * (D - x % D) := by
      rw [Nat.mul_sub_left_distrib]

/--
目的: `M := β ⌈/⌉ D`, `δ := M * D - β`, `q := x / D`, `r := x % D`
  に対し、`xδ < β(D-r)` と `⌊M * x / β⌋ = q` の同値を与える。
定義: 定理の主張内で `let M := β ⌈/⌉ D` と
  `let δ := M * D - β` を展開する。
入力/前提: `β D x : Nat`, `hD : 0 < D`, `hβ : 0 < β`。
出力: 両条件の論理同値。
役割: 上側 Barrett 型 1 段除算の厳密一致判定を、
  `xδ < β(D-r)` の形で直接使えるようにする。
-/
theorem approxQuot_eq_div_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    x * δ < β * (D - r) ↔ (M * x) / β = q := by
  dsimp
  calc
    x * delta β D < β * (D - x % D) ↔ upperY β D x < β :=
      (upperY_lt_iff_mul_delta_lt (β := β) (D := D) (x := x) hD).symm
    _ ↔ approxQuot β D x = x / D :=
      (approxQuot_eq_div_iff_upperY_lt (β := β) (D := D) (x := x) hD hβ).symm

/--
目的: 上側 magic number による商と剰余の同時復元条件を示す。
定義: `M := β ⌈/⌉ D`, `δ := M * D - β`, `q := x / D`, `r := x % D`,
  `y := M * x % β` とおく。
入力/前提: `β D x : Nat`, `hD : 0 < D`, `hβ : 0 < β`。
主張: `xδ < β` は
  `⌊M * x / β⌋ = q ∧ ⌊y * D / β⌋ = r`
  と同値である。
内容: 共通補助量 `upperY = qδ + M⁺r` を使い、
  `M * x = qβ + upperY` と `D * upperY = xδ + βr` を利用する。
  前向きでは `xδ < β` から商一致を得て `M * x % β = upperY` とし、
  `D * upperY` を `β` で割って剰余一致を導く。
  逆向きでは商一致から `upperY < β` と `M * x % β = upperY` を得て、
  剰余一致と `D * upperY = xδ + βr` から `(xδ) / β = 0` を回収する。
証明: `approxQuot_eq_div_iff_upperY_lt`,
  `upperY_lt_iff_mul_delta_lt`, `upperY_eq`, `upperY_scaled_eq` に帰着する。
役割: 上側 magic number による 1 段 `divmod` の正しさを、
  商と剰余を同時に含む形で参照できる最終定理である。
-/
theorem approxDivRem_iff_mul_delta_lt
    {β D x : Nat} (hD : 0 < D) (hβ : 0 < β) :
    let M := β ⌈/⌉ D
    let δ := M * D - β
    let q := x / D
    let r := x % D
    let y := (M * x) % β
    x * δ < β ↔ ((M * x) / β = q ∧ ((y * D) / β = r)) := by
  dsimp
  set M := magicUp β D
  set δ := delta β D
  set q := x / D
  set r := x % D
  set y := (M * x) % β
  set y0 := upperY β D x
  have hr : r < D := by
    dsimp [r]
    exact Nat.mod_lt _ hD
  have hMx0 : M * x = q * β + y0 := by
    simpa [M, q, y0] using upperY_eq (β := β) (D := D) (x := x) hD
  have hDy0 : D * y0 = x * δ + β * r := by
    simpa [δ, r, y0] using upperY_scaled_eq (β := β) (D := D) (x := x) hD
  have hquot0 : (M * x) / β = q ↔ y0 < β := by
    simpa [M, q, y0] using
      approxQuot_eq_div_iff_upperY_lt (β := β) (D := D) (x := x) hD hβ
  have hquot : (M * x) / β = q ↔ x * δ < β * (D - r) := by
    calc
      (M * x) / β = q ↔ y0 < β := hquot0
      _ ↔ x * δ < β * (D - r) := by
        simpa [δ, r, y0] using
          upperY_lt_iff_mul_delta_lt (β := β) (D := D) (x := x) hD
  constructor
  · intro hxδ
    have hxδ' : x * δ < β * (D - r) := by
      have hDr1 : 1 ≤ D - r := by omega
      calc
        x * δ < β := hxδ
        _ ≤ β * (D - r) := by simpa using Nat.mul_le_mul_left β hDr1
    have hqEq : (M * x) / β = q := hquot.mpr hxδ'
    have hxδ_delta : x * δ < β := by
      simpa [δ, delta, M, magicUp] using hxδ
    have hy0_lt : y0 < β := hquot0.mp hqEq
    have hy : y = y0 := by
      calc
        y = (q * β + y0) % β := by
          dsimp [y]
          rw [hMx0]
        _ = (y0 + q * β) % β := by rw [Nat.add_comm]
        _ = (y0 + β * q) % β := by rw [Nat.mul_comm q β]
        _ = y0 % β := Nat.add_mul_mod_self_left y0 β q
        _ = y0 := Nat.mod_eq_of_lt hy0_lt
    have hrEq : (y * D) / β = r := by
      calc
        (y * D) / β = (y0 * D) / β := by rw [hy]
        _ = (x * δ + r * β) / β := by
          rw [Nat.mul_comm y0 D, hDy0, Nat.mul_comm β r]
        _ = (x * δ) / β + r := Nat.add_mul_div_right (x * δ) r hβ
        _ = r := by rw [Nat.div_eq_of_lt hxδ_delta, zero_add]
    exact ⟨hqEq, hrEq⟩
  · rintro ⟨hqEq, hrEq⟩
    have hy0_lt : y0 < β := hquot0.mp hqEq
    have hy : y = y0 := by
      calc
        y = (q * β + y0) % β := by
          dsimp [y]
          rw [hMx0]
        _ = (y0 + q * β) % β := by rw [Nat.add_comm]
        _ = (y0 + β * q) % β := by rw [Nat.mul_comm q β]
        _ = y0 % β := Nat.add_mul_mod_self_left y0 β q
        _ = y0 := Nat.mod_eq_of_lt hy0_lt
    have hdiv_formula : (y * D) / β = r + (x * δ) / β := by
      calc
        (y * D) / β = (y0 * D) / β := by rw [hy]
        _ = (x * δ + r * β) / β := by
          rw [Nat.mul_comm y0 D, hDy0, Nat.mul_comm β r]
        _ = (x * δ + r * β) / β := by rw [Nat.add_comm]
        _ = r + (x * δ) / β := by
          rw [Nat.add_mul_div_right (x * δ) r hβ, Nat.add_comm]
    have hxδ_div0 : (x * δ) / β = 0 := by
      have htmp : r + (x * δ) / β = r + 0 := by
        calc
          r + (x * δ) / β = (y * D) / β := by symm; exact hdiv_formula
          _ = r := hrEq
          _ = r + 0 := by simp
      exact Nat.add_left_cancel htmp
    rcases (Nat.div_eq_zero_iff).1 hxδ_div0 with hβ0 | hxlt
    · exact False.elim (Nat.ne_of_gt hβ hβ0)
    · exact hxlt

/--
入力/前提: `N L D x : Nat`, `hD : 0 < D`, `hDlt : D < 2 ^ N`,
  `hδ : ((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L) ≤ 2 ^ L`,
  `hx : x < 2 ^ N`。
主張: `β := 2 ^ (N + L)` と `M := β ⌈/⌉ D` に対し、
  近似商 `(M * x) / β` と剰余復元 `(((M * x) % β) * D) / β` が正しい。
内容: `δ ≤ 2 ^ L` なら、任意の `N`-bit 入力 `x` で商と剰余を同時に復元できる。
証明: `xδ ≤ (2 ^ N - 1) * 2 ^ L < 2 ^ (N + L)` を示し、
  `approxDivRem_iff_mul_delta_lt` に渡す。
役割: `β = 2 ^ (N + L)` の固定基数で使う系であり、剰余式は次の資料の
  Theorem 1 と同じ形を与える。
参照: Daniel Lemire, Owen Kaser, Nathan Kurz,
  "Faster Remainder by Direct Computation: Applications to Compilers and
  Software Libraries", arXiv:1902.01961, 2019.
  URL: https://arxiv.org/abs/1902.01961
-/
theorem approxDivRem_pow2_of_delta_le
    {N L D x : Nat} (hD : 0 < D) (hDlt : D < 2 ^ N)
    (hδ : ((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L) ≤ 2 ^ L)
    (hx : x < 2 ^ N) :
    let β := 2 ^ (N + L)
    let M := β ⌈/⌉ D
    let y := (M * x) % β
    (M * x) / β = x / D ∧ (y * D) / β = x % D := by
  have hNpos : 0 < 2 ^ N := lt_trans hD hDlt
  dsimp
  have hβ : 0 < 2 ^ (N + L) := Nat.two_pow_pos _
  have hxle : x ≤ 2 ^ N - 1 := by
    simpa [Nat.pred_eq_sub_one] using Nat.le_pred_of_lt hx
  have hpred : 2 ^ N - 1 < 2 ^ N := by
    exact Nat.sub_lt hNpos (by decide)
  have hxδ : x * (((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L)) < 2 ^ (N + L) := by
    calc
      x * (((2 ^ (N + L)) ⌈/⌉ D) * D - 2 ^ (N + L)) ≤ (2 ^ N - 1) * 2 ^ L := by
        exact Nat.mul_le_mul hxle hδ
      _ < 2 ^ N * 2 ^ L := by
        exact Nat.mul_lt_mul_of_pos_right hpred (Nat.two_pow_pos L)
      _ = 2 ^ (N + L) := by
        rw [Nat.pow_add]
  have hmain :=
    (approxDivRem_iff_mul_delta_lt
      (β := 2 ^ (N + L)) (D := D) (x := x) hD hβ).mp hxδ
  exact hmain

end ExactDivRem