.

2 files changed, 821 insertions, 0 deletions

diff --git a/src/factor.c b/src/factor.c
new file mode 100644
index 000000000..967799e1c
--- /dev/null
+++ b/src/factor.c
@@ -0,0 +1,86 @@
+/* factor -- print factors of n.  lose if n > 2^32.
+   Copyright (C) 1986 Free Software Foundation, Inc.
+
+   This program is free software; you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 2, or (at your option)
+   any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program; if not, write to the Free Software
+   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+/* Written by Paul Rubin <phr@ocf.berkeley.edu>.  */
+
+#include <stdio.h>
+
+void do_stdin ();
+void factor ();
+
+void
+main (argc, argv)
+     int argc;
+     char **argv;
+{
+  if (argc == 1)
+    do_stdin ();
+  else if (argc == 2)
+    factor ((unsigned) atoi (argv[1]));
+  else
+    {
+      fprintf (stderr, "Usage: %s [number]\n", argv[0]);
+      exit (1);
+    }
+  exit (0);
+}
+
+void
+factor (n0)
+     unsigned long n0;
+{
+  register unsigned long n = n0, d;
+
+  if (n < 1)
+    return;
+
+  while (n % 2 == 0)
+    {
+      printf ("\t2\n");
+      n /= 2;
+    }
+  /* The exit condition in the following loop is correct because
+     any time it is tested one of these 3 conditions holds:
+      (1) d divides n
+      (2) n is prime
+      (3) n is composite but has no factors less than d.
+     If (1) or (2) obviously the right thing happens.
+     If (3), then since n is composite it is >= d^2. */
+  for (d = 3; d * d <= n; d += 2)
+    {
+      while (n % d == 0)
+	{
+	  printf ("\t%d\n", d);
+	  n /= d;
+	}
+    }
+  if (n != 1 || n0 == 1)
+    printf ("\t%d\n", n);
+}
+
+void
+do_stdin ()
+{
+  char buf[1000];
+
+  for (;;)
+    {
+      if (fgets (buf, sizeof buf, stdin) == 0)
+	exit (0);
+      factor ((unsigned long) atoi (buf));
+    }
+}
diff --git a/src/spline.c b/src/spline.c
new file mode 100644
index 000000000..fb6a8627a
--- /dev/null
+++ b/src/spline.c
@@ -0,0 +1,735 @@
+/* spline.c -- Do spline interpolation. */
+
+/******************************************************************************
+	David L. Fox
+	2140 Braun Dr.
+	Golden, CO 80401
+
+This program has been placed in the public domain by its author.
+
+Version Date		Change
+1.1	17 Dec 1985	Modify getdata() to realloc() more memory if needed.
+1.0	14 May 1985
+
+spline [options] [file]
+
+Spline reads pairs of numbers from file (or the standard input, if file is missing).
+The pairs are interpreted as abscissas and ordinates of a function.  The output of
+spline consists of similar pairs generated from the input by interpolation with
+cubic splines. (See R. W. Hamming, Numerical Methods for Scientists and Engineers,
+2nd ed., pp. 349ff.)  There are sufficiently many points in the output to appear
+smooth when plotted (e.g., by graph(1)).  The output points are approximately evenly
+spaced and include the input points.
+
+There may be one or more options of the form: -o [argument [arg2] ].
+
+The available options are:
+
+-a	Generate abscissa values automatically.  The input consists of a list of
+	ordinates.  The generated abscissas are evenly spaced with spacing given by
+	the argument, or 1 if the next argument is not a number.
+
+-f	Set the first derivative of the spline function at the left and right end 
+	points to the first and second arguments following -f.  If only one numerical
+	argument follows -f then that value is used for both left and right end points.
+
+-k	The argument of the k option is the value of the constant used to calculate
+	the boundary values by y''[0] = ky''[1] and y''[n] = ky''[n-1].  The default
+	value is k = 0.
+
+-m	The argument gives the maximum number of input data points.  The default is 1000.
+	If the input contains more than this many points additional memory will be
+	allocated.  This option may be used to slightly increase efficiency for large
+	data sets or reduce the amount of memory required for small data sets.
+
+-n	The number of output points is given by the argument.  The default value is 100.
+
+-p	The splines used for interpolation are forced to be periodic, i.e. y'[0] = y'[n].
+	The first and last ordinates should be equal.
+
+-s	Set the second derivative of the spline function at the left and right end 
+	points to the first and second arguments following -s.  If only one numerical
+	argument follows -s then that value is used for both left and right end points.
+
+-x	The argument (and arg2, if present) are the lower (and upper) limits on the
+	abscissa values in the output.  If the x option is not specified these limits
+	are calculated from the data.  Automatic abscissas start at the lower limit
+	(default 0).
+
+The data need not be monotonic in x but all the x values must be distinct.
+Non-numeric data in the input is ignored.
+******************************************************************************/
+
+#include	<a:stdio.h>
+#include	<ctype.h>
+
+/* The constant DOUBLE is a compile time switch.
+	If #define DOUBLE appears here double pecision variables are
+	used to store all data and parameters.
+	Otherwise, float variables are used for the data.
+	For most smoothing and and interpolation applications single
+	precision is more than adequate.
+	Double precision is used to solve the system of linear equations
+	in either case.
+*/
+/* #define DOUBLE	*/
+
+#ifdef	DOUBLE
+#define	double	real; 		/* Type used for data storage. */
+#define	IFMT	"%F"		/* Input format. */
+#define	OFMT	"%18.14g %18.14g\n"	/* Output format. */
+#else
+#define	float	real;		/* Type used for data storage. */
+#define	IFMT	"%f"		/* Input format. */
+#define	OFMT	"%8.5g %8.5g\n"	/* Output format. */
+#endif
+
+/* Numerical constants: These may be machine and/or precision dependent. */
+#define	HUGE	(1.e38)		/* Largest floating point number. */
+#define	EPS	1.e-5		/* Test for zero when solving linear equations. */
+
+/* Default parameters */
+#define	NPTS	1000
+#define	NINT	100
+#define	DEFSTEP	1.	/* Default step size for automatic abcissas. */
+#define	DEFBVK	0.	/* Default boundary value constant. */
+
+/* Boundary condition types. */
+#define	EXTRAP	0	/* Extrapolate second derivative:
+			   y''(0) = k * y''(1), y''(n) = k * y''(n-1) */
+#define	FDERIV	1	/* Fixed first derivatives y'(0) and y'(n). */
+#define	SDERIV	2	/* Fixed second derivatives y''(0) and y''(n). */
+#define	PERIOD	3	/* Periodic:  derivatives equal at end points. */
+
+/* Token types for command line processing. */
+#define	OPTION	1
+#define	NUMBER	2
+#define	OTHER	3
+
+/* Define error numbers. */
+#define	MEMERR	1
+#define	NODATA	2
+#define	BADOPT	4
+#define	BADFILE	5
+#define	XTRAARG	6
+#define	XTRAPTS	7
+#define	SINGERR	8
+#define	DUPX	9
+#define	BADBC	10
+#define	RANGERR	11
+
+/* Constants and flags are global. */
+int aflag = FALSE;	/* Automatic abscissa flag. */
+real step = DEFSTEP;	/* Spacing for automatic abscissas. */
+int bdcode = EXTRAP;	/* Type of boundary conditions:
+			0 extrapolate f'' with coeficient k.
+			1 first derivatives given
+			2 second derivatives given
+			3 periodic */
+real leftd = 0.,	/* Boundary values of derivatives. */
+   rightd = 0.;
+real k = DEFBVK;	/* Boundary value constant. */
+int rflag = 0;		/* 0: take range from data, 1: min. given, 2: min. & max. given. */
+real xmin = HUGE, xmax;	/* Output range. */
+unsigned nint = NINT;	/* Number of intervals in output. */
+unsigned maxpts = NPTS;	/* Maximun number of points. */
+unsigned nknots = 0;	/* Number of input points. */
+int sflag = FALSE;	/* If TRUE data must be sorted. */
+char *datafile;		/* Name of data file */
+
+int xcompare();		/* Function to compare x values of two data points. */
+
+struct pair {
+	real x, y;
+	} *data;	/* Points to list of data points. */
+
+struct bandm {
+	double diag, upper, rhs;
+	} *eqns;	/* Points to elements of band matrix used to calculate
+			coefficients of the splines. */
+
+main(argc, argv)
+int argc;
+char **argv;
+{
+	setup(argc, argv);
+	if (nknots == 1) {	/* Cannot interpolate with one data point.  */
+		printf(OFMT,data->x, data->y);
+		exit(0);
+	}
+	if (sflag)	/* Sort data if needed. */
+		qsort(data, nknots, sizeof(struct pair), xcompare);
+	calcspline();
+	interpolate();
+}
+
+/* xcompare -- Compare abcissas of two data points (for qsort). */
+xcompare(arg1, arg2)
+struct pair *arg1, *arg2;
+{
+	if (arg1->x > arg2->x)
+		return(1);
+	else if (arg1->x < arg2->x)
+		return(-1);
+	else
+		return(0);
+}
+
+/* error -- Print error message and abort fatal errors. */
+error(errno, auxmsg)
+int errno;
+char *auxmsg;
+{	char *msg;
+	int fatal, usemsg;
+	static char *usage = 
+	"usage: spline [options] [file]\noptions:\n",
+	*options = "-a spacing\n-k const\n-n intervals\n-m points\n-p\n-x xmin xmax\n";
+
+	fatal = TRUE;	/* Default is fatal error. */
+	usemsg = FALSE;	/* Default no usage message. */
+	fprintf(stderr, "spline: ");
+	switch(errno) {
+	case MEMERR:
+		msg = "not enough memory for %u data points\n";
+		break;
+	case NODATA:
+		msg = "data file %s empty\n";
+		break;
+	case BADOPT:
+		msg = "unknown option: %c\n";
+		usemsg = TRUE;
+		break;
+	case BADFILE:
+		msg = "cannot open file: %s\n";
+		break;
+	case XTRAARG:
+		msg = "extra argument ignored: %s\n";
+		fatal = FALSE;
+		usemsg = TRUE;
+		break;
+	case XTRAPTS:
+		fatal = FALSE;
+		msg = "%s";
+		break;
+	case DUPX:
+		msg = "duplicate abcissa value: %s\n";
+		break;
+	case RANGERR:
+		msg = "xmax < xmin not allowed %s\n";
+		break;
+	/* The following errors "can't happen." */
+	/* If they occur some sort of bug is indicated. */
+	case SINGERR:
+		msg = "singular matrix encountered %s\n";
+		break;
+	case BADBC:
+		msg = "internal error: bad boundary value code %d\n";
+		break;
+	default:
+		fprintf(stderr, "unknown error number: %d\n", errno);
+		exit(1);
+	}
+	fprintf(stderr, msg, auxmsg);
+	if (usemsg)
+		fprintf(stderr,"%s%s", usage, options);
+	if (fatal)
+		exit(1);
+}
+
+/* setup -- Initalize all constants and read data. */
+setup(argc, argv)
+int argc;
+char **argv;
+{	char *malloc();
+	FILE fdinp,		/* Source of input. */
+		doarg();
+
+	fdinp = doarg(argc, argv);	/* Process command line arguments. */
+
+	/* Allocate memory for data and band matrix of coefficients. */
+	if ((data = malloc(maxpts*sizeof(struct pair))) == NULL) {
+		error(MEMERR, (char *)maxpts);
+	}
+
+	getdata(fdinp);		/* Read data from fdinp. */
+	if (fdinp != stdin)
+		fclose(fdinp);	/* Close input data file. */
+	if (nknots == 0) {
+		error(NODATA, datafile);
+	}
+	/* Allocate memory for calculation of spline coefficients. */
+	if ((eqns = malloc((nknots+1)*sizeof(struct bandm))) == NULL) {
+		error(MEMERR, (char *)nknots);
+	}
+}
+
+/* doarg -- Process arguments. */
+FILE
+doarg(argc, argv)
+int argc;
+char **argv;
+{	int i, type;
+	double atof();
+	char *s, str[15];
+	FILE fdinp;
+
+	s = argv[i=1];
+	type = gettok(&i, &s, argv, str);
+	do {
+		if (type == OPTION) {
+			switch(*str) {
+			case 'a':	/* Automatic abscissa. */
+				aflag = TRUE;
+				rflag = rflag < 2 ? 1 : rflag;
+				if (xmin == HUGE)	/* Initialize xmin, if needed. */
+					xmin = 0.;
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					if ((step = atof(str)) <= 0.)
+						error(RANGERR,"");
+					type = gettok(&i, &s, argv, str);
+				}
+				break;
+			case 'f':	/* Fix first derivative at boundaries. */
+				bdcode = FDERIV;
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					leftd = atof(str);
+					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+						rightd = atof(str);
+						type = gettok(&i, &s, argv, str);
+					}
+					else {
+						rightd = leftd;
+					}
+				}
+				break;
+			case 'k':	/* Set boundary value constant. */
+				bdcode = EXTRAP;
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					k = atof(str);
+					type = gettok(&i, &s, argv, str);
+				}
+				break;
+			case 'm':	/* Set number of intervals in output. */
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					maxpts = (unsigned)atoi(str);
+					type = gettok(&i, &s, argv, str);
+				}
+				break;
+			case 'n':	/* Set number of intervals in output. */
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					nint = (unsigned)atoi(str);
+					type = gettok(&i, &s, argv, str);
+				}
+				break;
+			case 'p':	/* Require periodic interpolation function. */
+				bdcode = PERIOD;
+				type = gettok(&i, &s, argv, str);
+				break;
+			case 's':	/* Fix second derivative at boundaries. */
+				bdcode = SDERIV;
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					leftd = atof(str);
+					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+						rightd = atof(str);
+						type = gettok(&i, &s, argv, str);
+					}
+					else {
+						rightd = leftd;
+					}
+				}
+				break;
+			case 'x':	/* Set range of x values. */
+				rflag = 1;
+				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+					xmin = atof(str);
+					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
+						xmax = atof(str);
+						rflag = 2;
+						type = gettok(&i, &s, argv, str);
+						if (xmax <= xmin)
+							error(RANGERR, "");
+					}
+				}
+				break;
+			default:
+				error(BADOPT, (char *)argv[i][1]);
+			}
+		}
+		else {
+			if (argc > i) {
+				datafile = argv[i];
+				if ((fdinp = fopen(argv[i++], "r")) == NULL) {
+					error(BADFILE, argv[i-1]);
+				}
+				if (argc > i)
+					error(XTRAARG, argv[i]);
+			}
+			else
+				fdinp = stdin;
+			break;
+		}
+	} while (i < argc);
+	return fdinp;
+}
+
+/* gettok -- Get one token from command line, return type. */
+gettok(indexp, locp, argv, str)
+int *indexp;	/* Pointer to index in argv array. */
+char **locp;	/* Pointer to current location in argv[*indexp]. */
+char **argv;
+char *str;
+{	char *s;
+	char *strcpy(), *strchr();
+	int type;
+
+	s = *locp;
+	while (isspace(*s) || *s == ',')
+		++s;
+	if (*s == '\0')		/* Look at next element in argv. */
+		s = argv[++*indexp];
+	if (*s == '-' && isalpha(s[1])) {
+		/* Found an option. */
+		*str = *++s;
+		str[1] = '\0';
+		++s;	
+		type = OPTION;
+	}
+	else if (isnumber(s)) {
+		while (isnumber(s))
+			*str++ = *s++;
+		*str = '\0';
+		type = NUMBER;
+	}
+	else {
+		strcpy(str, s);
+		s = strchr(s, '\0');
+		type = OTHER;
+	}
+	*locp = s;
+	return(type);
+}
+
+/* isnumber -- Return TRUE if argument is the ASCII representation of a number. */
+isnumber(string)
+char *string;
+{
+	if (isdigit(*string) || 
+	   *string == '.'    ||
+	   *string == '+'   ||
+	   (*string == '-' && (isdigit(string[1]) || string[1] == '.')))
+		return(TRUE);
+	else
+		return(FALSE);
+}
+
+/* getdata -- Read data points from fdinp. */
+getdata(fdinp)
+FILE fdinp;
+{	int n, i;
+	real lastx, min, max;
+	struct pair *dp;
+	char msg[60], *realloc();
+
+	nknots = 0;
+	lastx = -HUGE;
+	dp = data;	/* Point to head of list of data. */
+	min = HUGE;
+	max = -HUGE;
+	do {
+		if (aflag) {	/* Generate abcissa.  */
+			dp->x = xmin + nknots*step;
+			n = 0;
+		}
+		else {		/* Read abcissa. */
+			while ((n = (fscanf(fdinp,IFMT,&dp->x))) == 0)
+				;	/* Skip non-numeric input. */
+		}
+		if (n == 1) {
+			if (min > dp->x)
+				min = dp->x;
+			if (max < dp->x)
+				max = dp->x;
+			if (lastx > dp->x) {		/* Check for monotonicity. */
+				sflag = TRUE;	/* Data must be sorted. */
+			}
+			lastx = dp->x;
+		}
+		/* Read ordinate. */
+		while ((n = (fscanf(fdinp, IFMT, &dp->y))) == 0)
+			;	/* Skip non-numeric input. */
+		if (++nknots >= maxpts) {		/* Too many points, allocate more memory. */
+			if ((data = realloc(data, (maxpts *= 2)*sizeof(struct pair))) == NULL) {
+				error(MEMERR, (char *)maxpts);
+			}
+			dp = data + nknots;
+			sprintf(msg, "more than %d input points, more memory allocated\n",
+			   maxpts/2);
+			error(XTRAPTS,msg);
+		}
+		else
+			++dp;
+	} while (n == 1);
+	--nknots;
+	if (aflag) {	/* Compute maximum ordinate. */
+		max = xmin + (nknots-1)*step;
+	}
+	if (rflag < 2) {
+		xmax = max;
+		if (rflag < 1)
+			xmin = min;
+	}
+}
+
+/* calcspline -- Calculate coeficients of spline functions. */
+calcspline()
+{
+	calccoef();
+	solvband();
+}
+
+/* calccoef -- Calculate coefficients of linear equations determining spline functions. */
+calccoef()
+{	int i, j;
+	struct bandm *ptr, tmp1, tmp2;
+	double h, h1;
+	char str[10];
+
+	ptr = eqns;
+	/* Initialize first row of matrix. */
+	if ((h1 = data[1].x - data[0].x) == 0.) {
+		sprintf(str, "%8.5g", data[0].x);
+		error(DUPX, str);
+	}
+	switch(bdcode) {	/* First equation depends on boundary conditions. */
+	case EXTRAP:
+		ptr->upper = -k;
+		ptr->diag = 1.;
+		ptr->rhs = 0.;
+		break;
+	case FDERIV:		/* Fixed first derivatives at ends. */
+		ptr->upper = 1.;
+		ptr->diag = 2.;
+		h = data[1].x - data[0].x;
+		ptr->rhs = 6.*((data[1].y - data[0].y)/(h*h) - leftd/h);
+		break;
+	case SDERIV:
+		ptr->upper = 0.;
+		ptr->diag = 1.;
+		ptr->rhs = leftd;
+		break;
+	case PERIOD:		/* Periodic splines. */
+		ptr->upper = 1.;
+		ptr->diag = 2.;
+		break;
+	default:
+		error(BADBC, (char *) bdcode);
+	}
+	++ptr;
+
+	/* Initialize rows 1 to n-1, sub-diagonal band is assumed to be 1. */
+	for (i=1; i 	< nknots-1; ++i, ++ptr) {
+		h = h1;
+		if ((h1 = data[i+1].x - data[i].x) == 0.) {
+			sprintf(str, "%8.5g", data[i].x);
+			error(DUPX, str);
+		}
+		ptr->diag = 2.*(h + h1)/h;
+		ptr->upper = h1/h;
+		ptr->rhs = 6.*((data[i+1].y-data[i].y)/(h*h1) -
+		   (data[i].y - data[i-1].y)/(h*h));
+	}
+	/* Initialize last row. */
+	switch(bdcode) {
+	case EXTRAP:		/* Standard case. */
+		ptr->diag = 1.;
+		ptr->upper = -k;	/* Upper holds actual sub-diagonal value. */
+		ptr->rhs = 0.;
+		break;
+	case FDERIV:		/* Fixed first derivatives at ends. */
+		ptr->upper = 1.;
+		ptr->diag = 2.;
+		h = data[nknots-1].x - data[nknots-2].x;
+		ptr->rhs = 6.*((data[nknots-2].y - data[nknots-1].y)/(h*h) + rightd/h);
+		break;
+	case SDERIV:
+		ptr->upper = 0.;
+		ptr->diag = 1.;
+		ptr->rhs = rightd;
+		break;
+	case PERIOD:	/* Use periodic boundary conditions. */
+		/* First and last row are not in tridiagonal form. */
+		h = data[1].x - data[0].x;
+		h1 = data[nknots-1].x - data[nknots-2].x;
+		ptr->diag = 1.;
+		ptr->upper = 0.;
+		tmp1.diag = -1.;
+		tmp1.upper = 0;
+		tmp1.rhs = 0.;
+		tmp2.diag = 2.*h1/h;
+		tmp2.upper = h1/h;
+		tmp2.rhs = 6.*((data[1].y - data[0].y)/(h*h) -
+			   (data[nknots-1].y - data[nknots-2].y)/(h1*h));
+		/* Transform periodic boundary equations to tri-diagonal form. */
+		for (i = 1; i < nknots - 1; ++i) {
+			tmp1.upper /= tmp1.diag;
+			tmp1.rhs /= tmp1.diag;
+			ptr->diag /= tmp1.diag;
+			ptr->upper /= tmp1.diag;
+			tmp1.diag = tmp1.upper - eqns[i].diag;
+			tmp1.upper = -eqns[i].upper;
+			tmp1.rhs -= eqns[i].rhs;
+			tmp2.upper /= tmp2.diag;
+			tmp2.rhs /= tmp2.diag;
+			eqns->diag /= tmp2.diag;
+			eqns->upper /= tmp2.diag;
+			tmp2.diag = tmp2.upper - eqns[nknots-1-i].diag/eqns[nknots-1-i].upper;
+			tmp2.upper = -1./eqns[nknots-1-i].upper;
+			tmp2.rhs -= eqns[nknots-1-i].rhs/eqns[nknots-1-i].upper;
+		}
+		/* Add in remaining terms of boundary condition equation. */
+		ptr->upper += tmp1.diag;
+		ptr->diag += tmp1.upper;
+		ptr->rhs = tmp1.rhs;
+		eqns->diag += tmp2.upper;
+		eqns->upper += tmp2.diag;
+		eqns->rhs = tmp2.rhs;
+		break;
+	default:
+		error(BADBC, (char *) bdcode);
+	}
+}
+
+/* solvband -- Solve band matrix for spline functions. */
+solvband()
+{	int i, flag;
+	struct bandm *ptr;
+	double k1;
+	double fabs();
+	int fcompare();
+
+	ptr = eqns;
+	flag = FALSE;
+	/* Make a pass to triangularize matrix. */
+	for (i=1; i < nknots - 1; ++i, ++ptr) {
+		if (fabs(ptr->diag) < EPS) {
+		/* Near zero on diagonal, pivot. */
+			if (fabs(ptr->upper) < EPS)
+				error(SINGERR, "");
+			flag = TRUE;
+			ptr->diag = i;		/* Keep row index in diag.
+						Actual value of diag is always 1. */
+			if (i == nknots - 2) {
+				flag = 2;
+				/* Exchange next to last and last rows. */
+				k1 = ptr->rhs/ptr->upper;
+				if (fabs((ptr+1)->upper) < EPS)
+					error(SINGERR, "");
+				ptr->rhs = (ptr+1)->rhs/(ptr+1)->upper;
+				ptr->upper = (ptr+1)->diag/(ptr+1)->upper;
+				(ptr+1)->upper = 0.;
+				(ptr+1)->rhs = k1;
+			}
+			else {
+				ptr->rhs = (ptr+1)->rhs - (k1 = ptr->rhs/ptr->upper)*(ptr+1)->diag;
+				ptr->upper = (ptr+1)->upper;	/* This isn't super-diagonal element
+								but rather one to its right. 
+								Must watch for this when 
+								back substituting. */
+				(++ptr)->diag = i-1;
+				++i;
+				ptr->upper = 0;
+				ptr->rhs = k1;
+				(ptr+1)->rhs -= ptr->rhs;
+			}
+		}
+		else {
+			ptr->upper /= ptr->diag;
+			ptr->rhs /= ptr->diag;
+			ptr->diag = i-1;		/* Used to reorder solution if needed. */
+			(ptr+1)->diag -= ptr->upper;
+			(ptr+1)->rhs -= ptr->rhs;
+		}
+	}
+	/* Last row is a special case. */
+	if (flag != 2) {
+		/* If flag == 2 last row is already computed. */
+		ptr->upper /= ptr->diag;
+		ptr->rhs /= ptr->diag;
+		ptr->diag = ptr - eqns;
+		++ptr;
+		ptr->diag -= (ptr-1)->upper * ptr->upper;
+		ptr->rhs -= (ptr-1)->rhs * ptr->upper;
+		ptr->rhs /= ptr->diag;
+		ptr->diag = ptr - eqns;
+	}
+
+	/* Now make a pass back substituting for solution. */
+	--ptr;
+	for ( ; ptr >= eqns; --ptr) {
+		if ((int)ptr->diag != ptr - eqns) {
+			/* This row and one above have been exchanged in a pivot. */
+			--ptr;
+			ptr->rhs -= ptr->upper * (ptr+2)->rhs;
+		}
+		else
+			ptr->rhs -= ptr->upper * (ptr+1)->rhs;
+	}
+	if (flag) {	/* Undo reordering done by pivoting. */
+		qsort(eqns, nknots, sizeof(struct bandm), fcompare);
+	}
+}
+
+/* fcompare -- Compare two floating point numbers. */
+fcompare(arg1, arg2)
+real *arg1, *arg2;
+{
+	if (arg1 > arg2)
+		return(1);
+	else if (arg1 < arg2)
+		return(-1);
+	else
+		return(0);
+}
+
+/* interpolate -- Do spline interpolation. */
+interpolate()
+{	int i;
+	struct pair *dp;
+	struct bandm *ep;
+	real h, xi, yi, hi, xu, xl, limit;
+	
+	h = (xmax - xmin)/nint;
+	ep = eqns;
+	dp = data + 1;
+	for (xi = xmin; xi < xmax + 0.25*h; xi += h) {
+		while (dp->x < xi && dp < data + nknots - 1) {
+			++dp;	/* Skip to correct interval. */
+			++ep;
+		}
+		if (dp < data + nknots - 1 && dp->x < xmax)
+			limit = dp->x;
+		else
+			limit = xmax;
+		for ( ; xi < limit - 0.25*h; xi += h) {
+			/* Do interpolation. */
+			hi = dp->x - (dp-1)->x;
+			xu = dp->x - xi;
+			xl = xi - (dp-1)->x;
+			yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
+			     (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
+			printf(OFMT, xi, yi);
+		}
+		if (limit != dp->x) {	/* Interpolate. */
+			hi = dp->x - (dp-1)->x;
+			xu = dp->x - xmax;
+			xl = xmax - (dp-1)->x;
+			yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
+			     (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
+			printf(OFMT, xmax, yi);
+		}
+		else {		/* Print knot. */
+			printf(OFMT, dp->x, dp->y);
+			xi = dp->x;
+		}
+	}
+}