Sqrt 2 by Rotsor

module Theorem where

open import Data.Nat
open import Data.Bool
open import Data.Empty
open import Data.Sum
open import Data.Product
open import Relation.Nullary
open import Relation.Binary.PropositionalEquality using (_≡_; refl; cong; cong₂; trans; sym; subst)
open import Data.Nat.Properties as NatProp
open import Relation.Binary


mutual
  data Even : ℕ → Set where
    zero : Even zero
    even : ∀ {n} → Odd n → Even (suc n)

  data Odd : ℕ → Set where
    odd : ∀ {n} → Even n → Odd (suc n)

2^_ : ℕ → ℕ
2^ zero = 1
2^ (suc n) = 2 * 2^ n

bit2ℕ : Bool → ℕ
bit2ℕ true = 1
bit2ℕ false = 0

lemma : ∀ {m n} → m ≤ n → m ≤ suc n
lemma leq = ≤′⇒≤ (≤′-step (≤⇒≤′ leq))

theorem : ∀ x y → x + (suc y) ≡ suc (x + y)
theorem zero y = refl
theorem (suc x) y = cong suc (theorem x y)

div2 : ∀ n → n ≡ 0 ⊎ ∃₂ λ m r → 2 * m + bit2ℕ r ≡ n × m < n
div2 (suc n) with div2 n
div2 (suc .0) | inj₁ refl = inj₂ (0 , true , refl , s≤s z≤n)
div2 (suc _) | inj₂ (m , false , eq , le) = 
  inj₂ (m , true , 
       trans (theorem (m + (m + 0)) 0) (cong suc eq) , 
       lemma le)
div2 (suc _) | inj₂ (m , true , eq , le) = 
  inj₂ (suc m , false , cong suc ( trans (cong (λ z → z + 0) eq'') eq') , s≤s le) where
  eq' = trans (sym (theorem (m + (m + 0)) 0)) eq
  eq'' = theorem m (m + 0)
div2 zero = inj₁ refl

wf-rec : (P : ℕ → Set) → 
  ((n : ℕ) → P n ⊎ (∃ (λ n' → n' < n × (P n' → P n))))
  → ∀ n → P n
wf-rec P pg n = go n n ≤′-refl where
  go : ∀ m n → n ≤′ m → P n
  go .(suc m) n (≤′-step {m} le) = go m n le
  go n .n ≤′-refl with pg n 
  go n .n ≤′-refl | inj₁ p = p
  go (suc m) .(suc m) ≤′-refl | inj₂ (n' , s≤s lesss , implies) = 
    implies (go m n' (≤⇒≤′ lesss))
  go 0 .0 ≤′-refl | inj₂ (n' , () , implies)

Twos : ℕ → ℕ → Set
Twos n k = ∃ λ x → Odd x × (2^ k * x ≡ n)

Decomposition = λ n → n ≡ 0 ⊎ ∃ (Twos n)

import Algebra
module S = Algebra.CommutativeSemiring (NatProp.commutativeSemiring)

lemma00 : ∀ {m} → Even m → Even (2 + m)
lemma00 e = even (odd e)

cong' : ∀ {A B : Set} (f : A → B) → ∀ {x y} → x ≡ y → f x ≡ f y
cong' f  {x} {y} eq = subst (λ z → f x ≡ f z) eq refl

{- cong P eq : P x ≡ P y

subst id (cong P eq) : P x → P y -}

2*-even : ∀ m → Even (2 * m)
2*-even zero = zero
2*-even (suc m) with lemma00 (2*-even m)
... | e =  subst Even (sym (proj₁ S.distrib 2 1 m)) e

lemma1 : ∀ m n → 2 * m + 1 ≡ n → Odd n
lemma1 m n eq = 
  subst Odd (trans (S.+-comm 1 (2 * m)) eq) 
  (odd (2*-even m))

pow2-lemma : ∀ m → 2^ suc m ≡ 2 * 2^ m
pow2-lemma m = refl

ggggg : ∀ {n} → Decomposition n → Decomposition (2 * n)
ggggg (inj₁ refl) = inj₁ refl
ggggg (inj₂ (m , k , oddk , eq)) = inj₂ (suc m , k , oddk , 
  trans (S.*-assoc 2 (2^ m) k) 
  (cong (_*_ 2) eq))

ggg : ((n : ℕ) → Decomposition n ⊎ (∃ (λ n' → n' < n × (Decomposition n' → Decomposition n))))
ggg n with div2 n
ggg n | inj₁ p = inj₁ (inj₁ p)
ggg n | inj₂ (m , true , eq , leq) = 
  inj₁ (inj₂ (0 , n , lemma1 m n eq , proj₂ S.+-identity n))
ggg n | inj₂ (m , false , eq , leq) = 
  inj₂ (m , leq , (λ tm → subst Decomposition (trans (sym (proj₂ S.+-identity _)) eq) (ggggg {m} tm)))

decompose : ∀ n → Decomposition n
decompose = wf-rec Decomposition ggg  

suc-inj : ∀ n m → suc n ≡ suc m → n ≡ m
suc-inj .m m refl = refl

2+-inj : ∀ n m → 2 + n ≡ 2 + m → n ≡ m
2+-inj zero zero eq = refl
2+-inj zero (suc m) ()
2+-inj (suc n) zero ()
2+-inj (suc n) (suc m) eq = cong suc (2+-inj _ _ (suc-inj (2 + n) (2 + m) eq))

2*-inj : ∀ n m → 2 * n ≡ 2 * m → n ≡ m
2*-inj zero zero eq = refl
2*-inj zero (suc m) ()
2*-inj (suc n) zero ()
2*-inj (suc n) (suc m) eq = cong suc 
  (2*-inj n m (2+-inj (2 * n) (2 * m) 
   (trans (sym (proj₁ S.distrib 2 1 n)) 
   (trans eq (proj₁ S.distrib 2 1 m)))))

sqrt-lemma : ∀ i j → i ≤ j → ∀ n m → 2^ i * n ≡ 2^ j * m →
  (∃ λ k → n ≡ 2^ k * m × i + k ≡ j)
sqrt-lemma zero j z≤n n m eq = 
  j , (trans (sym (proj₂ S.+-identity _)) eq) , refl
sqrt-lemma (suc i) (suc j) (s≤s leq) n m eq 
  = 
  let 
    (a , b , c) = sqrt-lemma i j leq n m 
      (2*-inj _ _ (trans (sym (S.*-assoc 2 (2^ i) n)) (trans eq (S.*-assoc 2 (2^ j) m))))
  in (a , b , cong suc c)

even+even-even : ∀ m n → Even m → Even n →
  Even (m + n)
even+even-even .0 n zero en = en
even+even-even ._ n (even (odd x)) en = 
  even (odd (even+even-even _ _ x en))

*even-even : ∀ n m → Even n → Even (n * m)
*even-even .zero _ zero = zero
*even-even (suc (suc n)) m (even (odd x)) = 
  subst Even (S.+-assoc m m (n * m)) 
   (even+even-even (m + m) (n * m)
     (subst Even (trans (sym (S.+-assoc m m 0)) (proj₂ S.+-identity (m + m)) ) (2*-even m))
     (*even-even n m x)
    )

square-odd : ∀ k → Odd k → Odd (k * k)
square-odd zero ()
square-odd (suc k) (odd k-even) = 
  odd (subst Even (trans (proj₁ S.distrib k 1 (1 + k))
    (cong (λ z → z + k * suc k) (proj₂ S.*-identity k))) 
      (*even-even k (2 + k) k-even))

module Solver = NatProp.SemiringSolver

foo : ∀ a b → (a * b) * (a * b) ≡ (a * a) * (b * b)
foo = solve 2 (λ a b → 
  (a :* b) :* (a :* b)
   :=
  (a :* a) :* (b :* b))
  refl 
  where
   open Solver

2^*-lemma : ∀ m → 2^ (m * 2) ≡ 2^ m * 2^ m
2^*-lemma zero = refl
2^*-lemma (suc m) = trans (
  trans (sym (S.*-assoc 2 2 (2^ (m * 2)))) 
    (cong (_*_ (2 * 2)) (2^*-lemma m))
  ) (sym (foo 2 (2^ m)))

odd-not-even : ∀ {n} → Odd n → Even n → ⊥
odd-not-even (odd x) (even x₁) = odd-not-even x₁ x

qwe : ∀ n i → Odd n → Twos n i → i ≡ 0
qwe n zero oddd (m , odd-m , eq) = refl
qwe n (suc i) odd-n (m , odd-m , eq) = ⊥-elim (odd-not-even odd-n 
  (subst Even (trans (sym (S.*-assoc 2 (2^ i) m)) eq) (2*-even (2^ i * m)))) 

deco-unique≤ : ∀ n → (i j : ℕ) → i ≤ j → Twos n i → Twos n j → i ≡ j
deco-unique≤ n i j leq (mi , mi-odd , eq-i) (mj , mj-odd , eq-j)
  with sqrt-lemma i j leq mi mj (trans eq-i (sym eq-j))
... | k , eq , kk with qwe mi k mi-odd (mj , mj-odd , sym eq)
deco-unique≤ n i j leq (mi , mi-odd , eq-i) (mj , mj-odd , eq-j) | .0 , eq , kk | refl = 
  trans (sym (proj₂ S.+-identity i)) kk

deco-unique : ∀ m → (i j : ℕ) → Twos m i → Twos m j → i ≡ j
deco-unique m i j twos-i twos-j with DecTotalOrder.total Data.Nat.decTotalOrder i j
deco-unique m i j twos-i twos-j | inj₁ x = deco-unique≤ m i j x twos-i twos-j
deco-unique m i j twos-i twos-j | inj₂ y = sym (deco-unique≤ m j i y twos-j twos-i )

deco2 : ∀ n → n ≡ 0 ⊎ ∃ λ m → Twos (n * n) (2 * m)
deco2 n with decompose n
deco2 n | inj₁ x = inj₁ x
deco2 n | inj₂ (m , k , oddk , eq) = 
  inj₂ (m , k * k , square-odd _ oddk , 
    trans (trans (cong (λ z → z * (k * k)) (trans (cong 2^_ (S.*-comm 2 m)) (2^*-lemma m))) (sym (foo (2^ m) k))) (cong₂ _*_ eq eq))

2*twos : ∀ n {i} → Twos n i → Twos (2 * n) (suc i)
2*twos n {i} (x , odd-x , eq) = x , odd-x , trans (S.*-assoc 2 (2^ i) x) (cong (_*_ 2) eq)

sqrt2 : (n m : ℕ) → n * n ≡ 2 * m * m → m ≡ 0
sqrt2 n m eq with deco2 n | deco2 m
sqrt2 .0 zero eq | inj₁ refl | _ = refl
sqrt2 .0 (suc m) () | inj₁ refl | _
sqrt2 n m eq | inj₂ y | inj₁ x = x
sqrt2 n m eq | inj₂ (i , n2-twos) | inj₂ (j , m2-twos) = ⊥-elim contradiction where
  2mm-twos : Twos (n * n) (suc (2 * j))
  2mm-twos = subst (λ z → Twos z (suc (2 * j))) (trans (sym (S.*-assoc 2 m m)) (sym eq)) (2*twos (m * m) {2 * j} m2-twos)

  bullshit : 2 * i ≡ suc (2 * j)
  bullshit = deco-unique (n * n) (2 * i) (suc (2 * j)) n2-twos 2mm-twos

  contradiction : ⊥
  contradiction = odd-not-even (subst Odd (sym bullshit) (odd (2*-even j))) (2*-even i)