diff options
author | Jim Meyering <jim@meyering.net> | 2003-11-04 14:42:54 +0000 |
---|---|---|
committer | Jim Meyering <jim@meyering.net> | 2003-11-04 14:42:54 +0000 |
commit | 4f081364fba6b4af6450d088bfda5ce070a3cb38 (patch) | |
tree | 46a5580c66d9941bbe41be4458b8c58efdcb7f22 | |
parent | 0dd12bdc0334f550c4ca45356511aea107f28d0b (diff) | |
download | coreutils-4f081364fba6b4af6450d088bfda5ce070a3cb38.tar.xz |
.
-rw-r--r-- | src/spline.c | 735 |
1 files changed, 0 insertions, 735 deletions
diff --git a/src/spline.c b/src/spline.c deleted file mode 100644 index 0189dec8f..000000000 --- a/src/spline.c +++ /dev/null @@ -1,735 +0,0 @@ -/* 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 (is_number(s)) { - while (is_number(s)) - *str++ = *s++; - *str = '\0'; - type = NUMBER; - } - else { - strcpy(str, s); - s = strchr(s, '\0'); - type = OTHER; - } - *locp = s; - return(type); -} - -/* is_number -- Return TRUE if argument is the ASCII representation of a number. */ -is_number(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; - } - } -} |