Lagrange Interpolation

(defun make-poly (xys)
  (case (length xys)
    (0 (error "must have some points."))
    (1 (destructuring-bind ((x1 . y1)) xys
       (declare (ignore x1))
       (constantly y1)))
    (2 (destructuring-bind ((x1 . y1)(x2 . y2)) xys
       (lambda(x)
         (+ (/ (* (- x x2) y1)
           (- x1 x2))
        (/ (* (- x x1) y2)
           (- x2 x1))))))
    (3 (destructuring-bind ((x1 . y1)
                (x2 . y2)
                (x3 . y3))
       xys
     (lambda(x)
       (+ (/ (* (- x x2)(- x x3) y1)
         (- x1 x2)(- x1 x3))
          (/ (* (- x x1)(- x x3) y2)
         (- x2 x1)(- x2 x3))
          (/ (* (- x x1)(- x x2) y3)
         (- x3 x1)(- x3 x2))))))
    (otherwise (lambda(x)
     (loop for (xi . yi) in xys
           for i from 1
           sum (* yi (fancy-product xi xys x i 1)))))))

(defun fancy-product (xi points x omit count)
  ;; neutral element for multiplication is 1
  (if (endp points) 1 ;all done
      (if (= count omit)
      ;; skip this one
      (fancy-product xi (cdr points) x omit (+ count 1))
      ;; do all the rest
      (destructuring-bind ((xj . yj) . rest-of-points) points
        (declare (ignore yj))
        (* (/ (- x xj)(- xi xj))
           (fancy-product xi rest-of-points x omit (+ count 1)))))))