# 0.999... = 1 -- a machine-checked proof, in one file.
#
# Bend is a Python-shaped language where a TYPE is a claim and a PROGRAM of
# that type is its proof. There are no tactics: every proof below is an
# ordinary recursive function, and `bend reals.bend --check` verifies them.
#
# The plan, in four steps:
#
#   1. Integers.   Pos{n} is n, NegS{n} is -(1+n). One term per integer.
#   2. Rationals.  Q{num, den} means num / (1 + den), so the denominator is
#                  never zero. Different terms may denote the same number.
#   3. Reals.      A real is a RECIPE -- a little piece of data saying how
#                  the number was built -- and `Real.at(x, n)` is the
#                  interpreter that reads the recipe and hands back the
#                  n-th rational approximation, good to within 2^-n.
#   4. Equality.   Two reals are equal when, for EVERY precision n, their
#                  n-th approximations sit within 2^-n + 2^-n of each other.
#                  Nothing ever has to compare infinitely many digits.
#
# Then 0.999... is the recipe RNines{}, whose n-th approximation is
# 1 - 10^-n, and the theorem `nines_is_one` says it is equal to 1. The
# whole numerical content is one inequality: the gap 10^-n never exceeds
# the allowance 2^-n + 2^-n, because 2^n <= 10^n.

import Base
import ./nat.bend as N

# ===================================================================
# 0. Two facts about naturals that lib/nat.bend does not ship
# ===================================================================

# N.Nat.Le(a, b) is a COMPUTED proposition: it reduces to Unit when a <= b
# and to Empty when it does not. So a proof of it is data -- but in Bend
# every variable is spent when used, and this one is needed twice below.
# le_copy re-derives a fresh pair by walking the two numbers down together.
def le_copy(a: Nat, b: Nat, h: N.Nat.Le(a, b), -P: Type, k: N.Nat.Le(a, b) -> N.Nat.Le(a, b) -> P) -> P:
  match a:
    case 0n:
      k(Unit{}, Unit{})
    case 1n+p:
      match b:
        case 0n:
          match h:
        case 1n+q:
          le_copy(p, q, h, P, k)

# x <= y  =>  c*x <= c*y
def le_mul_l(+c: Nat, +x: Nat, +y: Nat, h: N.Nat.Le(x, y)) -> N.Nat.Le(Nat.mul(c, x), Nat.mul(c, y)):
  match c:
    case 0n:
      Unit{}
    case 1n+p:
      le_copy(x, y, h, N.Nat.Le(Nat.mul(1n+p, x), Nat.mul(1n+p, y)), h1 => h2 =>
        N.Nat.le_add_both(x, y, Nat.mul(p, x), Nat.mul(p, y), h1, le_mul_l(p, x, y, h2)))

# x <= y  =>  x*c <= y*c
def le_mul_r(x: Nat, y: Nat, +c: Nat, h: N.Nat.Le(x, y)) -> N.Nat.Le(Nat.mul(x, c), Nat.mul(y, c)):
  match x y:
    case 0n 0n:
      Unit{}
    case 0n 1n+yp:
      Unit{}
    case 1n+xp 0n:
      match h:
    case 1n+xp 1n+yp:
      N.Nat.le_add_mono_l(c, Nat.mul(xp, c), Nat.mul(yp, c), le_mul_r(xp, yp, c, h))

# ===================================================================
# 1. The integers
# ===================================================================
#
# Pos{n} is n and NegS{n} is -(1+n), so every integer has exactly ONE
# representation and `==` on Int is honest equality.

type Int is Data:
  Pos{n: Nat}
  NegS{n: Nat}

# -n, for a natural n
def Int.neg_nat(n: Nat) -> Int:
  match n:
    case 0n:
      Pos{0n}
    case 1n+p:
      NegS{p}

# a - b for two naturals, as an integer: count both down until one runs out
def Int.sub_nat(a: Nat, b: Nat) -> Int:
  match a b:
    case 0n 0n:
      Pos{0n}
    case 0n 1n+bp:
      NegS{bp}
    case 1n+ap 0n:
      Pos{1n+ap}
    case 1n+ap 1n+bp:
      Int.sub_nat(ap, bp)

def Int.neg(a: Int) -> Int:
  match a:
    case Pos{n}:
      Int.neg_nat(n)
    case NegS{n}:
      Pos{1n+n}

def Int.add(a: Int, b: Int) -> Int:
  match a b:
    case Pos{x} Pos{y}:
      Pos{Nat.add(x, y)}
    case Pos{x} NegS{y}:
      Int.sub_nat(x, 1n+y)
    case NegS{x} Pos{y}:
      Int.sub_nat(y, 1n+x)
    case NegS{x} NegS{y}:
      NegS{1n+Nat.add(x, y)}

def Int.mul(a: Int, b: Int) -> Int:
  match a b:
    case Pos{x} Pos{y}:
      Pos{Nat.mul(x, y)}
    case Pos{x} NegS{y}:
      Int.neg_nat(Nat.mul(x, 1n+y))
    case NegS{x} Pos{y}:
      Int.neg_nat(Nat.mul(1n+x, y))
    case NegS{x} NegS{y}:
      Pos{Nat.mul(1n+x, 1n+y)}

def Int.abs(a: Int) -> Nat:
  match a:
    case Pos{n}:
      n
    case NegS{n}:
      1n+n

# a <= b, again as a computed proposition
def Int.Le(a: Int, b: Int) -> Type:
  match a b:
    case Pos{x} Pos{y}:
      N.Nat.Le(x, y)
    case Pos{x} NegS{y}:
      Empty
    case NegS{x} Pos{y}:
      Unit
    case NegS{x} NegS{y}:
      N.Nat.Le(y, x)

def Int.sub_nat_zero(x: Nat) -> {Int.sub_nat(x, 0n) == Pos{x} : Int}:
  match x:
    case 0n:
      {==}
    case 1n+xp:
      {==}

# ===================================================================
# 2. The rationals
# ===================================================================
#
# Q{num, den} means num / (1 + den). Storing "denominator minus one"
# is what makes division by zero unrepresentable rather than merely
# forbidden.

type Q is Data:
  Q{num: Int, den: Nat}

def Q.of_nat(n: Nat) -> Q:
  Q{Pos{n}, 0n}

def Q.num(q: Q) -> Int:
  match q:
    case Q{n, d}:
      n

def Q.den_int(q: Q) -> Int:
  match q:
    case Q{n, d}:
      Pos{1n+d}

# 2/4 and 1/2 are different TERMS but the same number, so equality of
# rationals is this relation -- cross multiplication -- not `==`.
# (The proof below never needs it; it is here because "what equality
# means" is the whole subject of this file.)
def Q.Eq(+x: Q, +y: Q) -> Type:
  {Int.mul(Q.num(x), Q.den_int(y)) == Int.mul(Q.num(y), Q.den_int(x)) : Int}

# the denominator of a product, in "minus one" form. The point of writing
# it this way: 1 + Q.dmul(b, d) and Nat.mul(1n+b, 1n+d) unfold to the SAME
# term, so no lemma is ever needed to move between the two.
def Q.dmul(b: Nat, +d: Nat) -> Nat:
  Nat.add(d, Nat.mul(b, 1n+d))

# a/(1+e1) + b/(1+e2) = (a(1+e2) + b(1+e1)) / ((1+e1)(1+e2))
def Q.add(a: Q, b: Q) -> Q:
  match a b:
    case Q{n1, d1} Q{n2, d2}:
      +e1 = d1
      +e2 = d2
      Q{Int.add(Int.mul(n1, Pos{1n+e2}), Int.mul(n2, Pos{1n+e1})), Q.dmul(e1, e2)}

def Q.neg(a: Q) -> Q:
  match a:
    case Q{n, d}:
      Q{Int.neg(n), d}

def Q.sub(a: Q, b: Q) -> Q:
  Q.add(a, Q.neg(b))

def Q.abs(a: Q) -> Q:
  match a:
    case Q{n, d}:
      Q{Pos{Int.abs(n)}, d}

# a/(1+d1) <= b/(1+d2) is decided by cross multiplication
def Q.Le(a: Q, b: Q) -> Type:
  match a b:
    case Q{n1, d1} Q{n2, d2}:
      Int.Le(Int.mul(n1, Pos{1n+d2}), Int.mul(n2, Pos{1n+d1}))

# ===================================================================
# 3. Powers, and the two rationals the theorem is about
# ===================================================================
#
# The checker evaluates Nat in unary, so a literal like 1024 is out of
# reach: powers are always built by recursion, never written down. Both
# functions return the power MINUS ONE, because that is exactly what the
# `den` field of Q wants.

# two_pow_m1(n) = 2^n - 1, so 1n+two_pow_m1(n) is 2^n
def two_pow_m1(n: Nat) -> Nat:
  match n:
    case 0n:
      0n
    case 1n+p:
      1n+Nat.double(two_pow_m1(p))

# ten_pow_m1(n) = 10^n - 1, so 1n+ten_pow_m1(n) is 10^n
def ten_pow_m1(n: Nat) -> Nat:
  match n:
    case 0n:
      0n
    case 1n+p:
      9n+Nat.mul(10n, ten_pow_m1(p))

# eps(n) = 2^-n : the error allowed at precision n
def eps(n: Nat) -> Q:
  Q{Pos{1n}, two_pow_m1(n)}

# ten_pow_inv(n) = 10^-n
def ten_pow_inv(n: Nat) -> Q:
  Q{Pos{1n}, ten_pow_m1(n)}

# ===================================================================
# 4. The reals: recipes, and their interpreter
# ===================================================================
#
# A real number is not a function and not a stream -- it is DATA: a small
# recipe saying how the number was built. Here two constructors suffice:
#
#   RQ{q}     -- the rational q, exactly
#   RNines{}  -- 0.999...
#
# (the full library also has sums, products, sqrt 2 and pi; they are
# dropped here because this theorem does not use them.)
#
# Storing the recipe rather than an approximation function is what makes
# everything below possible: data can be COPIED, so a proof may ask the
# same real for as many different precisions as it likes.

type Real is Data:
  RQ{q: Q}
  RNines{}

def Real.of_nat(n: Nat) -> Real:
  RQ{Q.of_nat(n)}

# THE INTERPRETER. Real.at(x, n) reads the recipe x and returns its n-th
# rational approximation. A rational is its own approximation at every
# precision; 0.999... rounded to n decimal places is 0.99...9, which is
# exactly 1 - 10^-n.
def Real.at(x: Real, n: Nat) -> Q:
  match x:
    case RQ{q}:
      q
    case RNines{}:
      Q.sub(Q.of_nat(1n), ten_pow_inv(n))

# WHAT EQUALITY MEANS. Two reals are equal when, at every precision n,
# their approximations differ by at most 2^-n + 2^-n -- each is allowed
# its own 2^-n of slack. This is a claim about all n at once, so it is
# proved, never computed; that is why 0.999... = 1 is a theorem.
def Real.Eq(x: Real, y: Real) -> Type:
  @n: Nat -> Q.Le(Q.abs(Q.sub(Real.at(x, n), Real.at(y, n))), Q.add(eps(n), eps(n)))

# ===================================================================
# 5. The proof
# ===================================================================
#
# Step 1 (nines_spec_eq): 1 - 10^-n, computed by Q.sub, comes out as the
#   term Q{Pos{10^n - 1}, 10^n - 1} -- i.e. (10^n - 1)/10^n. This is pure
#   bookkeeping: the arithmetic has already happened, we are only naming
#   the result.
#
# Step 2 (nines_one_num): the numerator of the difference
#   (1 - 10^-n) - 1 collapses to -1, so the difference is exactly -10^-n
#   and its absolute value is 10^-n.
#
# Step 3 (pow_le): 2^n <= 10^n, by induction. Doubling is at most
#   multiplying by ten.
#
# Step 4 (nines_one_est): therefore 10^-n <= 2^-n <= 2^-n + 2^-n.
#   Cross-multiplied and cleared of ones, that is the single inequality
#   2^n * 2^n <= (2^n + 2^n) * 10^n, which follows from step 3.
#
# The gap never exceeds the allowance, at any precision. So the two
# numbers are equal.

# (10^n - 1) / 10^n, which is what "0.99...9 with n nines" is
def nines_at(+k: Nat) -> Q:
  Q{Pos{ten_pow_m1(k)}, ten_pow_m1(k)}

def mkq(+z: Nat) -> Q:
  Q{Pos{z}, z}

# STEP 1: the interpreter's answer really is (10^n - 1)/10^n
def nines_spec_eq(+k: Nat) -> {Q.sub(Q.of_nat(1n), ten_pow_inv(k)) == nines_at(k) : Q}:
  %Equal.sym(Int, Int.sub_nat(Nat.add(ten_pow_m1(k), 0n), 0n), Pos{Nat.add(ten_pow_m1(k), 0n)}, Int.sub_nat_zero(Nat.add(ten_pow_m1(k), 0n))) : {Q{_, Nat.add(ten_pow_m1(k), 0n)} == nines_at(k) : Q}
  Equal.cong(Nat, Q, z => mkq(z), Nat.add(ten_pow_m1(k), 0n), ten_pow_m1(k), N.Nat.add_zero(ten_pow_m1(k)))

# 10^(p+1) = 10 * 10^p
def ten_succ(+p: Nat) -> {1n+ten_pow_m1(1n+p) == Nat.mul(10n, 1n+ten_pow_m1(p)) : Nat}:
  Equal.sym(Nat, Nat.mul(10n, 1n+ten_pow_m1(p)), Nat.add(10n, Nat.mul(10n, ten_pow_m1(p))), N.Nat.mul_succ(10n, ten_pow_m1(p)))

# STEP 3: 2^k <= 10^k. Induction: 2^(p+1) = 2 * 2^p <= 2 * 10^p <= 10 * 10^p.
def pow_le(+k: Nat) -> N.Nat.Le(1n+two_pow_m1(k), 1n+ten_pow_m1(k)):
  match k:
    case 0n:
      Unit{}
    case 1n+p:
      %N.Nat.two_mul(1n+two_pow_m1(p)) : N.Nat.Le(_, 1n+ten_pow_m1(1n+p))
      %Equal.sym(Nat, 1n+ten_pow_m1(1n+p), Nat.mul(10n, 1n+ten_pow_m1(p)), ten_succ(p)) : N.Nat.Le(Nat.mul(2n, 1n+two_pow_m1(p)), _)
      N.Nat.le_trans(Nat.mul(2n, 1n+two_pow_m1(p)), Nat.mul(2n, 1n+ten_pow_m1(p)), Nat.mul(10n, 1n+ten_pow_m1(p)),
        le_mul_l(2n, 1n+two_pow_m1(p), 1n+ten_pow_m1(p), pow_le(p)),
        le_mul_r(2n, 10n, 1n+ten_pow_m1(p), Unit{}))

# x + (-y) = x - y, for naturals cast into Int
def pos_add_neg(x: Nat, y: Nat) -> {Int.add(Pos{x}, Int.neg(Pos{y})) == Int.sub_nat(x, y) : Int}:
  match x y:
    case 0n 0n:
      {==}
    case 0n 1n+q:
      {==}
    case 1n+p 0n:
      Equal.cong(Nat, Int, z => Pos{1n+z}, Nat.add(p, 0n), p, N.Nat.add_zero(p))
    case 1n+p 1n+q:
      {==}

# A - (A+1) = -1
def sub_succ(A: Nat) -> {Int.sub_nat(A, 1n+A) == NegS{0n} : Int}:
  match A:
    case 0n:
      {==}
    case 1n+p:
      sub_succ(p)

def sub_one(+A: Nat) -> {Int.sub_nat(Nat.mul(A, 1n), Nat.mul(1n, 1n+A)) == NegS{0n} : Int}:
  %Equal.sym(Nat, Nat.mul(A, 1n), A, N.Nat.mul_one(A)) : {Int.sub_nat(_, Nat.mul(1n, 1n+A)) == NegS{0n} : Int}
  %Equal.sym(Nat, Nat.mul(1n, 1n+A), 1n+A, N.Nat.one_mul(1n+A)) : {Int.sub_nat(A, _) == NegS{0n} : Int}
  sub_succ(A)

# (-1) * b = -(1 * b): both sides compute to the same term
def neg_one_mul(+A: Nat) -> {Int.mul(Int.neg(Pos{1n}), Pos{1n+A}) == Int.neg(Int.mul(Pos{1n}, Pos{1n+A})) : Int}:
  {==}

# STEP 2: the numerator of (1 - 10^-n) - 1 is exactly -1.
# With A = 10^n - 1 this is A*1 + (-1)*(1+A) = -1.
def nines_one_num(+A: Nat) -> {Int.add(Int.mul(Pos{A}, Pos{1n+0n}), Int.mul(Int.neg(Pos{1n}), Pos{1n+A})) == NegS{0n} : Int}:
  %Equal.sym(Int, Int.mul(Int.neg(Pos{1n}), Pos{1n+A}), Int.neg(Int.mul(Pos{1n}, Pos{1n+A})), neg_one_mul(A)) : {Int.add(Pos{Nat.mul(A, 1n)}, _) == NegS{0n} : Int}
  Equal.trans(Int, Int.add(Pos{Nat.mul(A, 1n)}, Int.neg(Pos{Nat.mul(1n, 1n+A)})), Int.sub_nat(Nat.mul(A, 1n), Nat.mul(1n, 1n+A)), NegS{0n},
    pos_add_neg(Nat.mul(A, 1n), Nat.mul(1n, 1n+A)),
    sub_one(A))

# STEP 4: the estimate, cross-multiplied. With A = 10^n and Sn+1 = 2^n:
#   2^n * 2^n <= (2^n + 2^n) * 10^n
# because 2^n <= 10^n (h) and x <= x + x. The `%` lines only strip the
# multiplications by one that Q.Le left behind.
def nines_one_est(+A: Nat, +Sn: Nat, h: N.Nat.Le(1n+Sn, 1n+A)) -> N.Nat.Le(Nat.mul(1n, Nat.mul(1n+Sn, 1n+Sn)), Nat.mul(Nat.add(Nat.mul(1n, 1n+Sn), Nat.mul(1n, 1n+Sn)), Nat.mul(1n+A, 1n))):
  %Equal.sym(Nat, Nat.mul(1n, Nat.mul(1n+Sn, 1n+Sn)), Nat.mul(1n+Sn, 1n+Sn), N.Nat.one_mul(Nat.mul(1n+Sn, 1n+Sn))) : N.Nat.Le(_, Nat.mul(Nat.add(Nat.mul(1n, 1n+Sn), Nat.mul(1n, 1n+Sn)), Nat.mul(1n+A, 1n)))
  %Equal.sym(Nat, Nat.mul(1n, 1n+Sn), 1n+Sn, N.Nat.one_mul(1n+Sn)) : N.Nat.Le(Nat.mul(1n+Sn, 1n+Sn), Nat.mul(Nat.add(_, Nat.mul(1n, 1n+Sn)), Nat.mul(1n+A, 1n)))
  %Equal.sym(Nat, Nat.mul(1n, 1n+Sn), 1n+Sn, N.Nat.one_mul(1n+Sn)) : N.Nat.Le(Nat.mul(1n+Sn, 1n+Sn), Nat.mul(Nat.add(1n+Sn, _), Nat.mul(1n+A, 1n)))
  %Equal.sym(Nat, Nat.mul(1n+A, 1n), 1n+A, N.Nat.mul_one(1n+A)) : N.Nat.Le(Nat.mul(1n+Sn, 1n+Sn), Nat.mul(Nat.add(1n+Sn, 1n+Sn), _))
  %Equal.sym(Nat, Nat.mul(Nat.add(1n+Sn, 1n+Sn), 1n+A), Nat.add(Nat.mul(1n+Sn, 1n+A), Nat.mul(1n+Sn, 1n+A)), N.Nat.add_mul(1n+Sn, 1n+Sn, 1n+A)) : N.Nat.Le(Nat.mul(1n+Sn, 1n+Sn), _)
  N.Nat.le_trans(Nat.mul(1n+Sn, 1n+Sn), Nat.mul(1n+Sn, 1n+A), Nat.add(Nat.mul(1n+Sn, 1n+A), Nat.mul(1n+Sn, 1n+A)),
    le_mul_l(1n+Sn, 1n+Sn, 1n+A, h),
    N.Nat.le_add_r(Nat.mul(1n+Sn, 1n+A), Nat.mul(1n+Sn, 1n+A)))

# |0.99...9 - 1| <= 2^-n + 2^-n : steps 2, 3 and 4 put together
def nines_one_at(+n: Nat) -> Q.Le(Q.abs(Q.sub(nines_at(n), Q.of_nat(1n))), Q.add(eps(n), eps(n))):
  %Equal.sym(Int, Int.add(Int.mul(Pos{ten_pow_m1(n)}, Pos{1n+0n}), Int.mul(Int.neg(Pos{1n}), Pos{1n+ten_pow_m1(n)})), NegS{0n}, nines_one_num(ten_pow_m1(n))) : Q.Le(Q.abs(Q{_, Q.dmul(ten_pow_m1(n), 0n)}), Q.add(eps(n), eps(n)))
  nines_one_est(ten_pow_m1(n), two_pow_m1(n), pow_le(n))

# the same statement, phrased in the interpreter's own words
def nines_is_one_at(+n: Nat) -> Q.Le(Q.abs(Q.sub(Q.sub(Q.of_nat(1n), ten_pow_inv(n)), Q.of_nat(1n))), Q.add(eps(n), eps(n))):
  %Equal.sym(Q, Q.sub(Q.of_nat(1n), ten_pow_inv(n)), nines_at(n), nines_spec_eq(n)) : Q.Le(Q.abs(Q.sub(_, Q.of_nat(1n))), Q.add(eps(n), eps(n)))
  nines_one_at(n)

# ===================================================================
# 0.999... = 1
# ===================================================================
#
# Real.Eq unfolds to "for every n, ..."; the proof is the function that
# takes a precision n and hands back the estimate at that precision.

def nines_is_one() -> Real.Eq(RNines{}, Real.of_nat(1n)):
  n => nines_is_one_at(n)