(use srfi-1 srfi-4 bvsp-spline)

(define machine-epsilon
  (let loop ((e 1.0))
     (if (= 1.0 (+ e 1.0))
         (* 2 e)
         (loop (/ e 2)))))

(print "machine-epsilon = " machine-epsilon)

;;  Input data

(let (
      (x (list->f64vector (list-tabulate 10 (lambda (i) i))))
      (y (list->f64vector (list-tabulate 10 (lambda (i) (* i i)))))
      )

  (print "x = " x)
  (print "y = " y)

  ;;  Set the degree and the class of continuity of the spline.
  (let ((n 3) (k 1))
    
    (let-values (((d d2 constr errc diagn)
		  (compute n k x y)))

      (print "errc = " errc)

      (assert (zero? errc))

      (let* ((xn (f64vector-length x))
             (xpn (* 10 xn))
	     (dxp (/ (- (f64vector-ref x (- xn 1)) (f64vector-ref x 0)) xpn))
	     (xp  (list->f64vector (list-tabulate 101 (lambda (i) (+ 1 (* i dxp)))))))
	
	(let-values (((y0tab y1tab y2tab erre) (evaluate n k x y d d2 xp errc)))
	  
	  (print "erre = " erre)
	  
	  (for-each (lambda (i v)
                      (let ((idx (- i 1)))
                        (printf "v = ~A y0tab(~A) = ~A~%" v (f64vector-ref xp idx) (f64vector-ref y0tab idx))
                        (printf "abs (v - y0tab(~A)) = ~A~%" (f64vector-ref xp idx) (abs (- v (f64vector-ref y0tab idx))))
                        (assert (<= (abs (- (f64vector-ref y0tab idx) v)) machine-epsilon))))
		    '(1 51 101)
		    '(1 30.25 100))


	  )))
    ))