+static char *sccsid = "@(#)spline.c 4.1 (Berkeley) %G%";
+#include <stdio.h>
+
+#define NP 1000
+#define INF 1.e37
+
+struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y;
+float *diag, *r;
+float dx = 1.;
+float ni = 100.;
+int n;
+int auta;
+int periodic;
+float konst = 0.0;
+float zero = 0.;
+
+/* Spline fit technique
+let x,y be vectors of abscissas and ordinates
+ h be vector of differences h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
+ y" be vector of 2nd derivs of approx function
+If the points are numbered 0,1,2,...,n+1 then y" satisfies
+(R W Hamming, Numerical Methods for Engineers and Scientists,
+2nd Ed, p349ff)
+ h\e9i\e8y"\b\e9i-1\e9\e8\e8+2(h\e9i\e8+h\e9i+1\e8)y"\b\e9i\e8+h\e9i+1\e8y"\b\e9i+1\e8
+
+ = 6[(y\e9i+1\e8-y\e9i\e8)/h\e9i+1\e8-(y\e9i\e8-y\e9i-1\e8)/h\e9i\e8] i=1,2,...,n
+
+where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
+This is a symmetric tridiagonal system of the form
+
+ | a\e91\e8 h\e92\e8 | |y"\b\e91\e8| |b\e91\e8|
+ | h\e92\e8 a\e92\e8 h\e93\e8 | |y"\b\e92\e8| |b\e92\e8|
+ | h\e93\e8 a\e93\e8 h\e94\e8 | |y"\b\e93\e8| = |b\e93\e8|
+ | . | | .| | .|
+ | . | | .| | .|
+It can be triangularized into
+ | d\e91\e8 h\e92\e8 | |y"\b\e91\e8| |r\e91\e8|
+ | d\e92\e8 h\e93\e8 | |y"\b\e92\e8| |r\e92\e8|
+ | d\e93\e8 h\e94\e8 | |y"\b\e93\e8| = |r\e93\e8|
+ | . | | .| | .|
+ | . | | .| | .|
+where
+ d\e91\e8 = a\e91\e8
+
+ r\e90\e8 = 0
+
+ d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8 1<i<\b_n
+
+ r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8 1<\b_i<\b_n
+
+the back solution is
+ y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8
+
+ y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8 1<\b_i<n
+
+superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
+recalculated backward by the formulas
+
+ d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8) 1<i<\b_n
+
+ r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8 1<i<\b_n
+
+unhappily it turns out that the recursion forward for d
+is quite strongly geometrically convergent--and is wildly
+unstable going backward.
+There's similar trouble with r, so the intermediate
+results must be kept.
+
+Note that n-1 in the program below plays the role of n+1 in the theory
+
+Other boundary conditions\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b_________________________
+
+The boundary conditions are easily generalized to handle
+
+ y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b" = ky\e9n\e8\b"
+
+for some constant k. The above analysis was for k = 0;
+k = 1 fits parabolas perfectly as well as stright lines;
+k = 1/2 has been recommended as somehow pleasant.
+
+All that is necessary is to add h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.
+
+
+Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________
+
+To do this, add 1 more row and column thus
+
+ | a\e91\e8 h\e92\e8 h\e91\e8 | |y\e91\e8\b"| |b\e91\e8|
+ | h\e92\e8 a\e92\e8 h\e93\e8 | |y\e92\e8\b"| |b\e92\e8|
+ | h\e93\e8 a\e94\e8 h\e94\e8 | |y\e93\e8\b"| |b\e93\e8|
+ | | | .| = | .|
+ | . | | .| | .|
+ | h\e91\e8 h\e90\e8 a\e90\e8 | | .| | .|
+
+where h\e90\e8=\b_ h\e9n+1\e8
+
+The same diagonalization procedure works, except for
+the effect of the 2 corner elements. Let s\e9i\e8 be the part
+of the last element in the i\e8th\e9 "diagonalized" row that
+arises from the extra top corner element.
+
+ s\e91\e8 = h\e91\e8
+
+ s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8 2<\b_i<\b_n+1
+
+After "diagonalizing", the lower corner element remains.
+Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
+as the bottom element to its left is eliminated
+
+ t\e91\e8 = h\e91\e8
+
+ t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8
+
+Evidently t\e9i\e8 = s\e9i\e8.
+Elimination along the bottom row
+introduces further corrections to the bottom right element
+and to the last element of the right hand side.
+Call these corrections u and v.
+
+ u\e91\e8 = v\e91\e8 = 0
+
+ u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8
+
+ v\e9i\e8 = v\e9i-1\e8-r\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8 2<\b_i<\b_n+1
+
+The back solution is now obtained as follows
+
+ y"\b\e9n+1\e8 = (r\e9n+1\e8+v\e9n+1\e8)/(d\e9n+1\e8+s\e9n+1\e8+t\e9n+1\e8+u\e9n+1\e8)
+
+ y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8*y\e9i+1\e8-s\e9i\e8*y\e9n+1\e8)/d\e9i\e8 1<\b_i<\b_n
+
+Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula
+
+ y = y\e9i\e8x\e9+\e8 + y\e9i+1\e8x\e9-\e8 -(h\e82\e9\b\e9i+1\e8/6)[y"\b\e9i\e8(x\e9+\e8-x\e9+\e8\e8\b3\e9)+y"\b\e9i+1\e8(x\e9-\e8-x\e9-\e8\b\e83\e9)]
+where
+ x\e9+\e8 = x\e9i+1\e8-x
+
+ x\e9-\e8 = x-x\e9i\e8
+*/
+
+float
+rhs(i){
+ int i_;
+ double zz;
+ i_ = i==n-1?0:i;
+ zz = (y.val[i]-y.val[i-1])/(x.val[i]-x.val[i-1]);
+ return(6*((y.val[i_+1]-y.val[i_])/(x.val[i+1]-x.val[i]) - zz));
+}
+
+spline(){
+ float d,s,u,v,hi,hi1;
+ float h;
+ float D2yi,D2yi1,D2yn1,x0,x1,yy,a;
+ int end;
+ float corr;
+ int i,j,m;
+ if(n<3) return(0);
+ if(periodic) konst = 0;
+ d = 1;
+ r[0] = 0;
+ s = periodic?-1:0;
+ for(i=0;++i<n-!periodic;){ /* triangularize */
+ hi = x.val[i]-x.val[i-1];
+ hi1 = i==n-1?x.val[1]-x.val[0]:
+ x.val[i+1]-x.val[i];
+ if(hi1*hi<=0) return(0);
+ u = i==1?zero:u-s*s/d;
+ v = i==1?zero:v-s*r[i-1]/d;
+ r[i] = rhs(i)-hi*r[i-1]/d;
+ s = -hi*s/d;
+ a = 2*(hi+hi1);
+ if(i==1) a += konst*hi;
+ if(i==n-2) a += konst*hi1;
+ diag[i] = d = i==1? a:
+ a - hi*hi/d;
+ }
+ D2yi = D2yn1 = 0;
+ for(i=n-!periodic;--i>=0;){ /* back substitute */
+ end = i==n-1;
+ hi1 = end?x.val[1]-x.val[0]:
+ x.val[i+1]-x.val[i];
+ D2yi1 = D2yi;
+ if(i>0){
+ hi = x.val[i]-x.val[i-1];
+ corr = end?2*s+u:zero;
+ D2yi = (end*v+r[i]-hi1*D2yi1-s*D2yn1)/
+ (diag[i]+corr);
+ if(end) D2yn1 = D2yi;
+ if(i>1){
+ a = 2*(hi+hi1);
+ if(i==1) a += konst*hi;
+ if(i==n-2) a += konst*hi1;
+ d = diag[i-1];
+ s = -s*d/hi;
+ }}
+ else D2yi = D2yn1;
+ if(!periodic) {
+ if(i==0) D2yi = konst*D2yi1;
+ if(i==n-2) D2yi1 = konst*D2yi;
+ }
+ if(end) continue;
+ m = hi1>0?ni:-ni;
+ m = 1.001*m*hi1/(x.ub-x.lb);
+ if(m<=0) m = 1;
+ h = hi1/m;
+ for(j=m;j>0||i==0&&j==0;j--){ /* interpolate */
+ x0 = (m-j)*h/hi1;
+ x1 = j*h/hi1;
+ yy = D2yi*(x0-x0*x0*x0)+D2yi1*(x1-x1*x1*x1);
+ yy = y.val[i]*x0+y.val[i+1]*x1 -hi1*hi1*yy/6;
+ printf("%f ",x.val[i]+j*h);
+ printf("%f\n",yy);
+ }
+ }
+ return(1);
+ }
+readin() {
+ for(n=0;n<NP;n++){
+ if(auta) x.val[n] = n*dx+x.lb;
+ else if(!getfloat(&x.val[n])) break;
+ if(!getfloat(&y.val[n])) break; } }
+
+getfloat(p)
+ float *p;{
+ char buf[30];
+ register c;
+ int i;
+ extern double atof();
+ for(;;){
+ c = getchar();
+ if (c==EOF) {
+ *buf = '\0';
+ return(0);
+ }
+ *buf = c;
+ switch(*buf){
+ case ' ':
+ case '\t':
+ case '\n':
+ continue;}
+ break;}
+ for(i=1;i<30;i++){
+ c = getchar();
+ if (c==EOF) {
+ buf[i] = '\0';
+ break;
+ }
+ buf[i] = c;
+ if('0'<=c && c<='9') continue;
+ switch(c) {
+ case '.':
+ case '+':
+ case '-':
+ case 'E':
+ case 'e':
+ continue;}
+ break; }
+ buf[i] = ' ';
+ *p = atof(buf);
+ return(1); }
+
+getlim(p)
+ struct proj *p; {
+ int i;
+ for(i=0;i<n;i++) {
+ if(!p->lbf && p->lb>(p->val[i])) p->lb = p->val[i];
+ if(!p->ubf && p->ub<(p->val[i])) p->ub = p->val[i]; }
+ }
+
+
+main(argc,argv)
+ char *argv[];{
+ extern char *malloc();
+ int i;
+ x.lbf = x.ubf = y.lbf = y.ubf = 0;
+ x.lb = INF;
+ x.ub = -INF;
+ y.lb = INF;
+ y.ub = -INF;
+ while(--argc > 0) {
+ argv++;
+again: switch(argv[0][0]) {
+ case '-':
+ argv[0]++;
+ goto again;
+ case 'a':
+ auta = 1;
+ numb(&dx,&argc,&argv);
+ break;
+ case 'k':
+ numb(&konst,&argc,&argv);
+ break;
+ case 'n':
+ numb(&ni,&argc,&argv);
+ break;
+ case 'p':
+ periodic = 1;
+ break;
+ case 'x':
+ if(!numb(&x.lb,&argc,&argv)) break;
+ x.lbf = 1;
+ if(!numb(&x.ub,&argc,&argv)) break;
+ x.ubf = 1;
+ break;
+ default:
+ fprintf(stderr, "Bad agrument\n");
+ exit(1);
+ }
+ }
+ if(auta&&!x.lbf) x.lb = 0;
+ readin();
+ getlim(&x);
+ getlim(&y);
+ i = (n+1)*sizeof(dx);
+ diag = (float *)malloc((unsigned)i);
+ r = (float *)malloc((unsigned)i);
+ if(r==NULL||!spline()) for(i=0;i<n;i++){
+ printf("%f ",x.val[i]);
+ printf("%f\n",y.val[i]); }
+}
+numb(np,argcp,argvp)
+ int *argcp;
+ float *np;
+ char ***argvp;{
+ double atof();
+ char c;
+ if(*argcp<=1) return(0);
+ c = (*argvp)[1][0];
+ if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0);
+ *np = atof((*argvp)[1]);
+ (*argcp)--;
+ (*argvp)++;
+ return(1); }
+
projects / unix-history / commitdiff