bridge/class-master/depend/pari/src/language/sumiter.c
/* Fornberg interpolation algorithm for finite differences coefficients
* using N+1 equidistant grid points around 0 [ assume N even >= M ].
* Compute \delta[m]_{N,nu} for all derivation orders m = 0..M such that
*   h^m * f^{(m)}(0) = \sum_{nu = 0}^n delta[m]_{N,nu}  f(a_nu) + O(h^{N-m+1}),
* for step size h.
* Return a = [0,-1,1...,-N2,N2] and vector of vectors d: d[m+1][nu+1]
* = (M!/m!) * delta[m]_{N,nu}, nu = 0..N */
Fornberg interpolation algorithm for finite differences coefficients
using N+1 equidistant grid points around 0 [ assume N even >= M ].
Fornberg interpolation algorithm for finite differences coefficients using N + 1 equidistant grid points around 0 [ assume N even > = M ] .
static void
FD(long M, long N, GEN *pd, GEN *pa)
{
  GEN d, a, b, W, Wp, t, F, Mfact;
  long N2, m, nu, i;
  if (odd(N)) N++; 
  N2 = N>>1;
  F = cgetg(N+2, t_VEC);
  a = cgetg(N+2, t_VEC);
  b = cgetg(N2+1, t_VEC);
  gel(a,1) = gen_0;
  for (i = 1; i <= N2; i++)
  {
    gel(a,2*i)   = utoineg(i);
    gel(a,2*i+1) = utoipos(i);
    gel(b,i) = sqru(i);
  }
  Mfact = mpfact(M);
  W = roots_to_pol(b, 0);
  Wp = ZX_deriv(W);
  t = gel(W,2); 
  t = diviiexact(t, Mfact);
  gel(F,1) = RgX_Rg_div(RgX_inflate(W,2), t);
  for (i = 1; i <= N2; i++)
  { 
    GEN r, t = mulii(shifti(gel(b,i),1), poleval(Wp, gel(b,i)));
    GEN U, S, T;
    U = RgX_inflate(RgX_div_by_X_x(W, gel(b,i), &r), 2);
    U = RgX_shift_shallow(U, 1);
    U = RgXn_red_shallow(U, M+1); 
    t = diviiexact(t, Mfact);
    U = RgX_Rg_div(U, t);
    S = RgX_shift_shallow(U,1);
    T = RgX_Rg_mul(U, gel(a,2*i+1));
    gel(F,2*i)   = RgX_sub(S, T);
    gel(F,2*i+1) = RgX_add(S, T);
  }
  d = cgetg(M+2, t_VEC);
  for (m = 0; m <= M; m++)
  {
    GEN v = cgetg(N+2, t_VEC); 
    for (nu = 0; nu <= N; nu++) gel(v, nu+1) = gmael(F, nu+1, m+2);
    gel(d,m+1) = v;
  }
  *pd = d;
  *pa = a;
}
static void FD ( long M , long N , GEN * pd , GEN * pa ) { GEN d , a , b , W , Wp , t , F , Mfact ; long N2 , m , nu , i ; if ( odd ( N ) ) N ++ ; N2 = N >> 1 ; F = cgetg ( N + 2 , t_VEC ) ; a = cgetg ( N + 2 , t_VEC ) ; b = cgetg ( N2 + 1 , t_VEC ) ; gel ( a , 1 ) = gen_0 ; for ( i = 1 ; i <= N2 ; i ++ ) { gel ( a , 2 * i ) = utoineg ( i ) ; gel ( a , 2 * i + 1 ) = utoipos ( i ) ; gel ( b , i ) = sqru ( i ) ; } Mfact = mpfact ( M ) ; W = roots_to_pol ( b , 0 ) ; Wp = ZX_deriv ( W ) ; t = gel ( W , 2 ) ; t = diviiexact ( t , Mfact ) ; gel ( F , 1 ) = RgX_Rg_div ( RgX_inflate ( W , 2 ) , t ) ; for ( i = 1 ; i <= N2 ; i ++ ) { GEN r , t = mulii ( shifti ( gel ( b , i ) , 1 ) , poleval ( Wp , gel ( b , i ) ) ) ; GEN U , S , T ; U = RgX_inflate ( RgX_div_by_X_x ( W , gel ( b , i ) , & r ) , 2 ) ; U = RgX_shift_shallow ( U , 1 ) ; U = RgXn_red_shallow ( U , M + 1 ) ; t = diviiexact ( t , Mfact ) ; U = RgX_Rg_div ( U , t ) ; S = RgX_shift_shallow ( U , 1 ) ; T = RgX_Rg_mul ( U , gel ( a , 2 * i + 1 ) ) ; gel ( F , 2 * i ) = RgX_sub ( S , T ) ; gel ( F , 2 * i + 1 ) = RgX_add ( S , T ) ; } d = cgetg ( M + 2 , t_VEC ) ; for ( m = 0 ; m <= M ; m ++ ) { GEN v = cgetg ( N + 2 , t_VEC ) ; for ( nu = 0 ; nu <= N ; nu ++ ) gel ( v , nu + 1 ) = gmael ( F , nu + 1 , m + 2 ) ; gel ( d , m + 1 ) = v ; } * pd = d ; * pa = a ; }
Fornberg interpolation algorithm for finite differences coefficients
using N+1 equidistant grid points around 0 [ assume N even >= M ].
Fornberg interpolation algorithm for finite differences coefficients using N + 1 equidistant grid points around 0 [ assume N even > = M ] .
/* make it even */ /* w = \prod (X - a[i]) = x W(x^2) */ /* w'(0) */ /* t = w'(a_{2i}) = w'(a_{2i+1}) */ /* higher terms not needed */ /* F[i] = M! w(X) / ((X-a[i])w'(a[i])) + O(X^(M+1)) */ /* coeff(F[nu],X^m) */