(declaim (inline 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) 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 (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)))
		    (* (- (aref nodes i) u) (b-spline i (- m 1) u nodes))
		    0.0d0)))
        (declare (type double-float b1 b2))
	(/ (+ b1 b2) (- (aref nodes i) (aref nodes (- i m)))))))

(defparameter number-of-nodes 11)
(defparameter nodes (make-array number-of-nodes :element-type 'double-float :initial-element 0.0d0))
(defparameter min-u 0.0d0)
(defparameter max-u 1.0d0)
(defparameter du (/ (- max-u min-u) (- number-of-nodes 1)))

(let ((u min-u))
  (loop for i from 0 below number-of-nodes
     do
       (format t "u = ~a~%" u)
       (setf (aref nodes i) u)
       (incf u du)))

(defparameter number-of-steps 10000000)

(time
 (let ((u min-u)
       (du (/ (- max-u min-u) (- number-of-steps 1))))
   (declare
    (type fixnum number-of-steps)
    (type double-float u du)
    (type (simple-array double-float) nodes))
   (loop for i fixnum from 0 below number-of-steps
      do
	(b-spline 4 4 u nodes)
	(incf u du))))