A276321 List of B-spline interpolation matrices M(n,i,j) of orders n >= 1 read by rows (0 <= i < n, 0 <= j < n).

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, -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, -10, 10, -5, 1, 5, -20, 30, -20, 5, 0, -10, 20, 0, -20, 10, 0, 10

Offset

1,7

Comments

B-spline interpolation matrices of orders n >= 1. Each matrix M(i, j) is of dimensions n X n.

Links

T .C. I. Niessink, Matrices for n = 1..31

J. I. Craig, Notes on B-spline development, (spring 2004), 7.

T. C. I. Niessink, C source code, bigint (n >= 1).

T. C. I. Niessink, C source code, int (1 <= n <= 12).

T. C. I. Niessink, C source code, long long (1 <= n <= 20).

Wikipedia, B-spline.

Formula

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)!).

Example

The B-spline interpolation matrices for orders n = 1..4 are:
--
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
--
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:
--
a0 = -1 * p[i] + 1 * p[i+1]
a1 =  1 * p[i] + 0 * p[i+1]
y = (u*a0 + a1) / (n-1)!
--
where u is in the [0, 1) range. And here is another example for order n = 4 (cubic interpolation):
--
a0 = -1 * p[i] +  3 * p[i+1] + -3 * p[i+2] + 1 * p[i+3]
a1 =  3 * p[i] + -6 * p[i+1] +  3 * p[i+2] + 0 * p[i+3]
a2 = -3 * p[i] +  0 * p[i+1] +  3 * p[i+2] + 0 * p[i+3]
a3 =  1 * p[i] +  4 * p[i+1] +  1 * p[i+2] + 0 * p[i+3]
y = (u^3*a0 + u^2*a1 + u*a2 + a3) / (n-1)!
--
You can optimize the algorithm by dividing all numbers in the matrix by (n-1)!, and then rewriting y as:
--
y = u*(u*(u*a0 + a1) + a2) + a3

Prog

(C)
int fact(int n)
{
  int y = 1;
  int i;
  for (i = 1; i <= n; ++i) y *= i;
  return y;
}
int ipow(int n, int p)
{
  int y = 1;
  int i;
  for (i = 0; i < p; ++i) y *= n;
  return y;
}
int C(int n, int k)
{
  return fact(n) / (fact(k) * fact(n - k));
}
int M(int n, int i, int j)
{
  int y = 0;
  int m;
  for (m = j; m < n; ++m)
  {
    y += ipow(n - (m + 1), i) * ipow(-1, m - j) * C(n, m - j);
  }
  return C(n - 1, i) * y;
}
(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));
doM(n) = for (i=0, n-1, for (j=0, n-1, print1(M(n, i, j), ", ")); print());
tabf(nn) = for (n=1, nn, doM(n)); \\ Michel Marcus, Sep 06 2016

Crossrefs

Cf. A007318.

Sequence in context: A174950 A159906 A263444 * A152196 A024375 A025075

Adjacent sequences: A276318 A276319 A276320 * A276322 A276323 A276324

Keyword

sign,tabf

Author

Theo Niessink, Aug 30 2016

Status

approved