~/bend-docscommunity

reals.bend source

reals.bend on the hub · documented module

# 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 Baseimport ./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*ydef 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*cdef 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 ndef 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 outdef 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 propositiondef 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}:      ndef 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 multiplicationdef 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^ndef 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^ndef 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 ndef eps(n: Nat) -> Q:  Q{Pos{1n}, two_pow_m1(n)}# ten_pow_inv(n) = 10^-ndef 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" isdef 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^ndef 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^pdef 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 Intdef 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) = -1def 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 termdef 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 togetherdef 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 wordsdef 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)