|
%I #21 Nov 21 2020 14:09:30
%S 1,-1,1,1,0,1,-2,1,-2,2,0,1,1,0,-1,3,-3,1,3,-6,3,0,-3,0,3,0,1,4,1,0,1,
%T -4,6,-4,1,-4,12,-12,4,0,6,-6,-6,6,0,-4,-12,12,4,0,1,11,11,1,0,-1,5,
%U -10,10,-5,1,5,-20,30,-20,5,0,-10,20,0,-20,10,0,10
%N List of B-spline interpolation matrices M(n,i,j) of orders n >= 1 read by rows (0 <= i < n, 0 <= j < n).
%C B-spline interpolation matrices of orders n >= 1. Each matrix M(i, j) is of dimensions n X n.
%H T .C. I. Niessink, <a href="/A276321/b276321.txt">Matrices for n = 1..31</a>
%H J. I. Craig, <a href="https://web.archive.org/web/20170712165119/http://soliton.ae.gatech.edu/people/jcraig/classes/ae4375/notes/b-splines-04.pdf">Notes on B-spline development</a>, (spring 2004), 7.
%H T. C. I. Niessink, <a href="/A276321/a276321_2.c.txt">C source code, bigint (n >= 1)</a>.
%H T. C. I. Niessink, <a href="/A276321/a276321.c.txt">C source code, int (1 <= n <= 12)</a>.
%H T. C. I. Niessink, <a href="/A276321/a276321_1.c.txt">C source code, long long (1 <= n <= 20)</a>.
%H Wikipedia, <a href="http://en.wikipedia.org/wiki/B-spline">B-spline</a>.
%F M(n, i, j) = C(n-1, i) * Sum_{m=j..n-1} (n-(m+1))^i * (-1)^(m-j) * C(n, m-j), where binomial coefficient C(n, k) = n! / (k!*(n-k)!).
%e The B-spline interpolation matrices for orders n = 1..4 are:
%e --
%e 1
%e --
%e -1 1
%e 1 0
%e --
%e 1 -2 1
%e -2 2 0
%e 1 1 0
%e --
%e -1 3 -3 1
%e 3 -6 3 0
%e -3 0 3 0
%e 1 4 1 0
%e --
%e These values can be used when implementing interpolation using b-splines. This is what the algorithm for order n = 2 (linear interpolation) looks like in pseudocode:
%e --
%e a0 = -1 * p[i] + 1 * p[i+1]
%e a1 = 1 * p[i] + 0 * p[i+1]
%e y = (u*a0 + a1) / (n-1)!
%e --
%e where u is in the [0, 1) range. And here is another example for order n = 4 (cubic interpolation):
%e --
%e a0 = -1 * p[i] + 3 * p[i+1] + -3 * p[i+2] + 1 * p[i+3]
%e a1 = 3 * p[i] + -6 * p[i+1] + 3 * p[i+2] + 0 * p[i+3]
%e a2 = -3 * p[i] + 0 * p[i+1] + 3 * p[i+2] + 0 * p[i+3]
%e a3 = 1 * p[i] + 4 * p[i+1] + 1 * p[i+2] + 0 * p[i+3]
%e y = (u^3*a0 + u^2*a1 + u*a2 + a3) / (n-1)!
%e --
%e You can optimize the algorithm by dividing all numbers in the matrix by (n-1)!, and then rewriting y as:
%e --
%e y = u*(u*(u*a0 + a1) + a2) + a3
%o (C)
%o int fact(int n)
%o {
%o int y = 1;
%o int i;
%o for (i = 1; i <= n; ++i) y *= i;
%o return y;
%o }
%o int ipow(int n, int p)
%o {
%o int y = 1;
%o int i;
%o for (i = 0; i < p; ++i) y *= n;
%o return y;
%o }
%o int C(int n, int k)
%o {
%o return fact(n) / (fact(k) * fact(n - k));
%o }
%o int M(int n, int i, int j)
%o {
%o int y = 0;
%o int m;
%o for (m = j; m < n; ++m)
%o {
%o y += ipow(n - (m + 1), i) * ipow(-1, m - j) * C(n, m - j);
%o }
%o return C(n - 1, i) * y;
%o }
%o (PARI) M(n, i, j) = binomial(n-1, i) * sum(m=j, n-1, (n-(m+1))^i * (-1)^(m-j) * binomial(n, m-j));
%o doM(n) = for (i=0, n-1, for (j=0, n-1, print1(M(n, i, j), ", ")); print());
%o tabf(nn) = for (n=1, nn, doM(n)); \\ _Michel Marcus_, Sep 06 2016
%Y Cf. A007318.
%K sign,tabf
%O 1,7
%A _Theo Niessink_, Aug 30 2016
|
