; $Id: extraction.scm 2156 2008-01-25 13:25:12Z schimans $

; (load "~/minlog/init.scm")
; (set! DOT-NOTATION #f)
; (set! COMMENT-FLAG #f)
; (libload "nat.scm")
; (libload "numbers.scm")
; (load "real.scm")
; (load "cont.scm")
; (set! COMMENT-FLAG #t)


(time (pp (nt (pt "(IntN 2#9)+(2*((IntP 7#9)-(IntN 2#9)))/4"))))
; "810#2916" 30 ms
(time (pp (nt (pt "(IntN 2#9)+((IntP 7#9)-(IntN 2#9))/2")))) ;"405#1458" 26 ms
(time (pp (nt (pt "((IntN 2#9)+(7#9))/2")))) ;"45#162" 15 ms
(time (pp (nt (pt "((IntN 2#9)+(7#9))")))) ;"45#81" 11 ms

; We deanimate what was left animated in cont.scm

; (deanimate "ApproxSplit")
; (deanimate "IVTAux")
; (deanimate "IVTcds")

; Now we check whether indeed all lets (i.e., cAp/cId) are not unfolded in
; the present rules.

; (display-pconst)

; We now animate all theorems, working from root to leaves

(animate "IVTApprox")
(animate "RealApprox")
(animate "IVTFinal")
(animate "IVTcds")
(animate "IVTAux")
(animate "ApproxSplit")

; (animate "GCDGInd")
; (animate "QRCorGInd")
; (animate "QRPos")
; (animate "QR")
; (animate "PosInt")

; For further speed-up we provide an external version of +

; We now want to view RatPlus as a program constant with external
; code, using add-external-code.  The external code - after evaluation
; and application to tsubst and the arguments - should give either the
; final value (e.g., the numeral for the sum) or else #f, in which
; case the rules are tried next on the arguments.

(define ratplus-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (sum (+ (/ numer1 denom1) (/ numer2 denom2)))
			  (numer (numerator sum))
			  (denom (denominator sum))
			  (numer-term (int-to-int-term numer))
			  (denom-term (make-numeric-term denom))
			  (constr (constr-name-to-constr "RatConstr"))
			  (term (mk-term-in-app-form
				 (make-term-in-const-form constr)
				 numer-term denom-term)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(define ratminus-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (diff (- (/ numer1 denom1) (/ numer2 denom2)))
			  (numer (numerator diff))
			  (denom (denominator diff))
			  (numer-term (int-to-int-term numer))
			  (denom-term (make-numeric-term denom))
			  (constr (constr-name-to-constr "RatConstr"))
			  (term (mk-term-in-app-form
				 (make-term-in-const-form constr)
				 numer-term denom-term)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(define rattimes-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (prod (* (/ numer1 denom1) (/ numer2 denom2)))
			  (numer (numerator prod))
			  (denom (denominator prod))
			  (numer-term (int-to-int-term numer))
			  (denom-term (make-numeric-term denom))
			  (constr (constr-name-to-constr "RatConstr"))
			  (term (mk-term-in-app-form
				 (make-term-in-const-form constr)
				 numer-term denom-term)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(define ratdiv-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (quot (/ (/ numer1 denom1) (/ numer2 denom2)))
			  (numer (numerator quot))
			  (denom (denominator quot))
			  (numer-term (int-to-int-term numer))
			  (denom-term (make-numeric-term denom))
			  (constr (constr-name-to-constr "RatConstr"))
			  (term (mk-term-in-app-form
				 (make-term-in-const-form constr)
				 numer-term denom-term)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(define ratlt-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (res (< (/ numer1 denom1) (/ numer2 denom2)))
			  (const (if res true-const false-const))
			  (term (make-term-in-const-form const)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(define ratle-code
  '(lambda (tsubst objs)
     (let ((val1 (nbe-object-to-value (car objs)))
	   (val2 (nbe-object-to-value (cadr objs))))
       (and (nbe-constr-value? val1) (nbe-constr-value? val2)
	    (let* ((args1 (nbe-constr-value-to-args val1))
		   (args2 (nbe-constr-value-to-args val2))
		   (vals1 (map nbe-object-to-value args1))
		   (vals2 (map nbe-object-to-value args2)))
	      (and (int-numeral-value? (car vals1))
		   (pos-numeral-value? (cadr vals1))
		   (int-numeral-value? (car vals2))
		   (pos-numeral-value? (cadr vals2))
		   (let* ((numer1 (int-numeral-value-to-number (car vals1)))
			  (denom1 (pos-numeral-value-to-number (cadr vals1)))
			  (numer2 (int-numeral-value-to-number (car vals2)))
			  (denom2 (pos-numeral-value-to-number (cadr vals2)))
			  (res (<= (/ numer1 denom1) (/ numer2 denom2)))
			  (const (if res true-const false-const))
			  (term (make-term-in-const-form const)))
		     (nbe-term-to-object
		      term (nbe-make-bindings '() '())))))))))

(add-external-code "RatPlus" ratplus-code)
(add-external-code "RatMinus" ratminus-code)
(add-external-code "RatTimes" rattimes-code)
(add-external-code "RatDiv" ratdiv-code)
(add-external-code "RatLt" ratlt-code)
(add-external-code "RatLe" ratle-code)

; Finally we animate Id, to enable numeric calculations

(animate "Id")

; Here is a list of the used theorems.  Their formulas can be printed
; using (pp "IVTApprox") und their extracted terms via
; (define IVTApprox-neterm
;   (nt (proof-to-extracted-term (theorem-name-to-proof "IVTApprox"))))
; (pp IVTApprox-neterm)

"RealApprox"
"ApproxSplit"
"IVTAux"
"IVTcds"
"IVTFinal"
"IVTApprox"

; Now the crucial test

(define a-sq-minus-two
  (pt "ContConstr 1 2([a0,n1]a0*a0-2)([k]Zero)([k]k+3)"))

(time (pp (nt
   (apply mk-term-in-app-form
	  (list (proof-to-extracted-term (theorem-name-to-proof "IVTApprox"))
		a-sq-minus-two
		(pt "IntN One") ;-1 is the modulus of increase
		(pt "IntZero")  ;1 <= b-a
		(pt "IntZero")  ;b-a <= 1
		(pt "20"))))))
; 17193534846817967675#12157665459056928801
; 767 ms
(exact->inexact (/ 17193534846817967675 12157665459056928801))
1.4142135186002784
(sqrt 2)
1.4142135623730951

; Check accuracy
; (define diff (- 1.4142135186002784 1.4142135623730951))
; diff
; -4.3772816704645834e-8
; (exact->inexact (expt 2 -20))
; 9.5367431640625e-7

(deanimate "Id")

(pp (nt (proof-to-extracted-term (theorem-name-to-proof "IVTApprox"))))

'(
[f0,k1,k2,k3,k4]
 left((cDC rat@@rat)(f0 doml@f0 domr)
      ([n5,cd6]
        [let cd7
          ((2#3)*left cd6+(1#3)*right cd6@(1#3)*left cd6+(2#3)*right cd6)
          [if (0<=
               (f0 approx left cd7
                (f0 uMod(IntS(IntS(IntS(IntS(IntS(IntS(k2+n5+k1))))))))+
                f0 approx right cd7
                (f0 uMod(IntS(IntS(IntS(IntS(IntS(IntS(k2+n5+k1)))))))))/
               2)
           (left cd6@right cd7)
           (left cd7@right cd6)]])
      (IntToNat(2*(k4+k3))))
)

(define sqrt-two-approx
  (nt (apply mk-term-in-app-form
	     (list (proof-to-extracted-term
		    (theorem-name-to-proof "IVTApprox"))
		   a-sq-minus-two
		   (pt "IntN One") ;-1 is the modulus of increase
		   (pt "IntZero")  ;1 <= b-a
		   (pt "IntZero")  ;b-a <= 1
		   ))))

(pp sqrt-two-approx)

'(
[k0]
 left((cDC rat@@rat)(1@2)
      ([n1,cd2]
        [let cd3
          ((2#3)*left cd2+(1#3)*right cd2@(1#3)*left cd2+(2#3)*right cd2)
          [if (0<=(left cd3*left cd3-2+(right cd3*right cd3-2))/2)
           (left cd2@right cd3)
           (left cd3@right cd2)]])
      (IntToNat(2*k0)))
)

(animate "Id")
(pp (nt (make-term-in-app-form sqrt-two-approx (pt "2"))))
; 107#81
(pp (nt (make-term-in-app-form sqrt-two-approx (pt "20"))))
; 17193534846817967675#12157665459056928801
(time (tag (nbe-normalize-term-without-eta
	    (make-term-in-app-form sqrt-two-approx (pt "20")))))
; 550 ms
(deanimate "Id")

; We now translate terms into scheme expressions, for faster
; evaluation (no conversions between internal and external numbers)

(term-to-expr sqrt-two-approx)

'(
(lambda (k0)
  (car (((|cDC| (cons 1 2))
          (lambda (n1)
            (lambda (cd2)
              (let ([cd3 (cons
                           (+ (* 2/3 (car cd2)) (* 1/3 (cdr cd2)))
                           (+ (* 1/3 (car cd2)) (* 2/3 (cdr cd2))))])
                (if (<= 0
                        (/ (+ (- (* (car cd3) (car cd3)) 2)
                              (- (* (cdr cd3) (cdr cd3)) 2))
                           2))
                    (cons (car cd2) (cdr cd3))
                    (cons (car cd3) (cdr cd2)))))))
         (|IntToNat| (* 2 k0)))))
)

(time ((ev (term-to-expr sqrt-two-approx)) 20))
; 8ms
; 1910392699673572643/1350851717672992089

(time ((ev (term-to-expr sqrt-two-approx)) 100))
; 136 ms
41737211713808721950509113461986613702889339109196103625535604673708288858253142530485267574435/29512665430652752148753480226197736314359272517043832886063884637676943433478020332709411004889

(time ((ev (term-to-expr sqrt-two-approx)) 300))
; 1387 ms
2944593304156165436102846247558257490730845085059145775348712785737552429558941055472664500523847993653960875075500489392461830532348741253430285578096821615417701491158209086792184369992128090401780684332213746112424235660933353732326055138766537198666286440104173743582833475176351331/2082141893205326654083779991150902602700941003443642395329656664801323440350862630969568906052114539645303398663539990042118787521457672342793285135263403898153882623763114393917433013110956461871522162788143751759237923280744039682511207437298831530097535001606799426410247097767236889

; Same for "Inv"

(deanimate "IVTApprox")
(deanimate "RealApprox")
(deanimate "IVTFinal")
(deanimate "IVTcds")
(deanimate "IVTAux")
(deanimate "ApproxSplit")

; We now animate all theorems, working from root to leaves

(animate "InvApprox")
(animate "RealApprox")
(animate "Inv")
(animate "IVTcds")
(animate "IVTAux")
(animate "ApproxSplit")

; We also need to animate "AC" "IP" with identities:

(animate "AC" (pt "[alpha1=>alpha2]alpha1=>alpha2"))
(animate "IP" (pt "[alpha]alpha"))

(define inv-approx-eterm
  (nt (proof-to-extracted-term (theorem-name-to-proof "InvApprox"))))
(pp inv-approx-eterm)

'(
[f0,k1,k2,k3,a4,a5,a6,k7]
 left((cDC rat@@rat)(f0 doml@f0 domr)
      ([n8,cd9]
        [let cd10
          ((2#3)*left cd9+(1#3)*right cd9@(1#3)*left cd9+(2#3)*right cd9)
          [if (0<=
               (f0 approx left cd10
                (f0 uMod(IntS(IntS(IntS(IntS(IntS(IntS(k2+n8+k1))))))))-
                a6+
                (f0 approx right cd10
                 (f0 uMod(IntS(IntS(IntS(IntS(IntS(IntS(k2+n8+k1))))))))-
                 a6))/
               2)
           (left cd9@right cd10)
           (left cd10@right cd9)]])
      (IntToNat(2*f0 uModCont(IntS(IntS(IntS(IntS(k7+k1)))))+k3+k3)))
)

(define sq (pt "ContConstr 1 2([a0,n1]a0*a0)([k]Zero)([k]k+3)"))

(define inv-sq-approx
  (nt (apply mk-term-in-app-form
	     (list (proof-to-extracted-term
		    (theorem-name-to-proof "InvApprox"))
		   sq ;continuous function to be inverted
  		   (pt "IntN One") ;uniform lower bound on the slope
		   (pt "IntZero") (pt "IntZero") ;bounds for b-a
		   (pt "1") (pt "4") ;interval in range
		   ))))

(pp inv-sq-approx)

'(
[a0,k1]
 left((cDC rat@@rat)(1@2)
      ([n2,cd3]
        [let cd4
          ((2#3)*left cd3+(1#3)*right cd3@(1#3)*left cd3+(2#3)*right cd3)
          [if (0<=(left cd4*left cd4-a0+(right cd4*right cd4-a0))/2)
           (left cd3@right cd4)
           (left cd4@right cd3)]])
      (IntToNat(2*IntS(IntS(IntS(IntS(IntS(IntS k1))))))))
)

; Finally we animate Id, to enable numeric calculations

(animate "Id")

(time
(pp (nbe-normalize-term-without-eta
      (mk-term-in-app-form
       inv-sq-approx
       (pt "3") ;argument of inverted function
       (pt "20") ;error bound (number of binary digits)
       ))))

; 3730307366945298869534434#2153693963075557766310747 in 803 ms

; (exact->inexact (/ 3730307366945298869534434 2153693963075557766310747))
; 1.7320508070785863
; (sqrt 3)
; 1.7320508075688772
; Difference at the 10th decimal digit

; We now translate terms into scheme expressions, for faster
; evaluation (no conversions between internal and external numbers)

(term-to-expr inv-sq-approx)

'(
(lambda (a0)
  (lambda (k1)
    (car (((|cDC| (cons 1 2))
            (lambda (n2)
              (lambda (cd3)
                (let ([cd4 (cons
                             (+ (* 2/3 (car cd3)) (* 1/3 (cdr cd3)))
                             (+ (* 1/3 (car cd3)) (* 2/3 (cdr cd3))))])
                  (if (<= 0
                          (/ (+ (- (* (car cd4) (car cd4)) a0)
                                (- (* (cdr cd4) (cdr cd4)) a0))
                             2))
                      (cons (car cd3) (cdr cd4))
                      (cons (car cd4) (cdr cd3)))))))
           (|IntToNat|
             (* 2 (+ (+ (+ (+ (+ (+ k1 1) 1) 1) 1) 1) 1)))))))
)

(time (((ev (term-to-expr inv-sq-approx)) 3) 20))
; 13 ms

(time (((ev (term-to-expr inv-sq-approx)) 3) 100))
; 134 ms

(time (((ev (term-to-expr inv-sq-approx)) 3) 200))
; 560 ms