(declaim (inline b-spline))
(declaim (ftype (function (fixnum fixnum double-float (simple-array double-float)) double-float) b-spline))
(defun b-spline (i m u nodes)
  (declare
   (type fixnum i)
   (type fixnum m)
   (type double-float u)
   (type (simple-array double-float 1) nodes)
   (optimize (debug 0) (safety 0) (speed 3)))
  (if (= m 1)
      (/ 1.0d0 (- (aref nodes i) (aref nodes (- i 1))))
      (let ((b1 (if (and (>= u (aref nodes (- i m))) (< u (aref nodes (- i 1))))
		    (* (- u (the double-float (aref nodes (- i m)))) (b-spline (- i 1) (- m 1) u nodes))
		    0.0d0))
	    (b2 (if (and (>= u (aref nodes (- i m -1))) (< u (aref nodes i)))
		    (* (- (the double-float (aref nodes i)) u) (b-spline i (- m 1) u nodes))
		    0.0d0)))
        (declare (type double-float b1 b2))
	(the double-float (/ (+ b1 b2) (- (aref nodes i) (aref nodes (- i m))))))))

(defmacro def-b-spline (name nodes-p number-of-nodes)
  `(let ((nodes ,nodes-p)
	 (b-spline-number-of-nodes ,number-of-nodes))
     (declare
      (type (simple-array double-float 1) nodes)
      (type fixnum b-spline-number-of-nodes)
      (optimize (debug 0) (safety 0) (speed 3)))
     (declaim (inline ,name))
     (declaim (ftype (function (fixnum fixnum double-float) double-float) ,name))
     (defun ,name(i m u)
       (declare
	(type fixnum i)
	(type fixnum m)
	(type double-float u)
	(optimize (debug 0) (safety 0) (speed 3)))
       (if (= m 1)
	   (/ 1.0d0 (- (aref nodes i) (aref nodes (- i 1))))
	   (let ((b1 (if (and (>= u (aref nodes (- i m))) (< u (aref nodes (- i 1))))
			 (* (- u (the double-float (aref nodes (- i m)))) (,name (- i 1) (- m 1) u))
			 0.0d0))
		 (b2 (if (and (>= u (aref nodes (- i m -1))) (< u (aref nodes i)))
			 (* (- (the double-float (aref nodes i)) u) (,name i (- m 1) u))
			 0.0d0)))
	     (declare (type double-float b1 b2))
	     (the double-float (/ (+ b1 b2) (- (aref nodes i) (aref nodes (- i m))))))))))