Although it is certainly possible to translate numerical algorithms into Haskell starting from Fortran/Matlab/C “everything is an array” style (and, using unboxed mutable vectors, performance will generally be not much worse) this is really missing the point about using a functional language. The underlying math is actually much closer to functional than to imperative programming, so the best thing is to start right there. Specifically, the recurrence formula

can be translated almost literally into Haskell, much better than into an imperative language:
baseFuncs :: [Double] -- ^ Knots, \(\{u_i\}_i\)
-> Int -- ^ Index \(i\) at which to evaluate
-> Int -- ^ Spline degree \(p\)
-> Double -- ^ Position \(u\) at which to evaluate
-> Double
baseFuncs us i 0 u
| u >= us!!i, u < us!!(i+1) = 1
| otherwise = 0
baseFuncs us i p u
= (u - us!!i)/(us!!(i+p) - us!!i) * baseFuncs us i (p-1) u
+ (us!!(i+p+1) - u)/(us!!(i+p+1) - us!!(i+1)) * baseFuncs us (i+1) (p-1) u
Unfortunately this will actually not be efficient though, for multiple reasons.
First, lists suck at random-access. A simple fix is to switch to unboxed (but pure) vectors. While we're at is let's wrap them in a newtype because the ui are supposed to be strictly increasing. Talking about types: the direct accesses are unsafe; we could fix this by bringing p
and the number of segments to the type level and only allowing indices i < n-p
, but I won't go into that here.
Also, it's awkward to pass us
and u
around all the way down the recursion, better just bind it once and then use a helper function to go down:
import Data.Vector.Unboxed (Vector, (!))
import qualified Data.Vector.Unboxed as VU
newtype Knots = Knots {getIncreasingKnotsSeq :: Vector Double}
baseFuncs :: Knots -- ^ \(\{u_i\}_i\)
-> Int -- ^ Index \(i\) at which to evaluate
-> Int -- ^ Spline degree \(p\)
-> Double -- ^ Position \(u\) at which to evaluate
-> Double
baseFuncs (Knots us) i₀ p₀ u = go i₀ p₀
where go i 0
| u >= us!i
, i>=VU.length us-1 || u < us!(i+1) = 1
| otherwise = 0
go i p
= (u - us!i)/(us!(i+p) - us!i) * go i (p-1)
+ (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go (i+1) (p-1)
The other thing that's not optimal is that we don't share the lower-level evaluations between neighbouring recursive calls. (The evaluation is effectively spanning a directed graph with p2⁄2 nodes, but we evaluate it as a tree with 2p nodes.) That's a massive inefficiency for large p, but actually quite harmless for the typical low-degree splines.
The way to avoid this inefficiency is to memoise. The C version does this explicitly with the N
array, but – this being Haskell – we can be lazy to save the effort of allocating the correct size by using a generic memoisation library, e.g. memo-trie:
import Data.MemoTrie (memo2)
baseFuncs (Knots us) i₀ p₀ u = go' i₀ p₀
where go i 0
| u >= us!i
, i>=VU.length us || u < us!(i+1) = 1
| otherwise = 0
go i p
= (u - us!i)/(us!(i+p) - us!i) * go' i (p-1)
+ (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go' (i+1) (p-1)
go' = memo2 go
That was the no-brains version (“just memoise the entire domain of go
”). As dfeuer remarks, it is easy enough to explicitly memoise only the regian that actually gets evaluated, and then we can again use an efficient unboxed vector:
baseFuncs (Knots us) i₀ p₀ u = VU.unsafeHead $ gol i₀ p₀
where gol i 0 = VU.generate (p₀+1) $ \j ->
if u >= us!(i+j)
&& (i+j>=VU.length us || u < us!(i+j+1))
then 1 else 0
gol i p = case gol i (p-1) of
res' -> VU.izipWith
(\j l r -> let i' = i+j
in (u - us!i')/(us!(i'+p) - us!i') * l
+ (us!(i'+p+1) - u)/(us!(i'+p+1) - us!(i'+1)) * r)
res' (VU.unsafeTail res')
(I can safely use unsafeHead
and unsafeTail
here, because at each recursion level the zipping reduces the length by 1, so at the top-level I still have p₀ - (p₀-1) = 1
elements left.)
This version should, I think, have the same asymptotics as the C version. With some more small improvements like precomputing the interval lengths and pre-checking that the arguments are in the allowed range so all accesses can be made unsafe
, it is probably very close in performance to the C version.
As – again – dfeuer remarks, it might not even be necessary to use vectors there because I just zip together the result. For this kind of stuff, GHC can be very good at optimising code even when using plain lists. But, I won't investigate the performance any further here.
The test I used to confirm it actually works:

https://gist.github.com/leftaroundabout/4fd6ef8642029607e1b222783b9d1c1e
STUArray
and theST
monad.