Como criar uma função spline recursiva em C ++

No momento, estou trabalhando em um método de solução de equações diferenciais chamado colocação de base-spline. Estou com problemas para criar um método para construir um spline de ordem arbitrária, com o relacionamento com o condição inicial e estou tendo problemas até para começar com esse problema, pois é recursivo ou poderia começar a partir do "top" ou do "bottom", e estou executando um tipo geral de bloco de gravadores é claro, onde eu não consigo entender o que eu preciso fazer.B 1 i (x)={

$$B^{k+1}_{i}(x)= \frac{x-x_i}{x_{k+i}-x_i}B^k_i + \frac{x_{k+i+1}-x}{x_{k+i+1}-x_{i+1}}B^k_{i+1}(x)$$ $$B^1_i(x)=\begin{cases} 1 & \; \text{for } \; x_i \leq x < x_{i+1} \\ 0 & \; \text{otherwise} \end{cases}$$

Posso recomendar a consulta do livro NURBS , que parece ser um texto clássico sobre esse assunto. O próprio algoritmo é apresentado na página 72 , está disponível para visualização online.

Eu apóio o livro NURBS e também gostaria de destacar que, embora essa fórmula recursiva seja tipicamente como se expressa a base do spline B (como em um artigo), não é como você implementa as funções da base. Como todas as funções de ordem superior são baseadas em ordem inferior, é possível calcular todas as funções diferentes zero de uma vez e reutilizar as avaliações de ordem inferior. É assim que o algoritmo de Piegl faz isso.$p+1$

Sinceramente, não sei o quão eficiente isso é, mas uma maneira de fazer isso é com os modelos c ++:

A ordem é k, t é a estrutura do nó e x é o valor que você deseja.

template <int k> 
real BSpline(real x, real *t)
{
    if (*t <= x && x < *(t+k))
    {
        real a = (x - *t) / (*(t+k-1) - *t);
        real b = (*(t+k) - x) / (*(t+k) - *(t+1));

        return a * BSpline<k-1>(x, t) + b * BSpline<k-1>(x, (t+1));
    }
    else
        return 0;
};

template <>
real BSpline<1>(real x, real *t)
{
    if (*t <= x && x < *(t+1))
        return 1.;
    else
        return 0.;
};