#include <stdio.h>

inline double
_bspline_ ( int i, int m, double t, double *u )
{
  register double b1, b2;
  if ( m == 1 )
    return 1.0/( u[i]-u[i-1] );

  b1 = ( t >= u[i-m] && t < u[i-1] ) ? ( t-u[i-m] ) *_bspline_ ( i-1,m-1, t, u) : 0.0;
  b2 = ( t >= u[i-m+1] && t < u[i] ) ? ( u[i]-t ) *_bspline_ ( i,m-1,t, u ) : 0.0;
  return ( b1+b2 ) / ( u[i]-u[i-m] );
}

inline double
_derivative_ ( int k, int i, int m, double t, double *u )
{
  register double d1, d2;
  if ( k >= m )
    return 0.0;
  if ( t < u[i-m] || t >= u[i] )
    return 0.0;
  if ( k == 0 )
    return _bspline_ ( i, m, t, u);

  d1 = ( t >= u[i-m] && t < u[i-1] ) ? k*_derivative_ ( k-1,i-1,m-1,t, u ) +
                                       ( t-u[i-m] ) *_derivative_ ( k,i-1,m-1,t, u ) : 0.0;

  d2 = ( t >= u[i-m+1] && t < u[i] ) ? -k*_derivative_ ( k-1,i,m-1,t, u ) +
                                       ( u[i]-t ) *_derivative_ ( k,i,m-1,t, u ) : 0.0;

  return ( d1+d2 )/( u[i]-u[i-m] );
}

int main(int argc, char **argv)
{
  int number_of_nodes = 11;
  long int grid_size = 10000000;
  double u, du;
  double u_min, u_max;

  double nodes[11];

  u_min = 0.0;
  u_max = 1.0;

  du = (u_max - u_min)/(number_of_nodes-1);

  for (int i = 0; i < number_of_nodes; ++i) {
    nodes[i] = u_min + du*i;
  }

  du = (u_max - u_min)/(grid_size-1);
  printf("%e\n", du);
  for (long int i = 0; i < grid_size; ++i) {
    u = u_min + du*i;
    _bspline_(4, 4, u, nodes);
  }

  return 0;
}