;; 2026-06-25.  heronb.scm.  Based on ss26/solheron.scm, wendler/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)
|#

;; First some general additions to rat.scm

;; RatLeZeroSquare
(set-goal "all a Zero<=a*a")
(cases)
(cases)
;; 3-5
(ng)
(search)
;; 4
(ng)
(search)
;; 5
(ng)
(search)
;; Proof finished.
;; (cp)
(save "RatLeZeroSquare")

;; 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")

;; 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 of general additions to rat.scm

;; We now follow Forster to obtain approximations of sqrt(a) (for 0<a)
;; as a Cauchy sequence with modulus.  The method dates back to Heron.

;; Consider a positive rational number a, which has the form p/q.
;; We define a sequence (a_n)_n approximating sqrt(a) from above

(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)")

(set-totality-goal "RatSqRtR")
(fold-alltotal)
(assume "p")
(fold-alltotal)
(assume "q")
(fold-alltotal)
(ind)
;; 7,8
(use "TotalVar")
;; 8
(assume "n" "IH")
(ng)
(use "RatTimesTotal")
(use "TotalVar")
(use "RatPlusTotal")
(use "IH")
(use "RatTimesTotal")
(use "TotalVar")
(use "RatUDivTotal")
(use "IH")
;; Proof finished.
;; (cp)
(save-totality)

;; By induction an n we prove 0<a_n.

;; RatLtZeroSqRtR
(set-goal "all p,q,n 0<RatSqRtR p q n")
(assume "p" "q")
(ind)
(use "Truth")
;; Step
(assume "n" "IH")
(ng)
(use "RatLtZeroTimes")
(use "Truth")
(use "RatLtZeroPlus")
(use "IH")
(use "RatLtZeroTimes")
(use "Truth")
(use "RatLtZeroUDiv")
(use "IH")
;; Proof finished.
;; (cp)
(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

;; RatSqRtRApproxlbAux
(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")
(simp "<-" "RatTimesAssoc")
(simp (pf "(b+c)*(1#2)=(1#2)*(b+c)"))
(use "Truth")
(use "RatTimesComm")
(use "Truth")
;; 4
(simp "<-" "RatTimesAssoc")
(simp "<-" "RatTimesAssoc")
(simprat "RatSqPlus")
(simprat "RatSqMinus")
(use "RatEqvTrans" (pt "(1#4)*(b*b+2*b*c+c*c)+(1#4)*(4* ~(b*c))"))
;; 14,15
(use "Truth")
;; 15
(simprat "<-" "RatTimesPlusDistr")
(use "RatTimesCompat")
(use "Truth")
;; ?^18:b*b+2*b*c+c*c+4* ~(b*c)==b*b+ ~(2*b*c)+c*c
(simp "<-" "RatPlusAssoc")
(simp "<-" "RatPlusAssoc")
(simp "<-" "RatPlusAssoc")
(use "RatPlusCompat")
(use "Truth")
;; ?^23:2*b*c+(c*c+4* ~(b*c))== ~(2*b*c)+c*c
(use "RatEqvTrans" (pt "2*b*c+(4* ~(b*c)+c*c)"))
(use "RatPlusCompat")
(use "Truth")
(simp "RatPlusComm")
(use "Truth")
;; ?^25:2*b*c+(4* ~(b*c)+c*c)== ~(2*b*c)+c*c
(ng)
;; ?^29:2*b*c+ ~(4*b*c)== ~(2*b*c)
(simp "<-" "RatTimesAssoc")
(simp "<-" "RatTimesAssoc")
(simp "<-" "RatTimesUMinusId")
(simp "<-" "RatTimesUMinusId")
(simprat "<-" "RatTimesPlusDistrLeft")
(use "Truth")
;; Proof finished.
;; (cp)
(save "RatSqRtRApproxLbAux")

;; RatSqRtRApproxLb
(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")
(simp "cDef")
(use "RatLtZeroTimes")
(simp "aDef")
(use "Truth")
(use "RatLtZeroUDiv")
(use "0<b")
;; 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")
(use "Truth")
(simp "RatAbsEq")
(use "0<b")
(use "RatLtToLe")
(use "0<b")
(use "RatTimesComm")
;; 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 "RatSqRtRApproxLbAux")
;; all b,c (1#2)*(b+c)*(1#2)*(b+c)+ ~(b*c)==(1#4)*(b+ ~c)*(b+ ~c)

(simprat "RatSqRtRApproxLbAux")
;; ?^59:0<=(1#4)*(b+ ~c)*(b+ ~c)
(simp "<-" "RatTimesAssoc")
(use "RatLeZeroTimes")
(use "Truth")
(use "Truth")
;; Proof finished.
;; (cp)
(save "RatSqRtRApproxLb")

;; RatSqRtRApproxLbCor
(set-goal "all p,q,n (p#q)<=RatSqRtR p q(Succ n)*RatSqRtR p q(Succ n)")
(assume "p" "q" "n")
(simp (pf "(p#q)eqd~ ~(p#q)"))
;; 3,4
(use "RatZeroLePlusToUMinusLe")
(use "RatSqRtRApproxLb")
(ng)
(use "InitEqD")
;; Proof finished.
;; (cp)
(save "RatSqRtRApproxLbCor")

;; This concludes the proof of a<=a_{n+1}^2.

;; Next we prove a_{n+2}<=a_{n+1}.

;; Useful theorems
;; (pp "RatAbsId") all a(0<=a -> abs a=a)
;; (pp "RatLtZeroSqRtR") all p,q,n 0<RatSqRtR p q n
;; (pp "RatLtToLe") all a,b(a<b -> a<=b)
;; (pp "RatSqRtRApproxLbCor")
;;   all p,q,n (p#q)<=RatSqRtR p q(Succ n)*RatSqRtR p q(Succ n)
;; (pp "RatTimesUDivR") all a(0<abs a -> a*RatUDiv a==1)
;; (pp "RatTimes0RewRule") all a a*1=a
;; (pp "RatLtZeroUDiv") all a(0<a -> 0<RatUDiv a)
;; (pp "RatDoubleEqv") all a a+a==2*a

;; To find theorems try for instance
;; (search-about "Rat" "Le" "Mon" "Times")

;; RatSqRtRDecrSucc
(set-goal "all p,q,n RatSqRtR p q(Succ(Succ n))<=RatSqRtR p q(Succ n)")
(assume "p" "q" "n")
(defnc "a" "p#q")
(defnc "a0" "RatSqRtR p q n")
(defnc "a1" "RatSqRtR p q(Succ n)")
(defnc "a2" "RatSqRtR p q(Succ(Succ n))")
(simp "<-" "a1Def")
(simp "<-" "a2Def")
;; ?^32:a2<=a1

(assert "0<abs a1")
(simp "a1Def")
...
;; Assertion proved.
(assume "0<|a1|")

(assert "a<=a1*a1")
(simp "aDef")
(simp "a1Def")
...
;; Assertion proved.
(assume "aBd")

(assert "a2=(1#2)*(a1+a*RatUDiv a1)")
(simp "aDef")
(simp "a1Def")
(simp "a2Def")
(ng #t)
(use "Truth")
;; Assertion proved.
(assume "a2Prop")

(simp "a2Prop")
(drop "a2Def" "a2Prop")

(use "RatLeTrans" (pt "(1#2)*(a1+a1*(a1*RatUDiv a1))"))
...
;; ?^62:a1+a1*a1*RatUDiv a1<=a1+a1
...
;; ?^55:(1#2)*(a1+a1*(a1*RatUDiv a1))<=a1
(simprat "RatTimesUDivR")
...
;; ?^76:(1#2)*(a1+a1)<=a1
(simprat "RatDoubleEqv")
...
;; Proof finished.
;; (cp)
(save "RatSqRtRDecrSucc")