{-# LANGUAGE ParallelListComp #-}-- |Reference implementation of B-Splines; very inefficient but \"obviously\"-- correct.moduleMath.Spline.BSpline.Reference(bases,basisFunctions,basisPolynomials,basisPolynomialsAt)whereimportMath.Spline.KnotsimportMath.Polynomial(Poly)importqualifiedMath.PolynomialasPolyindTrue=1indFalse=0bases::(Fractionala,Orda)=>Knotsa->a->[[a]]basesktsx=coxDeBoorinterpinitialktswhereinitial=[ind(t_j<=x&&x<t_jp1)|(t_j,t_jp1)<-knotSpanskts1]interpt_jd0b_nm1_jt_jpnp1d1b_nm1_jp1=(ifd0==0then0else(x-t_j)/d0)*b_nm1_j+(ifd1==0then0else(t_jpnp1-x)/d1)*b_nm1_jp1-- Alternate version constructing table of functions rather than computing-- table of valuesbasisFunctions::(Fractionala,Orda)=>Knotsa->[[a->a]]basisFunctionskts=coxDeBoorinterpinitialktswhereinitial=[\x->ind(t_j<=x&&x<t_jp1)|(t_j,t_jp1)<-knotSpanskts1]interpt_jd0b_nm1_jt_jpnp1d1b_nm1_jp1x=(ifd0==0then0else(x-t_j)/d0)*b_nm1_jx+(ifd1==0then0else(t_jpnp1-x)/d1)*b_nm1_jp1x-- compute all the basis polynomials for a knot vector, ordered by knot span.basisPolynomials::(Fractionala,Orda)=>Knotsa->[[[Polya]]]basisPolynomialskts|isEmptykts=[]|otherwise=[basisPolynomialsAtktskt|kt<-init(distinctKnotskts)]-- compute all the basis polynomials for the knot span containing a given location.basisPolynomialsAt::(Fractionala,Orda)=>Knotsa->a->[[Polya]]basisPolynomialsAtktsx=coxDeBoorinterpinitialktswhereindPolyTrue=Poly.oneindPolyFalse=Poly.zeroinitial=[indPoly(t_j<=x&&x<t_jp1)|(t_j,t_jp1)<-knotSpanskts1]interpt_jd0b_nm1_jt_jpnp1d1b_nm1_jp1=(ifd0==0thenPoly.zeroelse(Poly.x-Poly.constPolyt_j)/d0)*b_nm1_j+(ifd1==0thenPoly.zeroelse(Poly.constPolyt_jpnp1-Poly.x)/d1)*b_nm1_jp1whereinfixl6+,-p+q=Poly.addPolypqp-q=p+(Poly.negatePolyq)infixl7*,/p*q=Poly.multPolypqp/s=Poly.scalePoly(recips)p-- This is a straightforward implementation of the Cox-De Boor recursion scheme-- generalized in a slightly strange way; the initial vector is a parameter -- and the actual computation of the recursion step is a function parameter.-- The purpose is to allow the same recursion to be applied when computing basis-- function values and basis polynomials.coxDeBoorinterpinitialkts=tablewherets=knotsktstable=initial:[[interpt_jd0b_nm1_jt_jpnp1d1b_nm1_jp1|(b_nm1_j,b_nm1_jp1)<-spans1prevBasis|(d0,d1)<-spans1(spanDiffsnts)|(t_j,t_jpnp1)<-spans(n+1)ts]|prevBasis<-takeWhile(not.null)table|n<-[1..]]spans::Int->[a]->[(a,a)]spans=spansWith(,)spanDiffs::Numa=>Int->[a]->[a]spanDiffs=spansWithsubtractspansWithfnts=zipWithfts(dropnts)