login
The OEIS is supported by the many generous donors to the OEIS Foundation.

 

Logo
Hints
(Greetings from The On-Line Encyclopedia of Integer Sequences!)
A276321 List of B-spline interpolation matrices M(n,i,j) of orders n >= 1 read by rows (0 <= i < n, 0 <= j < n). 2
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 (list; graph; refs; listen; history; text; internal format)
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

Lookup | Welcome | Wiki | Register | Music | Plot 2 | Demos | Index | Browse | More | WebCam
Contribute new seq. or comment | Format | Style Sheet | Transforms | Superseeker | Recents
The OEIS Community | Maintained by The OEIS Foundation Inc.

License Agreements, Terms of Use, Privacy Policy. .

Last modified October 5 21:37 EDT 2022. Contains 357261 sequences. (Running on oeis4.)