ted@adelie.Adelie.COM (Ted Stefanik) (08/09/89)
Here is a program that allows you to fit a set of points to the closest
polynomial. It also has an imbedded matrix inverter/solver, which I couldn't
find elsewhere on my ULTRIX system.
Note that solving for a polynomial above about the fifth order can easily
cause floating point exceptions.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# Makefile
# fit.c
# fit.h
# mat.c
# This archive created: Tue Aug 8 18:49:44 1989
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'Makefile'" '(148 characters)'
if test -f 'Makefile'
then
echo shar: "will not over-write existing file 'Makefile'"
else
sed 's/^X#//' << \SHAR_EOF > 'Makefile'
X#D = -g
X#CFLAGS = $(D) $(FLOAT)
X#FLOAT =
X#MATH = -lm
X#
X#
X#fit: fit.o mat.o
X# cc fit.o mat.o $(D) $(FLOAT) -o fit $(MATH)
X#
X#fit.o: fit.h
X#mat.o: fit.h
SHAR_EOF
fi
echo shar: "extracting 'fit.c'" '(3032 characters)'
if test -f 'fit.c'
then
echo shar: "will not over-write existing file 'fit.c'"
else
sed 's/^X#//' << \SHAR_EOF > 'fit.c'
X##include "fit.h"
X#
X#static void ParseArgs();
X#static void Usage();
X#
X#static void ReadPts();
X#static void PolyFit();
X#static void Results();
X#static void Stats();
X#
X#
X#static points Points;
X#static int MaxPt;
X#static int PtCnt;
X#static int Order;
X#
X#static matrix Mat;
X#static double Det;
X#
X#
X#main(argc, argv)
X# int argc;
X# char **argv;
X#{
X# ParseArgs(argc, argv);
X#
X# ReadPts();
X# if (PtCnt <= 0)
X# Usage();
X# PolyFit();
X# Results();
X#
X# Stats();
X#}
X#
X#
X#static void ParseArgs(argc, argv)
X# int argc;
X# char **argv;
X#{
X# if (argc != 3)
X# Usage(argv[0]);
X#
X# MaxPt = atoi(argv[1]);
X# Order = atoi(argv[2]);
X#
X# if (MaxPt <= 0 || Order <= 0)
X# Usage(argv[0]);
X#}
X#
X#
X#static void Usage(prog)
X# char *prog;
X#{
X# fprintf(stderr,
X# "Usage: %s max-number-of-data-points order-of-polynomial\n", prog);
X# fprintf(stderr,
X# "Note: points are read from stdin (a coordinate point/line)\n");
X# exit(1);
X#}
X#
X#
X#static void ReadPts()
X#{
X# char line[1024];
X#
X# Points = (points)Calloc(MaxPt, sizeof(point));
X#
X# for (PtCnt = 0; PtCnt < MaxPt; PtCnt++)
X# {
X# if (gets(line) == NULL)
X# return;
X# sscanf(line, "%f %f", &Points[PtCnt].X, &Points[PtCnt].Y);
X# }
X#}
X#
X#
X#static void PolyFit()
X#{
X# register int i, p, r, c;
X# register matrow xvec;
X#
X# Mat = MatInit(Order + 1);
X# xvec = (matrow)Calloc(Order * 2 + 1, sizeof(matelm));
X#
X# for (p = 0; p < PtCnt; p++)
X# {
X# for (xvec[0] = 1.0, i = 1; i <= Order * 2; i++)
X# xvec[i] = xvec[i - 1] * Points[p].X;
X# for (r = 0; r < Order + 1; r++)
X# {
X# Mat[r][Order + 1] += Points[p].Y * xvec[r];
X# for (c = 0; c < Order + 1; c++)
X# Mat[r][c] += xvec[r + c];
X# }
X# }
X#
X# Det = MatInv(Order + 1, Mat);
X#}
X#
X#
X#static void Results()
X#{
X# register int i;
X#
X# printf("%5d = Number of data points input\n", PtCnt);
X# printf("%5d = Order of polynomial fit\n\n", Order);
X#
X##ifndef NO_DETERM
X# printf("%e = Determinate of matrix solution\n\n", Det);
X##endif
X#
X# printf("y = %f", Mat[0][Order + 1]);
X# for (i = 1; i < Order + 1; i++)
X# printf(" + %f x^%d", Mat[i][Order + 1], i);
X#}
X#
X#
X#static void Stats()
X#{
X# register int i, j;
X# double x, y;
X# double sum, sum2, minerr, maxerr;
X#
X# sum = sum2 = maxerr = 0;
X# minerr = HUGE;
X#
X# for (i = 0; i < PtCnt; i++)
X# {
X# y = Mat[0][Order + 1];
X# for (x = 1, j = 1; j < Order + 1; j++)
X# {
X# x *= Points[i].X;
X# y += x * Mat[j][Order + 1];
X# }
X#
X# y -= Points[i].Y;
X# if (y < 0)
X# y = -y;
X#
X# sum += y;
X# sum2 += y * y;
X#
X# if (y < minerr)
X# minerr = y;
X# if (y > maxerr)
X# maxerr = y;
X# }
X#
X# printf("\n\nFit statistics:\n");
X#
X# printf("%30.15f = Standard Error of Estimate\n", sqrt(sum2 / (PtCnt - 2)));
X# printf("%30.15f = Average Error\n", sum / PtCnt);
X# printf("%30.15f = Error's Standard Deviation\n",
X# sqrt((PtCnt * sum2 - sum * sum) / PtCnt / (PtCnt - 1)));
X# printf("%30.15f = Minimum Error\n", minerr);
X# printf("%30.15f = Maximum Error\n", maxerr);
X#}
SHAR_EOF
fi
echo shar: "extracting 'fit.h'" '(314 characters)'
if test -f 'fit.h'
then
echo shar: "will not over-write existing file 'fit.h'"
else
sed 's/^X#//' << \SHAR_EOF > 'fit.h'
X##include <stdio.h>
X##include <math.h>
X#
X#struct point
X#{
X# double X;
X# double Y;
X#};
X#
X#typedef struct point point;
X#
X#typedef point *points;
X#
X#typedef double matelm;
X#typedef matelm *matrow;
X#typedef matrow *matrix;
X#
X#typedef char *pointer;
X#
X#
X#extern char *calloc();
X#
X#matrix MatInit();
X#double MatInv();
X#pointer Calloc();
SHAR_EOF
fi
echo shar: "extracting 'mat.c'" '(1330 characters)'
if test -f 'mat.c'
then
echo shar: "will not over-write existing file 'mat.c'"
else
sed 's/^X#//' << \SHAR_EOF > 'mat.c'
X##include "fit.h"
X#
X#matrix MatInit(n)
X# register int n;
X#{
X# register int i;
X# register matrix mat;
X#
X# mat = (matrix)Calloc(n, sizeof(matrow));
X#
X# for (i = 0; i < n; i++)
X# mat[i] = (matrow)Calloc(n + 1, sizeof(matelm));
X#
X# return (mat);
X#}
X#
X#
X#double MatInv(n, mat)
X# register int n;
X# register matrix mat;
X#{
X# int pivRow;
X# double piv, det, div, pivTmp;
X# matrow tmpRow;
X# register int p, r, c;
X#
X# for (det = 1.0, p = 0; p < n; p++)
X# {
X# for (piv = mat[pivRow = p][p], r = p; r < n; r++)
X# if (((pivTmp = mat[r][p]) < 0 ? -pivTmp : pivTmp) > piv)
X# piv = mat[pivRow = r][p];
X#
X##ifndef NO_DETERM
X# if ((det *= piv) == 0.0)
X# break;
X##endif
X#
X# if (pivRow != p)
X# {
X##ifndef NO_DETERM
X# det = - det;
X##endif
X# tmpRow = mat[p];
X# mat[p] = mat[pivRow];
X# mat[pivRow] = tmpRow;
X# }
X#
X# for (c = p; c <= n; c++)
X# mat[p][c] /= piv;
X#
X# for (r = 0; r < n; r++)
X# if (r != p)
X# for (c = n; c >= p; c--)
X# mat[r][c] -= mat[r][p] * mat[p][c];
X# }
X#
X##ifndef NO_DETERM
X# return (det);
X##endif
X#}
X#
X#
X#pointer Calloc(n, s)
X# unsigned n, s;
X#{
X# register pointer ret;
X#
X# if ((ret = (pointer)calloc(n, s)) == NULL)
X# {
X# perror("Allocating matrix");
X# exit(1);
X# }
X#
X# return (ret);
X#}
SHAR_EOF
fi
exit 0
# End of shell archive
--
LIVE: Ted Stefanik, (617) 965-8480 x14
USPS: Adelie Corporation, 288 Walnut St., Newtonville, MA 02160
UUCP: ..!{wanginst!infinet | decvax!cca}!emacs!adelie!ted
ARPA: emacs!adelie!ted-unix.ARPA