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

;;  Input data

(let ((x (f64vector .000000000E+00
		    .628318500E+00
		    .125663700E+01
		    .188495600E+01
		    .251327400E+01
		    .314159300E+01
		    .376991100E+01
		    .439823000E+01
		    .502654800E+01
		    .565486700E+01
		    .628318500E+01))
      (y (f64vector .200000000E+00
		    .387785300E+00
		    .115105700E+01
		    .751056500E+00
		    .787785200E+00
		    -.200000100E+00
		    -.387785300E+00
		    -.115105700E+01
		    -.751056500E+00
		    -.787784900E+00
		    .200000200E+00)))

  ;;  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* ((xpn 100)
	     (dxp (/ (- (f64vector-ref x (- (f64vector-length x) 1)) (f64vector-ref x 0)) xpn))
	     (xp  (list->f64vector (list-tabulate 100 (lambda (i) (* 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)
		      (printf "y0tab(~A) = ~A~%" (f64vector-ref xp i) (f64vector-ref y0tab i))
		      (assert (< (abs (- (f64vector-ref y0tab i) v)) 1e-15)))
		    '(32 65 98)
		    '(0.746843749779158 -0.806157453777544 -0.0795180362289622))


	  )))
    ))