;; 2026-06-07.  heron.scm

#|
(load "~/git/minlog/init.scm")
(set! COMMENT-FLAG #f)
(libload "nat.scm")
(libload "list.scm")
(libload "pos.scm")
(libload "int.scm")
(libload "rat.scm")
;; (set! COMMENT-FLAG #t)
|#

(add-program-constant "RatSqRtR" (py "pos=>pos=>nat=>rat"))
(add-computation-rules
 "RatSqRtR p q Zero" "p#q"
 "RatSqRtR p q(Succ n)" "([b]((1#2)*(b+(p#q)*RatUDiv b)))(RatSqRtR p q n)")

;; Useful theorems:

;; (use "TotalVar") proves totality of a term all of whose program
;; constants are known to be total.

;; RatPlusTotal RatTimesTotal RatUDivTotal

(set-totality-goal "RatSqRtR")
(fold-alltotal)
(assume "p")
...
(save-totality)

;; Write a for p#q and a_n for RatSqRtR p q n.  By induction an n we
;; can prove 0<a_n.  In the proof we need some additions to lib/rat.scm

;; In rat.scm

;; RatLtZeroPlus
(set-goal "all a,b(0<a -> 0<b -> 0<a+b)")
(cases)
(cases)
;; 3-5
(assume "p" "q")
(cases)
(cases)
;; 8-10
(assume "p0" "q0" "Useless1" "Useless2")
(ng)
(use "Truth")
;; 9
(ng)
(search)
;; 10
(ng)
(assume "p0" "q0" "Useless" "Absurd")
(use "EfAtom")
(use "Absurd")
;; 4
(ng)
(assume "p" "b" "Absurd" "Useless")
(use "EfAtom")
(use "Absurd")
;; 5
(ng)
(assume "p" "q" "b" "Absurd" "Useless")
(use "EfAtom")
(use "Absurd")
;; Proof finished.
;; (cp)
(save "RatLtZeroPlus")

;; RatLtZeroTimes
(set-goal "all a,b(0<a -> 0<b -> 0<a*b)")
(cases)
(cases)
;; 3-5
(assume "p" "q")
(cases)
(cases)
;; 8-10
(assume "p0" "q0" "Useless1" "Useless2")
(ng)
(use "Truth")
;; 9
(ng)
(search)
;; 10
(ng)
(search)
;; 4
(ng)
(assume "p" "b" "Absurd" "Useless")
(use "EfAtom")
(use "Absurd")
;; 5
(ng)
(assume "p" "q" "b" "Absurd" "Useless")
(use "EfAtom")
(use "Absurd")
;; Proof finished.
;; (cp)
(save "RatLtZeroTimes")

;; RatLtZeroUDiv
(set-goal "all a(0<a -> 0<RatUDiv a)")
(cases)
(cases)
;; 3-5
(ng)
(search)
;; 4
(ng)
(search)
;; 5
(ng)
(search)
;; Proof finished.
;; (cp)
(save "RatLtZeroUDiv")

;; End rat.scm

;; RatLtZeroSqRtR
(set-goal "all p,q,n 0<RatSqRtR p q n")
(assume "p" "q")
(ind)
...
(save "RatLtZeroSqRtR")

;; Next we aim at a<=a_{n+1}^2.  Informal proof:
;; 0<= 1/4(a_n-a/a_n)^2
;;   = 1/4(a_n^2-2a+a^2/a_n^2)
;;   = 1/4(a_n^2+2a+a^2/a_n^2)-a
;;   = a_{n+1}^2-a

;; In the proof again we need some additions to lib/rat.scm

;; In rat.scm

;; RatSqPlus
(set-goal "all a,b (a+b)*(a+b)==a*a+2*a*b+b*b")
(assume "a" "b")
(simprat "RatTimesPlusDistr")
(simprat "RatTimesPlusDistrLeft")
(simprat "RatTimesPlusDistrLeft")
(simprat (pf "2*a*b==a*b+a*b"))
(ng)
(simp "RatTimesComm")
(use "Truth")
(simp "<-" "RatTimesAssoc")
(use "RatEqvSym")
(use "RatDoubleEqv")
;; Proof finished.
;; (cp)
(save "RatSqPlus")

;; RatSqMinus
(set-goal "all a,b (a+ ~b)*(a+ ~b)==a*a+ ~(2*a*b)+b*b")
(assume "a" "b")
(simprat "RatTimesPlusDistr")
(simprat "RatTimesPlusDistrLeft")
(simprat "RatTimesPlusDistrLeft")
(simprat (pf "2*a*b==a*b+a*b"))
(ng)
(simp "RatTimesComm")
(use "Truth")
(simp "<-" "RatTimesAssoc")
(use "RatEqvSym")
(use "RatDoubleEqv")
;; Proof finished.
;; (cp)
(save "RatSqMinus")

;; RatTimesUMinusId
(set-goal "all a,b ~a*b= ~(a*b)")
(assume "a")
(cases)
(cases)
(ng)
(search)
(ng)
(search)
(ng)
(search)
;; Proof finished.
;; (cp)
(save "RatTimesUMinusId")

;; RatAbsEq
(set-goal "all a(0<=a -> abs a=a)")
(cases)
(cases)
(ng)
(search)
(ng)
(search)
(ng)
(search)
;; Proof finished.
;; (cp)
(save "RatAbsEq")

;; RatLeZeroTimes
(set-goal "all a,b(0<=a -> 0<=b -> 0<=a*b)")
(cases)
(cases)
;; 3,4
(assume "p" "q")
(cases)
(cases)
;; 8,9
(assume "p0" "q0" "Useless1" "Useless2")
(ng)
(use "Truth")
;; 9
(ng)
(search)
;; 10
(ng)
(search)
;; 4
(ng)
(assume "p")
(cases)
(cases)
;; 18-20
(ng)
(search)
;; 19
(ng)
(search)
;; 20
(ng)
(search)
;; 5
(ng)
(assume "p" "q" "b" "Absurd" "Useless")
(use "EfAtom")
(use "Absurd")
;; Proof finished.
;; (cp)
(save "RatLeZeroTimes")

;; End rat.scm

;; Equality for rational numbers is written a==b.
(pp (nf (pf "(1#2)=(1#2)")))
;; T
(pp (nf (pf "(1#2)=(2#4)")))
;; F
(pp (nf (pf "(1#2)==(2#4)")))
;; T

;; Useful theorems:

(pp "RatEqvTrans")
;; all a,b,c(a==b -> b==c -> a==c)
(pp "RatTimesCompat")
;; all a,b,c,d(a==b -> c==d -> a*c==b*d)
(pp "RatTimesAssoc")
;; all a,b,c a*(b*c)=a*b*c
(pp "RatTimesPlusDistr")
;; all a,b,c a*(b+c)==a*b+a*c

;; RatSqRtApproxLbAux
(set-goal "all b,c (1#2)*(b+c)*(1#2)*(b+c)+ ~(b*c)==(1#4)*(b+ ~c)*(b+ ~c)")
(assume "b" "c")
(use "RatEqvTrans" (pt "(1#4)*(b+c)*(b+c)+ ~(b*c)"))
;; 3,4
(use "RatTimesCompat")
...
(save "RatSqRtApproxLbAux")

;; RatSqRtApproxLb
(set-goal "all p,q,n 0<=RatSqRtR p q(Succ n)*RatSqRtR p q(Succ n)+ ~(p#q)")
(assume "p" "q" "n")
(defnc "a" "p#q")
(simp "<-" "aDef")
(defnc "b" "RatSqRtR p q n")
(defnc "c" "a*RatUDiv b")

(assert "0<b")
(simp "bDef")
(use "RatLtZeroSqRtR")
;; Assertion proved.
(assume "0<b")

(assert "0<c")
...
;; Assertion proved.
(assume "0<c")

(assert "a==b*c")
(simp "cDef")
(simp "bDef")
(simp "aDef")
(ng)
(simp "<-" "aDef")
(simp "<-" "bDef")
;; ?^44:a==b*a*RatUDiv b
(simp (pf "b*a=a*b"))
(simp "<-" "RatTimesAssoc")
(simprat "RatTimesUDivR") ;simprat needed because we have == and not =
...
;; Assertion proved.
(assume "a==b*c")
;; ?^53:0<=RatSqRtR p q(Succ n)*RatSqRtR p q(Succ n)+ ~a
(ng)
;; ?^54:0<=
;;      (1#2)*(RatSqRtR p q n+(p#q)*RatUDiv(RatSqRtR p q n))*(1#2)*
;;      (RatSqRtR p q n+(p#q)*RatUDiv(RatSqRtR p q n))+ 
;;      ~a
(simp "<-" "bDef")
(simp "<-" "aDef")
(simp "<-" "cDef")
(simprat "a==b*c")
;; ?^58:0<=(1#2)*(b+c)*(1#2)*(b+c)+ ~(b*c)

;; (pp "RatSqRtApproxLbAux")
;; all b,c (1#2)*(b+c)*(1#2)*(b+c)+ ~(b*c)==(1#4)*(b+ ~c)*(b+ ~c)

(simprat "RatSqRtApproxLbAux")
;; ?^59:0<=(1#4)*(b+ ~c)*(b+ ~c)
...
(save "RatSqRtApproxLb")