(define (<_e e1 e2) #t)

(define (some p l)
 (and (not (null? l)) (or (p (car l)) (some p (cdr l)))))

(define (find-if p l)
 (let loop ((l l))
  (cond ((null? l) #f)
	((p (car l)) (car l))
	(else (loop (cdr l))))))

(define (remove-if p l)
 (let loop ((l l) (c '()))
  (cond ((null? l) (reverse c))
	((p (car l)) (loop (cdr l) c))
	(else (loop (cdr l) (cons (car l) c))))))

(define (removeq x l)
 (let loop ((l l) (c '()))
  (cond ((null? l) (reverse c))
	((eq? x (car l)) (loop (cdr l) c))
	(else (loop (cdr l) (cons (car l) c))))))

(define terms
 (let ((pair? pair?))
  (lambda (p)
   (if (and (pair? p) (eq? (car p) 'dual-number))
       (cadr p)
       (list (cons '() p))))))

(define (terms->dual-number terms)
 (cond ((null? terms) 0)
       ((and (null? (cdr terms)) (null? (car (car terms))))
	(cdr (car terms)))
       (else (list 'dual-number terms))))

(define (dual-number? p)
 (some (lambda (term) (not (null? (car term)))) (terms p)))

(define (dual-number e x x-prime)
 (terms->dual-number
  (append (terms x)
	  (map (lambda (term) (cons (cons e (car term)) (cdr term)))
	       (terms x-prime)))))

(define (epsilon p)
 (car (car (find-if (lambda (term) (not (null? (car term)))) (terms p)))))

(define (primal e p)
 (terms->dual-number
  ;; memq
  (remove-if (lambda (term) (memq e (car term))) (terms p))))

(define (perturbation e p)
 (terms->dual-number
  ;; removeq
  (map (lambda (term) (cons (removeq e (car term)) (cdr term)))
       ;; memq
       (remove-if (lambda (term) (not (memq e (car term)))) (terms p)))))

(define (generate-epsilon) (cons #f #f))

(define (lift-real->real f df/dx)
 (letrec ((self (lambda (p)
		 (if (dual-number? p)
		     (let ((e (epsilon p)))
		      (dual-number
		       e
		       (self (primal e p))
		       (* (df/dx (primal e p)) (perturbation e p))))
		     (f p)))))
  self))

(define (lift-real*real->real f df/dx1 df/dx2)
 (letrec ((self
	   (lambda (p1 p2)
	    (if (or (dual-number? p1)
		    (dual-number? p2))
		(let ((e (if (or (not (dual-number? p1))
				 (and (dual-number? p2)
				      (<_e (epsilon p1) (epsilon p2))))
			     (epsilon p2)
			     (epsilon p1))))
		 (dual-number
		  e
		  (self (primal e p1) (primal e p2))
		  (+ (* (df/dx1 (primal e p1) (primal e p2))
			(perturbation e p1))
		     (* (df/dx2 (primal e p1) (primal e p2))
			(perturbation e p2)))))
		(f p1 p2)))))
  self))

(define (primal* p)
 (if (dual-number? p) (primal* (primal (epsilon p) p)) p))

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

(define pair?
 (let ((pair? pair?))
  (lambda (x) (and (pair? x) (not (dual-number? 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))))

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

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

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

(define (derivative f)
 (lambda (x)
  (let ((e (generate-epsilon)))
   (perturbation e (f (dual-number e x 1))))))