/*  Routines for manipulating B-splines.  These are intended for use with
 *  the New S language described in Becker, Chambers, and Wilks (Wadsworth,
 *  1988).  
 *  Copyright 1989 by Douglas M. Bates.  Permission to use, copy,
 *  modify and distribute is freely granted to all.
 *  The routines are loosely based on the pseudo-code in Schumaker (Wiley,
 *  1981) and the CMLIB library DBSPLINES.
*/

#include <memory.h>
#include "Rdll.h"

#ifdef __STDC__
extern void *S_alloc();
#else
extern char *S_alloc();
#endif

REXPORT void spline_value(double*, double*, long*, long*, double*, long*,
        long*, double*);
REXPORT void spline_basis(double*, long*, long*, double*, long*, long*,
  double*, long*);
REXPORT void lin_interp(double*, double*, double*, double*, long*);


static double *ldel, *rdel;
static long orderm1;

static void diff_table(ti, x, ndiff)
     double *ti, x;
     register long ndiff;
    
{
  register double *r = rdel, *l = ldel, *dpt = ti;
  
  while (ndiff--) {
    *r++ = *dpt++ - x;
    *l++ = x - *--ti;
  }
}
      

static void basis_funcs(ti, x, b)
     double *ti, x, *b;
{
  long j, r;
  double saved, term;

  diff_table(ti, x, orderm1);
  b[0] = 1.;
  for (j = 1; j <= orderm1; j++) {
    saved = 0.;
    for (r = 0; r < j; r++) {
      term = b[r]/(rdel[r] + ldel[j - 1 - r]);
      b[r] = saved + rdel[r] * term;
      saved = ldel[j - 1 - r] * term;
    }
    b[j] = saved;
  }
}  

static double evaluate(ti, x, a, nder)
     double *ti, x, *a;
{
  register double *lpt, *rpt, *apt;
  register int inner;
  int outer = orderm1;

  while(nder--) {
    for(inner = outer, apt = a, lpt = ti - outer; inner--; apt++, lpt++)
      *apt = outer * (*(apt + 1) - *apt)/(*(lpt + outer) - *lpt);
    outer--;
  }
  diff_table(ti, x, (long) outer);
  while(outer--)
    for(apt = a, lpt = ldel + outer, rpt = rdel, inner = outer + 1;
        inner--; lpt--, rpt++, apt++)
      *apt = (*(apt + 1) * *lpt + *apt * *rpt)/(*rpt + *lpt);
  return(*a);
}  
  
void spline_value(knots, coeff, ncoeff, order, x, nx, deriv, y)
     double *knots, *coeff, *x, *y;
     long *ncoeff, *order, *nx, *deriv;
{
  long n = *nx;
  double *a, *last = knots + *ncoeff;

  a = (double *) S_alloc(*order, sizeof(double));
  orderm1 = *order - 1L;        /* allocate difference tables */
  rdel = (double *) S_alloc(orderm1, sizeof(double));
  ldel = (double *) S_alloc(orderm1, sizeof(double));

  knots += *order;              /* First *order knots must be < all x's */
  while(n--) {
    while(knots < last && *knots <= *x) {knots++; coeff++;}
    memcpy((char *) a, (char *) coeff, (int) *order * sizeof(double));
    *y++ = evaluate(knots, *x++, a, (int) *deriv);
  }
}

void spline_basis(knots, ncoeff, order, xvals, derivs, nx, basis, offsets)
     double *knots, *xvals, *basis;
     long *ncoeff, *order, *derivs, *nx, *offsets;
{                               /* evaluate the non-zero B-spline basis */
                                /* functions (or their derivatives) at */
                                /* xvals.  */
  long n = *nx, i, j;
  double *dpt, *coeff, *last = knots + *ncoeff;

  dpt = (knots += *order);      /* first order knots must be <= all xvals */
  orderm1 = *order - 1L;
  rdel = (double *) S_alloc(orderm1, sizeof(double));
  ldel = (double *) S_alloc(orderm1, sizeof(double));
  coeff = (double *) S_alloc(*order, sizeof(double));
  for( ; n--; xvals++, derivs++) {
    while(dpt < last && *dpt <= *xvals) dpt++;
    if (*derivs) {              /* slow method for derivatives */
      for(i = 0; i < *order; i++) {
        for(j = 0; j < *order; j++) coeff[j] = 0;
        coeff[i] = 1;
        *basis++ = evaluate(dpt, *xvals, coeff, (int) *derivs);
      }
    }
    else {
      basis_funcs(dpt, *xvals, basis); /* fast method for value */
      basis += *order;
    }
    *offsets++ = (long)(dpt - knots);
  }
}

void lin_interp(x, y, x0, y0, nvals)
     double *x, *y, *x0, *y0;
     long *nvals;
{
  long n = *nvals;

  while(n--) {
    while (*x < *x0) {x++; y++;}
    if (*x > *x0) *y0++ = *y + (*(y+1) - *y)*(*x0 - *x)/(*(x+1) - *x);
    else *y0++ = *y;
    x0++;
  }
}
