; Notes: http://integrals.wolfram.com/index.jsp
; Integral 1/√(a+bx+x^2) = log(2√(a+bx+x^2) + b + 2x)
(define a 3.6)
(define b -0.93)
(define (h x) (sqrt (+ a (* b x) (* x x))))
(define (f x) (/ (h x)))
(define (I x) (log (+ (* 2 (h x)) b x x)))
(define (D f) (let ((ep 0.0000001))
   (lambda (x) (/ (- (f (+ x ep)) (f (- x ep))) (* 2 ep)))))
(define g (D I))
; functions f and g should be near.

(define (sq x) (* x x))
(define (p1 x y z) (let* ((r2 (+ (sq x) (sq y) (sq z))))
  (log (/ (+ (sqrt (+ r2 (* -2 x) 1)) (- x) 1) (- (sqrt r2) x)))))

(define (p0 x y z) (/ (sqrt (+ (sq x) (sq y) (sq z)))))

(define ((poi p) x y z) (let ((e 0.00001)) (/ (- (* 6 (p x y z))
   (+ (p (+ x e) y z) (p (- x e) y z)
      (p x (+ y e) z) (p x (- y e) z)
      (p x y (+ z e)) (p x y (- z e)))) (* e e))))
; Two sticks from [0, 1/2] & [1/2, 1]
(define (pv x y z) (+ (p1 (* 2 x) (* 2 y) (* 2 z))
                (p1 (- (* 2 x) 1) (* 2 y) (* 2 z))))
; pv=p1
; Two sticks together from (-1, 0, 0) to (1, 0, 0)
(define (pu x y z) (+ (p1 x y z) (p1 (- x) y z)))
(define (px x y z) (p1 (/ (+ x 1) 2) (/ y 2) (/ z 2)))
(define (p3 L x y z) (let ((r (/ L))) (px (* x r) (* y r) (* z r))))
; px=pu
(define (pd x y z) (log (/
  (+ (sqrt (+ (/ (sq x) 4) (* -0.5 x) 1/4 (/ (sq y) 4) (/ (sq z) 4))) (* -0.5 x) 0.5)
  (+ (sqrt (+ (/ (sq x) 4) (*  0.5 x) 1/4 (/ (sq y) 4) (/ (sq z) 4))) (* -0.5 x) -0.5))))
(define (pd x y z) (log (/
  (- (sqrt (+ (sq (- x 1)) (sq y) (sq z))) (- x 1))
  (- (sqrt (+ (sq (+ x 1)) (sq y) (sq z))) (+ x 1)))))
(define (pd x y z) (let ((f (lambda (x) (- (sqrt (+ (sq x) (sq y) (sq z))) x))))
  (log (/ (f (- x 1)) (f (+ x 1))))))
(define (pd x y z) (let* ((yz (+ (sq y) (sq z)))
   (f (lambda (x) (- (sqrt (+ (sq x) yz)) x)))) (log (/ (f (- x 1)) (f (+ x 1))))))
; px=pu=pd
; The last definition of pd is a bit quicker than the previous, but not as pretty.
(define ((pg u d) x y z) (let* ((yz (+ (sq y) (sq z)))
   (f (lambda (x) (- (sqrt (+ (sq x) yz)) x)))) (log (/ (f (- x d)) (f (- x u))))))
(define ni (let* ((n 100) (del (+ 0.0 (/ n))))
   (lambda (f l p) (let ((dx (* del (- l f))))
    (* dx (let sum ((s 0)(x (+ f (/ dx 2)))(k n))
    (if (zero? k) s (sum (+ s (p x))(+ x dx)(- k 1)))))))))
(ni 2 5 (lambda (x) 2)) ; => 6.0
(ni 0 1 sq) ; => 1/3 (0.33332500000000054)

; following is the field of a brick with corners (±L, ±B, ±D).
; https://arxiv.org/pdf/1206.3857.pdf
(define ((brick L B D) X Y Z) (let ((I (lambda (x y z)
  (let ((r (sqrt (+ (* x x) (* y y) (* z z)))))(+
    (* y z (log (+ x r))) (* -1/2 (* x x) (atan (/ (* y z)(* x r))))
    (* x z (log (+ y r))) (* -1/2 (* y y) (atan (/ (* x z)(* y r))))
    (* x y (log (+ z r))) (* -1/2 (* z z) (atan (/ (* x y)(* z r)))))))))
  (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X))))
    (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r)))
       (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r)))))))

; br1 = brick
(define ((br1 L B D) X Y Z) (let ((I (lambda (x y z)
  (let ((xx (* x x)) (xy (* x y)) (xz (* x z)) (yy (* y y)) (yz (* y z)) (zz (* z z)))
  (let ((r (sqrt (+ xx yy zz))))
    (- (+ (* yz (log (+ x r))) (* xz (log (+ y r))) (* xy (log (+ z r)))) (* 1/2 (+
         (* xx (atan (/ yz (* x r))))
         (* yy (atan (/ xz (* y r))))
         (* zz (atan  (/ xy (* z r))))))))))))
  (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X))))
    (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r)))
       (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r)))))))

; This version works on the brick faces.
(define ((br2 L B D) X Y Z) (let ((I (lambda (x y z)
  (let* ((r (sqrt (+ (* x x) (* y y) (* z z))))
     (J (lambda (x y z) (+ (* y z (log (+ x r)))
         (if (zero? x) 0 (* -1/2 (* x x) (atan (/ (* y z)(* x r)))))))))
     (+ (J x y z) (J y z x) (J z x y))))))
  (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X))))
    (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r)))
       (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r)))))))

; potential at center of square face of long skinny stick: L by B by B.
; good optimization upon inlining.
(define (stke L B) (let ((I (lambda (x y z)
  (let* ((r (sqrt (+ (* x x) (* y y) (* z z))))
     (J (lambda (x y z) (+ (* y z (log (+ x r)))
         (if (zero? x) 0 (* -1/2 (* x x) (atan (/ (* y z)(* x r)))))))))
     (+ (J x y z) (J y z x) (J z x y))))))
  (let ( (p (- B)) (r (* -2 L)))
    (- (- (- (I B B 0) (I B B r)) (- (I B p 0) (I B p r)))
       (- (- (I p B 0) (I p B r)) (- (I p p 0) (I p p r)))))))

; Faster; less pretty:
(define (stke L B) (let* ((B2 (* B B)) (r (* 2 L)) (r2 (* r r))
     (R (sqrt (+ (* 2 B2) r2))) (Br (* B r)))
(* 4 (- (+ (* B2 (log (* B (sqrt 2)))) (* Br (log (+ B R))))
        (+ (* Br (log (- R B))) (* B2 (log (- R r)))
           (* 1/2 r2 (atan (/ B2 (* r R)))) (* B2 (atan (/ r R))))))))

(define (f x) (let ((b (brick 1/2 1/2 1/2)))
  (* 1/2 (+ (b 0 0 (+ 1/2 x)) (b 0 0 (- 1/2 x))))))
(f 0.001) => 1.7928098404628092

; tran moves field P by vector (p q r). New field is ((tran p q r) P)
(define ((tran p q r) P) (lambda (x y z) (P (- x p) (- y q) (- z r))))
; siz scales field linearly by d. Mass unchanged.
(define ((siz d) P) (let ((r (/ d)))
   (lambda (x y z) (* r (P (* r x) (* r y) (* r z))))))
; mm multiplies mass by d.
(define ((mm d) P) (lambda (x y z) (* d (P x y z))))

Two candidate creators of field for triangle.
(define ((p2c xp yp) x y z) (ni 0 1 (lambda (λ) (let ((b (- 1 λ)))
  (((tran (* λ xp) (* λ yp) 0) ((mm b) ((siz b) p1))) x y z)))))
(define ((pqc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ))(r (/ b)))
 ; (((tran (* λ xp) (* λ yp) 0) ((mm b) ((siz b) p1))) X Y Z)
 ; (((mm b) ((siz b) p1)) (- X (* λ xp)) (- Y (* λ yp)) Z)
 ; (* b (((siz b) p1) (- X (* λ xp)) (- Y (* λ yp)) Z))
 ; (* b (* r (p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z))))
   (p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z))
  ))))
(define (ET pxc) (list
; test case: triangle with verts: (0 0 0), (1 0 0), (0 1 0). 
  (let ((pt (pxc 0 1))) (cons (pt 1 0 1) (pt 0 1 1)))
; test case: triangle with verts: (0 0 0), (1 0 0), (1 1 0). 
  (let ((pt (pxc 1 1))) (cons (pt 0 0 1) (pt 1 1 -1)))
; Equilateral triangle.
  (let* ((t (/ (sqrt 3) 2))(te (pxc 1/2 t))) (list (te 1/2 (- (/ t 3) 1/2) 0)
    (te (/ (- 1 t) 2) (+ (/ t 3) 1/4) 0) (te (/ (+ t 1) 2) (+ (/ t 3) 1/4) 0)))
))
(ET p2c); => ((0.3966777995670832 . 0.3966825245760029)
; (0.3966777995670833 . 0.39668252457600284)
; (1.0135429731716272 1.0135867669194305 1.013586766919481))
(ET pqc); => ((0.3966777995670832 . 0.39668252457600284)
; (0.3966777995670833 . 0.39668252457600284)
; (1.0135429731716272 1.0135867669194303 1.013586766919481))

;(p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z))
;(let* ((r2 (+ (sq (* r (- X (* λ xp)))) (sq (* r (- Y (* λ yp)))) (sq (* r Z)))))
; (log (/ (+ (sqrt (+ r2 (* -2 (* r (- X (* λ xp)))) 1)) (- (* r (- X (* λ xp)))) 1) (- (sqrt r2) (* r (- X (* λ xp)))))))
;(let* ((r2 (* (sq r) (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))))
; (log (/ (+ (sqrt (+ r2 (* -2 (* r (- X (* λ xp)))) 1)) (- (* r (- X (* λ xp)))) 1) (- (sqrt r2) (* r (- X (* λ xp)))))))
;(let* ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1) (- (* r (sqrt s2)) (* r (- X (* λ xp)))))))
;(let* ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1) (* r (- (sqrt s2) (- X (* λ xp)))))))
;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (* b (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1)) (* b (* r (- (sqrt s2) (- X (* λ xp))))))))
;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (- X (* λ xp))) b) (* b (* r (- (sqrt s2) (- X (* λ xp))))))))
;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))
;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))
(define ((prc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ)))
 (let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))
  (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))
 ))))
; (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))
; (+ (* λ λ (+ (sq xp) (sq yp))) (* λ (- (+ (* 2 X xp) (* 2 Y yp)))) (+ (sq X) (sq Y) (sq Z)))

(define ((prc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ)) (r (/ b))) ; good
 (let ((s2
  (+ (* λ λ (+ (sq xp) (sq yp))) (* λ (- (+ (* 2 X xp) (* 2 Y yp)))) (+ (sq X) (sq Y) (sq Z)))
 ))
  (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))
 ))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
 (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0)))
  (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
 (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0)))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))))
   (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))))))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))))
 (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0)))
  (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))))))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))
       (n1 (- xp 1)) (n0 (- 1 X)))
  (ni 0 1 (lambda (λ) (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0))
                  (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp))))))))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))
       (n1 (- xp 1)) (n0 (- 1 X)))
  (ni 0 1 (lambda (λ) (- (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0)))
     (log (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp))))))))))
(define ((prc xp yp) X Y Z) ; good
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))
       (n1 (- xp 1)) (n0 (- 1 X)))
  (- (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0)))))
     (ni 0 1 (lambda (λ) (log (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp))))))))))

(define prc (let ((int (lambda (a b c d e) ; good
      (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e))))))))
   (lambda (xp yp) (lambda (X Y Z)
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))
       (n1 (- xp 1)) (n0 (- 1 X)))
  (- 
 ;    (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0)))))
     (int m2 m1 m0 n1 n0)
   ;  (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- (* λ xp) X)))))
     (int l2 l1 l0 xp (- X))
     )))))))

(define prc (let ((int (lambda (a b c d e) ; good
      (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e))))))))
   (lambda (xp yp) (lambda (X Y Z)
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))
       (n1 (- xp 1)) (n0 (- 1 X)))
  (- (int m2 m1 m0 n1 n0) (int l2 l1 l0 xp (- X)))))))))

(define prc (let ((int (lambda (a b c d e) ; good
      (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e))))))))
   (lambda (xp yp) (lambda (X Y Z)
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (- (int (+ l2 1 (* -2 xp)) (+ l1 (* 2 (+ xp X)) -2) (+ l0 1 (* -2 X)) (- xp 1) (- 1 X))
      (int l2 l1 l0 xp (- X))))))))

(define prc (let ((int (lambda (a b c d e) ; good
      (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e))))))))
   (lambda (xp yp) (lambda (X Y Z)
   (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
   (- (int (+ l2 (* -2 xp) 1) (+ l1 (* 2 (+ xp X)) -2) (+ l0 1 (* -2 X)) (- xp 1) (- 1 X))
      (int l2 l1 l0 xp (- X))))))))

; Log(Sqrt(ax^2 + bx + c) + dx + e)
; Log(Sqrt(-2x^2 + 4x +2) + 2x + 1)

(let ((pt (prc 1 1))) (cons (pt 0 0 1) (pt 1 1 -1)))
(ET prc)

(let ((xp 0.48)(yp 0.5)(X 3/2)(Y 0.6)(Z -0.4))
(let ((λ 1/2)(l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z))))
(let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0)))
(list (+ s2 (* b -2 (- X (* λ xp))) (sq b))
  (+ s2 (* b -2 (- X (* λ xp))) (sq (- 1 λ)))
  (+ s2 (* b -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ)))
  (+ (+ (* λ λ l2) (* λ l1) l0) (* (- 1 λ) -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ)))
  (+ (* λ λ l2) (* λ l1) l0 (* (- 1 λ) -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ)))
  (+ (* λ λ l2) (* λ l1) l0 (* -2 (- 1 λ) (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ)))
  (+ (* λ λ l2) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) (+ 1 (* -2 λ) (* λ λ)))
  (+ (* λ λ l2) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) 1 (* -2 λ) (* λ λ))
  (+ (* λ λ (+ l2 1)) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) 1 (* -2 λ))
  (+ (* λ λ (+ l2 1)) (* -2 xp λ λ) (* λ l1)  (* -2 (+ X (* λ (- (+ xp X))) )) 1 (* -2 λ) l0)
  (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ l1) (* -2 (+ X (* λ (- (+ xp X))) )) 1 (* -2 λ) l0)
  (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)))) (* -2 X) 1 (* -2 λ) l0)
  (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)) -2)) (* -2 X) 1 l0)
  (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)) -2)) (+ (* -2 X) 1 l0))
  (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))))
    (+ (* λ λ m2) (* λ m1) m0)))))) ; same