[unix-history] / usr / src / cmd / spline.c

#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); }
Commit	Line	Data
9c79dca0 BJ	1	#include <stdio.h>
	2
	3	#define NP 1000
	4	#define INF 1.e37
	5
	6	struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y;
	7	float diag, r;
	8	float dx = 1.;
	9	float ni = 100.;
	10	int n;
	11	int auta;
	12	int periodic;
	13	float konst = 0.0;
	14	float zero = 0.;
	15
	16	/* Spline fit technique
	17	let x,y be vectors of abscissas and ordinates
	18	h be vector of differences h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
	19	y" be vector of 2nd derivs of approx function
	20	If the points are numbered 0,1,2,...,n+1 then y" satisfies
	21	(R W Hamming, Numerical Methods for Engineers and Scientists,
	22	2nd Ed, p349ff)
	23	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
	24
	25	= 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
	26
	27	where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
	28	This is a symmetric tridiagonal system of the form
	29
	30	\| a\e91\e8 h\e92\e8 \| \|y"\b\e91\e8\| \|b\e91\e8\|
	31	\| h\e92\e8 a\e92\e8 h\e93\e8 \| \|y"\b\e92\e8\| \|b\e92\e8\|
	32	\| h\e93\e8 a\e93\e8 h\e94\e8 \| \|y"\b\e93\e8\| = \|b\e93\e8\|
	33	\| . \| \| .\| \| .\|
	34	\| . \| \| .\| \| .\|
	35	It can be triangularized into
	36	\| d\e91\e8 h\e92\e8 \| \|y"\b\e91\e8\| \|r\e91\e8\|
	37	\| d\e92\e8 h\e93\e8 \| \|y"\b\e92\e8\| \|r\e92\e8\|
	38	\| d\e93\e8 h\e94\e8 \| \|y"\b\e93\e8\| = \|r\e93\e8\|
	39	\| . \| \| .\| \| .\|
	40	\| . \| \| .\| \| .\|
	41	where
	42	d\e91\e8 = a\e91\e8
	43
	44	r\e90\e8 = 0
	45
	46	d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8 1<i<\b_n
	47
	48	r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8 1<\b_i<\b_n
	49
	50	the back solution is
	51	y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8
	52
	53	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8 1<\b_i<n
	54
	55	superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
	56	recalculated backward by the formulas
	57
	58	d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8) 1<i<\b_n
	59
	60	r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8 1<i<\b_n
	61
	62	unhappily it turns out that the recursion forward for d
	63	is quite strongly geometrically convergent--and is wildly
	64	unstable going backward.
65	There's similar trouble with r, so the intermediate
66	results must be kept.
67
68	Note that n-1 in the program below plays the role of n+1 in the theory
69
70	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_________________________
71
72	The boundary conditions are easily generalized to handle
73
74	y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b" = ky\e9n\e8\b"
75
76	for some constant k. The above analysis was for k = 0;
77	k = 1 fits parabolas perfectly as well as stright lines;
78	k = 1/2 has been recommended as somehow pleasant.
79
80	All that is necessary is to add h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.
81
82
83	Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________
84
85	To do this, add 1 more row and column thus
86
87	\| a\e91\e8 h\e92\e8 h\e91\e8 \| \|y\e91\e8\b"\| \|b\e91\e8\|
88	\| h\e92\e8 a\e92\e8 h\e93\e8 \| \|y\e92\e8\b"\| \|b\e92\e8\|
89	\| h\e93\e8 a\e94\e8 h\e94\e8 \| \|y\e93\e8\b"\| \|b\e93\e8\|
90	\| \| \| .\| = \| .\|
91	\| . \| \| .\| \| .\|
92	\| h\e91\e8 h\e90\e8 a\e90\e8 \| \| .\| \| .\|
93
94	where h\e90\e8=\b_ h\e9n+1\e8
95
96	The same diagonalization procedure works, except for
97	the effect of the 2 corner elements. Let s\e9i\e8 be the part
98	of the last element in the i\e8th\e9 "diagonalized" row that
99	arises from the extra top corner element.
100
101	s\e91\e8 = h\e91\e8
102
103	s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8 2<\b_i<\b_n+1
104
105	After "diagonalizing", the lower corner element remains.
106	Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
107	as the bottom element to its left is eliminated
108
109	t\e91\e8 = h\e91\e8
110
111	t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8
112
113	Evidently t\e9i\e8 = s\e9i\e8.
114	Elimination along the bottom row
115	introduces further corrections to the bottom right element
116	and to the last element of the right hand side.
117	Call these corrections u and v.
118
119	u\e91\e8 = v\e91\e8 = 0
120
121	u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8
122
123	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
124
125	The back solution is now obtained as follows
126
127	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)
128
129	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y\e9i+1\e8-s\e9i\e8y\e9n+1\e8)/d\e9i\e8 1<\b_i<\b_n
130
131	Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula
132
133	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)]
134	where
135	x\e9+\e8 = x\e9i+1\e8-x
136
137	x\e9-\e8 = x-x\e9i\e8
138	*/
139
140	float
141	rhs(i){
142	int i_;
143	double zz;
144	i_ = i==n-1?0:i;
145	zz = (y.val[i]-y.val[i-1])/(x.val[i]-x.val[i-1]);
146	return(6*((y.val[i_+1]-y.val[i_])/(x.val[i+1]-x.val[i]) - zz));
147	}
148
149	spline(){
150	float d,s,u,v,hi,hi1;
151	float h;
152	float D2yi,D2yi1,D2yn1,x0,x1,yy,a;
153	int end;
154	float corr;
155	int i,j,m;
156	if(n<3) return(0);
157	if(periodic) konst = 0;
158	d = 1;
159	r[0] = 0;
160	s = periodic?-1:0;
161	for(i=0;++i<n-!periodic;){ /* triangularize */
162	hi = x.val[i]-x.val[i-1];
163	hi1 = i==n-1?x.val[1]-x.val[0]:
164	x.val[i+1]-x.val[i];
165	if(hi1*hi<=0) return(0);
166	u = i==1?zero:u-s*s/d;
167	v = i==1?zero:v-s*r[i-1]/d;
168	r[i] = rhs(i)-hi*r[i-1]/d;
169	s = -hi*s/d;
170	a = 2*(hi+hi1);
171	if(i==1) a += konst*hi;
172	if(i==n-2) a += konst*hi1;
173	diag[i] = d = i==1? a:
174	a - hi*hi/d;
175	}
176	D2yi = D2yn1 = 0;
177	for(i=n-!periodic;--i>=0;){ /* back substitute */
178	end = i==n-1;
179	hi1 = end?x.val[1]-x.val[0]:
180	x.val[i+1]-x.val[i];
181	D2yi1 = D2yi;
182	if(i>0){
183	hi = x.val[i]-x.val[i-1];
184	corr = end?2*s+u:zero;
185	D2yi = (endv+r[i]-hi1D2yi1-s*D2yn1)/
186	(diag[i]+corr);
187	if(end) D2yn1 = D2yi;
188	if(i>1){
189	a = 2*(hi+hi1);
190	if(i==1) a += konst*hi;
191	if(i==n-2) a += konst*hi1;
192	d = diag[i-1];
193	s = -s*d/hi;
194	}}
195	else D2yi = D2yn1;
196	if(!periodic) {
197	if(i==0) D2yi = konst*D2yi1;
198	if(i==n-2) D2yi1 = konst*D2yi;
199	}
200	if(end) continue;
201	m = hi1>0?ni:-ni;
202	m = 1.001mhi1/(x.ub-x.lb);
203	if(m<=0) m = 1;
204	h = hi1/m;
205	for(j=m;j>0\|\|i==0&&j==0;j--){ /* interpolate */
206	x0 = (m-j)*h/hi1;
207	x1 = j*h/hi1;
208	yy = D2yi(x0-x0x0x0)+D2yi1(x1-x1x1x1);
209	yy = y.val[i]x0+y.val[i+1]x1 -hi1hi1yy/6;
210	printf("%f ",x.val[i]+j*h);
211	printf("%f\n",yy);
212	}
213	}
214	return(1);
215	}
216	readin() {
217	for(n=0;n<NP;n++){
218	if(auta) x.val[n] = n*dx+x.lb;
219	else if(!getfloat(&x.val[n])) break;
220	if(!getfloat(&y.val[n])) break; } }
221
222	getfloat(p)
223	float *p;{
224	char buf[30];
225	register c;
226	int i;
227	extern double atof();
228	for(;;){
229	c = getchar();
230	if (c==EOF) {
231	*buf = '\0';
232	return(0);
233	}
234	*buf = c;
235	switch(*buf){
236	case ' ':
237	case '\t':
238	case '\n':
239	continue;}
240	break;}
241	for(i=1;i<30;i++){
242	c = getchar();
243	if (c==EOF) {
244	buf[i] = '\0';
245	break;
246	}
247	buf[i] = c;
248	if('0'<=c && c<='9') continue;
249	switch(c) {
250	case '.':
251	case '+':
252	case '-':
253	case 'E':
254	case 'e':
255	continue;}
256	break; }
257	buf[i] = ' ';
258	*p = atof(buf);
259	return(1); }
260
261	getlim(p)
262	struct proj *p; {
263	int i;
264	for(i=0;i<n;i++) {
265	if(!p->lbf && p->lb>(p->val[i])) p->lb = p->val[i];
266	if(!p->ubf && p->ub<(p->val[i])) p->ub = p->val[i]; }
267	}
268
269
270	main(argc,argv)
271	char *argv[];{
272	extern char *malloc();
273	int i;
274	x.lbf = x.ubf = y.lbf = y.ubf = 0;
275	x.lb = INF;
276	x.ub = -INF;
277	y.lb = INF;
278	y.ub = -INF;
279	while(--argc > 0) {
280	argv++;
281	again: switch(argv[0][0]) {
282	case '-':
283	argv[0]++;
284	goto again;
285	case 'a':
286	auta = 1;
287	numb(&dx,&argc,&argv);
288	break;
289	case 'k':
290	numb(&konst,&argc,&argv);
291	break;
292	case 'n':
293	numb(&ni,&argc,&argv);
294	break;
295	case 'p':
296	periodic = 1;
297	break;
298	case 'x':
299	if(!numb(&x.lb,&argc,&argv)) break;
300	x.lbf = 1;
301	if(!numb(&x.ub,&argc,&argv)) break;
302	x.ubf = 1;
303	break;
304	default:
305	fprintf(stderr, "Bad agrument\n");
306	exit(1);
307	}
308	}
309	if(auta&&!x.lbf) x.lb = 0;
310	readin();
311	getlim(&x);
312	getlim(&y);
313	i = (n+1)*sizeof(dx);
314	diag = (float *)malloc((unsigned)i);
315	r = (float *)malloc((unsigned)i);
316	if(r==NULL\|\|!spline()) for(i=0;i<n;i++){
317	printf("%f ",x.val[i]);
318	printf("%f\n",y.val[i]); }
319	}
320	numb(np,argcp,argvp)
321	int *argcp;
322	float *np;
323	char ***argvp;{
324	double atof();
325	char c;
326	if(*argcp<=1) return(0);
327	c = (*argvp)[1][0];
328	if(!('0'<=c&&c<='9' \|\| c=='-' \|\| c== '.' )) return(0);
329	np = atof((argvp)[1]);
330	(*argcp)--;
331	(*argvp)++;
332	return(1); }
333