; Lagrange Coefficients
; p(x) = sum (i in S) (p(i) *
;    (product (j in S but not i) (x - j))/
;    (product (j in S but not i) (i - j)))

(lambda (rfi zero zer? one fa fn fm fi)
  (lambda (xl x) (map (lambda (i) (let v ((xl xl)(p one)(q one))
   (if (null? xl) (fm p (fi q))
    (let* ((j (car xl))(mj (fn j))(d (fa i mj)))
    (if (zer? d) (v (cdr xl) p q)
    (v (cdr xl) (fm p (fa x mj))(fm q d))))))) xl)))

; Demo
(define numbers (list '() 0 zero? 1 + - * /))
(define oi (apply (fileVal "LagCoef") numbers))
(define ip (cadddr (apply (fileVal "Matrix") numbers)))
(define la '(1 3 4))
(oi la 3) ; => (0 1 0)

(oi '(2 3) 2.5) ; => (0.5 0.5)
(oi '(0 1 2) 1/2) ; => (3/8 3/4 -1/8) ??
(oi '(0 1 2) 1.01) ; => (-0.004950000000000005 0.9999 0.005050000000000004)
(oi '(0 1 2) 3) ; => (1 -3 3)
(define (p x) (+ (* 3.2 x x) (* -2.1 x) 2.2))
(define P (let ((w (map p la))) (lambda (x) (ip w (oi la x)))))
(cons (p 2) (P 2)) ; => (10.8 . 10.8)

(define GF256 ((fileVal "GFpq") 2 8 "p"))
(define LC (apply (fileVal "LagCoef") GF256))
(define la (list #(1 0 0 1 0 1 1 0) #(0 1 0 1 0 0 1 1) #(0 0 1 0 0 0 0 0)))
(LC la #(0 1 0 1 0 0 1 1)) ; => (#0() #8(1 0) #0())