;;; Figure 1

(define <_e <)

(define =_e =)

(define linear-term?
 (let ((pair? pair?)) (lambda (p) (and (pair? p) (eq? (car p) 'linear-term)))))

(define (invariant1? e p)
 (or (not (linear-term? p))
     (and (<_e (epsilon p) e) (invariant1? (epsilon p) (caddr p)))))

(define (invariant2 e p)
 (if (and (linear-term? p) (<_e e (epsilon p)))
     (error "Invariant (2) violated")
     p))

(define-syntax linear-term
 (syntax-rules ()
  ;; Can't check whether x-prime is zero because that would break laziness.
  ;; We want the following invariants:
  ;;  1. x does not contain any (nested) linear terms over e'>=e. By (2) we
  ;;     only need to check the left spine.
  ;;  2. x-prime (after forcing) does not contain any (nested) linear terms
  ;;     over e'>e.
  ((linear-term e x x-prime)
   (if (invariant1? e x)
       (list 'linear-term e x (delay (invariant2 e x-prime)))
       (error "Invariant (1) violated")))))

(define epsilon cadr)

(define (r e p)
 (cond
  ;; Equation (13)
  ((not (linear-term? p)) p)
  ;; Equation (14)
  ;; This is an optimization. From the invariants, p cannot contain e.
  ((<_e (epsilon p) e) p)
  ;; Equation (15)
  ;; We optimize here. This would nominally be (r e (caddr p)). But from
  ;; invariant (1), (caddr p) cannot contain e. This just saves one level of
  ;; recursion since this would be caught by equations (13) or (14) on the next
  ;; recursive call.
  ((=_e (epsilon p) e) (caddr p))
  ;; Equation (16)
  ;; The linear-term invariants are upheld by the following argument.
  ;; (caddr p) cannot contain e'>=(epsilon p). Ditto for (r e (caddr p)). Thus
  ;; invariant (1) is upheld. (force (cadddr p)) cannot contain
  ;; e'>(epsilon p). Ditto for (r e (force (cadddr p))). Thus invariant (2) is
  ;; upheld.
  (else (linear-term (epsilon p) (r e (caddr p)) (r e (force (cadddr p)))))))

(define (q e p)
 (cond
  ;; Equation (17)
  ((not (linear-term? p)) 0)
  ;; Equation (18)
  ;; This is an optimization. From the invariants, p cannot contain e.
  ((<_e (epsilon p) e) 0)
  ;; Equation (19)
  ;; We optimize here. This would nominally be
  ;; (+ (q e (caddr p)) (force (cadddr p))) which must use the lifted +. But
  ;; from invariant (1), (caddr p) cannot contain e. Thus (q e (caddr p))=0.
  ((=_e (epsilon p) e) (force (cadddr p)))
  ;; Equation (20)
  ;; The linear-term invariants are upheld by the following argument.
  ;; (caddr p) cannot contain e'>=(epsilon p). Ditto for (q e (caddr p)). Thus
  ;; invariant (1) is upheld. (force (cadddr p)) cannot contain
  ;; e'>(epsilon p). Ditto for (q e (force (cadddr p))). Thus invariant (2) is
  ;; upheld.
  (else (linear-term (epsilon p) (q e (caddr p)) (q e (force (cadddr p)))))))

;;; This can use the non-lifted +.
(define generate-epsilon (let ((e 0)) (lambda () (set! e (+ e 1)) e)))

(define (univariate-e e i p)
 ;; This can use the non-lifted zero? and -.
 ;; This must use the lifted *.
 (cond
  ;; Equation (5)
  ((zero? i) (r e p))
  ;; Equation (6)
  (else (* i (univariate-e e (- i 1) (q e p))))))

(define (multivariate-e e i p)
 (cond
  ;; Equation (10)
  ((null? i) p)
  ;; Equation (11)
  (else (multivariate-e (cdr e) (cdr i) (univariate-e (car e) (car i) p)))))

;;; Figure 2

(define (c e i p)
 (cond
  ;; Equation (24)
  ;; This can use the non-lifted / and +.
  ((not (linear-term? p)) (/ p (+ i 1)))
  ;; Equation (25)
  ;; The linear-term invariants are upheld by the following argument.
  ;; (r (epsilon p) p) cannot contain e'>=(epsilon p). Ditto for
  ;; (c e i (r (epsilon p) p)). Thus invariant (1) is upheld. (q (epsilon p) p)
  ;; cannot contain e'>(epsilon p). Ditto for (c e (+ i 1) (q (epsilon p) p)).
  ;; Thus invariant (2) is upheld.
  ((=_e (epsilon p) e)
   (linear-term
    ;; This can use the non-lifted +.
    (epsilon p) (c e i (r (epsilon p) p)) (c e (+ i 1) (q (epsilon p) p))))
  ;; Equation (26)
  ;; The linear-term invariants are upheld by the following argument.
  ;; (r (epsilon p) p) cannot contain e'>=(epsilon p). Ditto for
  ;; (c e i (r (epsilon p) p)). Thus invariant (1) is upheld. (q (epsilon p) p)
  ;; cannot contain e'>(epsilon p). Ditto for (c e i (q (epsilon p) p)). Thus
  ;; invariant (2) is upheld.
  (else (linear-term
	 (epsilon p) (c e i (r (epsilon p) p)) (c e i (q (epsilon p) p))))))

(define (rename e1 e2 p)
 (cond
  ;; Equation (27)
  ((=_e e1 e2) p)
  ;; Equation (28)
  ((not (linear-term? p)) p)
  ;; Equation (29)
  ((<_e (epsilon p) e1) p)
  ;; Equation (30)
  ((=_e (epsilon p) e1)
   ;; This must use the lifted +.
   ;; needs work: To prove that the linear-term invariants are upheld.
   (linear-term
    e2 (r e2 (r e1 p)) (+ (q e2 (r e1 p)) (rename e1 e2 (q e1 p)))))
  (else (error "This case should never occur in this program."))))

(define (lift-real->real f df/dx)
 (letrec ((self
	   (lambda (p)
	    (cond
	     ;; Equation (31)
	     ((linear-term? p)
	      (let ((e (epsilon p)))
	       ;; This must use the lifted * and df/dx which precludes eta
	       ;; reduction.
	       ;; The linear-term invariants are upheld by the following
	       ;; argument. (r e p) cannot contain contain e'>=e. Ditto for
	       ;; (self (r e p)), because lifted computation does not introduce
	       ;; epsilons. Thus invariant (1) is upheld.
	       ;; needs work: To prove that invariant (2) is upheld.
	       ;; This upholds the induction that lifted computation does not
	       ;; introduce epsilons.
	       (linear-term
		e
		(self (r e p))
		(* (let ((e-prime (generate-epsilon)))
		    (rename
		     e-prime
		     e
		     (c e-prime
			0
			(df/dx (linear-term e-prime (r e p) (q e p))))))
		   (q e p)))))
	     (else (f p))))))
  self))

(define (lift-real*real->real f df/dx1 df/dx2)
 (letrec ((self
	   (lambda (p1 p2)
	    (cond
	     ;; Equation (32)
	     ((and
	       (linear-term? p1)
	       (or (not (linear-term? p2)) (<_e (epsilon p2) (epsilon p1))))
	      (let ((e1 (epsilon p1)))
	       ;; This must use the lifted * and df/dx1 which precludes eta
	       ;; reduction.
	       ;; The linear-term invariants are upheld by the following
	       ;; argument. p2 cannot contain e'>=e1. Ditto for (r e1 p1).
	       ;; Ditto for (self (r e1 p1) p2), because lifted computaion
	       ;; does not introduce epsilons. Thus invariant (1) is upheld.
	       ;; needs work: To prove that invariant (2) is upheld.
	       ;; This upholds the induction that lifted computation does not
	       ;; introduce epsilons.
	       (linear-term
		e1
		(self (r e1 p1) p2)
		(* (let ((e-prime (generate-epsilon)))
		    (rename
		     e-prime
		     e1
		     (c e-prime
			0
			(df/dx1
			 (linear-term e-prime (r e1 p1) (q e1 p1)) p2))))
		   (q e1 p1)))))
	     ;; Equation (33)
	     ((and
	       (linear-term? p2)
	       (or (not (linear-term? p1)) (<_e (epsilon p1) (epsilon p2))))
	      (let ((e2 (epsilon p2)))
	       ;; This must use the lifted * and df/dx2 which precludes eta
	       ;; reduction.
	       ;; The linear-term invariants are upheld by the following
	       ;; argument. p1 cannot contain e'>=e2. Ditto for (r e2 p2).
	       ;; Ditto for (self p1 (r e2 p2)), because lifted computaion
	       ;; does not introduce epsilons. Thus invariant (1) is upheld.
	       ;; needs work: To prove that invariant (2) is upheld.
	       ;; This upholds the induction that lifted computation does not
	       ;; introduce epsilons.
	       (linear-term
		e2
		(self p1 (r e2 p2))
		(* (let ((e-prime (generate-epsilon)))
		    (rename
		     e-prime
		     e2
		     (c e-prime
			0
			(df/dx2
			 p1 (linear-term e-prime (r e2 p2) (q e2 p2))))))
		   (q e2 p2)))))
	     ;; Equation (34)
	     ((and (linear-term? p1)
		   (linear-term? p2)
		   (=_e (epsilon p1) (epsilon p2)))
	      (let ((e (epsilon p1)) (e-prime (generate-epsilon)))
	       ;; This upholds the induction that lifted computation does not
	       ;; introduce epsilons.
	       (rename e-prime e (self p1 (rename e e-prime p2)))))
	     (else (f p1 p2))))))
  self))

(define (r* p) (if (linear-term? p) (r* (r (epsilon p) p)) p))

(define (lift-real^n->boolean f) (lambda ps (apply f (map r* ps))))

;;; Figure 3

(define pair?
 (let ((pair? pair?)) (lambda (x) (and (pair? x) (not (linear-term? x))))))

(define + (lift-real*real->real + (lambda (x1 x2) 1) (lambda (x1 x2) 1)))

(define - (lift-real*real->real - (lambda (x1 x2) 1) (lambda (x1 x2) -1)))

(define * (lift-real*real->real * (lambda (x1 x2) x2) (lambda (x1 x2) x1)))

(define /
 (lift-real*real->real
  / (lambda (x1 x2) (/ 1 x2)) (lambda (x1 x2) (- 0 (/ x1 (* x2 x2))))))

(define sqrt (lift-real->real sqrt (lambda (x) (/ 1 (* 2 (sqrt x))))))

(define exp (lift-real->real exp (lambda (x) (exp x))))	;Can't be eta reduced.

(define log (lift-real->real log (lambda (x) (/ 1 x))))

(define sin (lift-real->real sin (lambda (x) (cos x))))	;Can't be eta reduced.

(define cos (lift-real->real cos (lambda (x) (- 0 (sin x)))))

(define atan (lift-real*real->real
              atan
              (lambda (x1 x2) (/ (- 0 x2) (+ (* x1 x1) (* x2 x2))))
              (lambda (x1 x2) (/ x1 (+ (* x1 x1) (* x2 x2))))))

(define = (lift-real^n->boolean =))

(define < (lift-real^n->boolean <))

(define > (lift-real^n->boolean >))

(define <= (lift-real^n->boolean <=))

(define >= (lift-real^n->boolean >=))

(define zero? (lift-real^n->boolean zero?))

(define positive? (lift-real^n->boolean positive?))

(define negative? (lift-real^n->boolean negative?))

(define real? (lift-real^n->boolean real?))

;;; Figure 4

;;; Equation (2)
(define (derivative f)
 (lambda (c)
  ;; The linear-term invariants are trivially upheld.
  (let ((e (generate-epsilon))) (univariate-e e 1 (f (linear-term e c 1))))))

(define (ith-derivative-by-repetition i f)
 ;; This can use the non-lifted zero? and -.
 (cond
  ;; Equation (3)
  ((zero? i) f)
  ;; Equation (4)
  (else (ith-derivative-by-repetition (- i 1) (derivative f)))))

;;; Equation (7)
(define (ith-derivative-by-tower i f)
 (lambda (c)
  ;; The linear-term invariants are trivially upheld.
  (let ((e (generate-epsilon))) (univariate-e e i (f (linear-term e c 1))))))

(define (position-of-nonzero i)
 (cond ((null? i) #f)
       ;; This can use the non-lifted zero?.
       ((zero? (car i))
	(let ((position (position-of-nonzero (cdr i))))
	 ;; This can use the non-lifted +.
	 (if position (+ position 1) #f)))
       (else 0)))

(define (decrement-lth i l)
 ;; This can use the non-lifted zero?.
 (if (zero? l)
     ;; This can use the non-lifted -.
     (cons (- (car i) 1) (cdr i))
     ;; This can use the non-lifted -.
     (cons (car i) (decrement-lth (cdr i) (- l 1)))))

(define (list-replace-lth c l u)
 ;; This can use the non-lifted zero?.
 (if (zero? l)
     (cons u (cdr c))
     ;; This can use the non-lifted -.
     (cons (car c) (list-replace-lth (cdr c) (- l 1) u))))

(define (partial-derivative-by-repetition i f)
 ;; This can use the non-lifted zero?.
 (let ((l (position-of-nonzero i)))
  (cond
   ;; Equation (8)
   ((not l) f)
   ;; Equation (9)
   (else (partial-derivative-by-repetition
	  (decrement-lth i l)
	  (lambda (c)
	   ((derivative (lambda (u) (f (list-replace-lth c l u))))
	    (list-ref c l))))))))

;;; Equation (12)
(define (partial-derivative-by-tower i f)
 (lambda (c)
  (let ((e (map (lambda (cl) (generate-epsilon)) c)))
   ;; The linear-term invariants are trivially upheld.
   (multivariate-e e i (f (map (lambda (el cl) (linear-term el cl 1)) e c))))))

(define (f x) (+ (* x (* x (* x x))) (* 2 (* x (* x x)))))

(define (test1)
 ;; (135 162 144 84 24)
 (list ((ith-derivative-by-repetition 0 f) 3)
       ((ith-derivative-by-repetition 1 f) 3)
       ((ith-derivative-by-repetition 2 f) 3)
       ((ith-derivative-by-repetition 3 f) 3)
       ((ith-derivative-by-repetition 4 f) 3)))

(define (test2)
 ;; (135 162 144 84 24)
 (list ((ith-derivative-by-tower 0 f) 3)
       ((ith-derivative-by-tower 1 f) 3)
       ((ith-derivative-by-tower 2 f) 3)
       ((ith-derivative-by-tower 3 f) 3)
       ((ith-derivative-by-tower 4 f) 3)))

(define (g xy)
 (let ((x (list-ref xy 0)) (y (list-ref xy 1)))
  (+ (* x (* x (* x y))) (* x (* x (* y y))))))

(define (test3)
 ;; (60 72 32 54 36 8 18 24 8 6 4)
 (list ((partial-derivative-by-repetition (list 0 0) g) (list 2 3))
       ((partial-derivative-by-repetition (list 1 0) g) (list 2 3))
       ((partial-derivative-by-repetition (list 0 1) g) (list 2 3))
       ((partial-derivative-by-repetition (list 2 0) g) (list 2 3))
       ((partial-derivative-by-repetition (list 1 1) g) (list 2 3))
       ((partial-derivative-by-repetition (list 0 2) g) (list 2 3))
       ((partial-derivative-by-repetition (list 3 0) g) (list 2 3))
       ((partial-derivative-by-repetition (list 2 1) g) (list 2 3))
       ((partial-derivative-by-repetition (list 1 2) g) (list 2 3))
       ((partial-derivative-by-repetition (list 3 1) g) (list 2 3))
       ((partial-derivative-by-repetition (list 2 2) g) (list 2 3))))

(define (test4)
 ;; (60 72 32 54 36 8 18 24 8 6 4)
 (list ((partial-derivative-by-tower (list 0 0) g) (list 2 3))
       ((partial-derivative-by-tower (list 1 0) g) (list 2 3))
       ((partial-derivative-by-tower (list 0 1) g) (list 2 3))
       ((partial-derivative-by-tower (list 2 0) g) (list 2 3))
       ((partial-derivative-by-tower (list 1 1) g) (list 2 3))
       ((partial-derivative-by-tower (list 0 2) g) (list 2 3))
       ((partial-derivative-by-tower (list 3 0) g) (list 2 3))
       ((partial-derivative-by-tower (list 2 1) g) (list 2 3))
       ((partial-derivative-by-tower (list 1 2) g) (list 2 3))
       ((partial-derivative-by-tower (list 3 1) g) (list 2 3))
       ((partial-derivative-by-tower (list 2 2) g) (list 2 3))))

(define (test5)
 ;; (1 1 1 1 1)
 (list ((ith-derivative-by-repetition 0 exp) 0)
       ((ith-derivative-by-repetition 1 exp) 0)
       ((ith-derivative-by-repetition 2 exp) 0)
       ((ith-derivative-by-repetition 3 exp) 0)
       ((ith-derivative-by-repetition 4 exp) 0)))

(define (test6)
 ;; (1 1 1 1 1)
 (list ((ith-derivative-by-tower 0 exp) 0)
       ((ith-derivative-by-tower 1 exp) 0)
       ((ith-derivative-by-tower 2 exp) 0)
       ((ith-derivative-by-tower 3 exp) 0)
       ((ith-derivative-by-tower 4 exp) 0)))

(define (test7)
 ;; (0 1 0 -1 0 1 0 -1)
 (list ((ith-derivative-by-repetition 0 sin) 0)
       ((ith-derivative-by-repetition 1 sin) 0)
       ((ith-derivative-by-repetition 2 sin) 0)
       ((ith-derivative-by-repetition 3 sin) 0)
       ((ith-derivative-by-repetition 4 sin) 0)
       ((ith-derivative-by-repetition 5 sin) 0)
       ((ith-derivative-by-repetition 6 sin) 0)
       ((ith-derivative-by-repetition 7 sin) 0)))

(define (test8)
 ;; (0 1 0 -1 0 1 0 -1)
 (list ((ith-derivative-by-tower 0 sin) 0)
       ((ith-derivative-by-tower 1 sin) 0)
       ((ith-derivative-by-tower 2 sin) 0)
       ((ith-derivative-by-tower 3 sin) 0)
       ((ith-derivative-by-tower 4 sin) 0)
       ((ith-derivative-by-tower 5 sin) 0)
       ((ith-derivative-by-tower 6 sin) 0)
       ((ith-derivative-by-tower 7 sin) 0)))

(define (sqr x) (* x x))

(define (test9)
 ;; ((0 0) (0 0) (0 0))
 (list
  (list ((ith-derivative-by-tower 3 (lambda (x) (sqrt (* x x)))) 1)
	((ith-derivative-by-repetition 3 (lambda (x) (sqrt (* x x)))) 1))
  (list ((ith-derivative-by-tower 2 (lambda (x) (/ x x))) 1)
	((ith-derivative-by-repetition 2 (lambda (x) (/ x x))) 1))
  (list ((ith-derivative-by-tower 4 (lambda (x) (/ x x))) 1)
	((ith-derivative-by-repetition 4 (lambda (x) (/ x x))) 1))
  (list
   ((ith-derivative-by-tower 4 (lambda (x) (+ (sqr (sin x)) (sqr (cos x))))) 0)
   ((ith-derivative-by-repetition
     4 (lambda (x) (+ (sqr (sin x)) (sqr (cos x)))))
    0))))

;;; Local Variables:
;;; lisp-body-indent: 1
;;; End: