From 5492f4cb890460e54448d0601dd734bfaac2aaa6 Mon Sep 17 00:00:00 2001 From: arno Date: Sat, 13 May 2000 10:56:58 +0000 Subject: [PATCH] * src/frontend/interp.c: Broken up into individual files and moved into their own subdirectory: src/maths/poly. * src/maths/poly/.cvsignore src/maths/poly/Makefile.am src/maths/poly/interpolate.c src/maths/poly/interpolate.h src/maths/poly/poly.h src/maths/poly/polyderiv.c src/maths/poly/polyderiv.h src/maths/poly/polyeval.c src/maths/poly/polyeval.h src/maths/poly/polyfit.c src/maths/poly/polyfit.h: The resulting files. * src/Makefile.am src/maths/Makefile.am: Updates for the new library. --- src/Makefile.am | 2 + src/frontend/interp.c | 274 ----------------------------------- src/maths/Makefile.am | 2 +- src/maths/poly/.cvsignore | 3 + src/maths/poly/Makefile.am | 18 +++ src/maths/poly/interpolate.c | 122 ++++++++++++++++ src/maths/poly/interpolate.h | 8 + src/maths/poly/poly.h | 12 ++ src/maths/poly/polyderiv.c | 11 ++ src/maths/poly/polyderiv.h | 6 + src/maths/poly/polyeval.c | 20 +++ src/maths/poly/polyeval.h | 6 + src/maths/poly/polyfit.c | 116 +++++++++++++++ src/maths/poly/polyfit.h | 9 ++ 14 files changed, 334 insertions(+), 275 deletions(-) create mode 100644 src/maths/poly/.cvsignore create mode 100644 src/maths/poly/Makefile.am create mode 100644 src/maths/poly/interpolate.c create mode 100644 src/maths/poly/interpolate.h create mode 100644 src/maths/poly/poly.h create mode 100644 src/maths/poly/polyderiv.c create mode 100644 src/maths/poly/polyderiv.h create mode 100644 src/maths/poly/polyeval.c create mode 100644 src/maths/poly/polyeval.h create mode 100644 src/maths/poly/polyfit.c create mode 100644 src/maths/poly/polyfit.h diff --git a/src/Makefile.am b/src/Makefile.am index 9da736ad2..1da2f0ce4 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -75,6 +75,7 @@ ngspice_LDADD = \ hlp/libhlp.a \ circuit/libinp.a \ maths/cmaths/libcmaths.a \ + maths/poly/libpoly.a \ maths/ni/libni.a \ maths/sparse/libsparse.a \ misc/libmisc.a @@ -97,6 +98,7 @@ nutmeg_LDADD = \ parser/libparser.a \ hlp/libhlp.a \ maths/cmaths/libcmaths.a \ + maths/poly/libpoly.a \ misc/libmisc.a diff --git a/src/frontend/interp.c b/src/frontend/interp.c index 6ecb8eb1a..d24177ae4 100644 --- a/src/frontend/interp.c +++ b/src/frontend/interp.c @@ -14,271 +14,6 @@ Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group #include "interp.h" -/* Interpolate data from oscale to nscale. data is assumed to be olen long, - * ndata will be nlen long. Returns FALSE if the scales are too strange - * to deal with. Note that we are guaranteed that either both scales are - * strictly increasing or both are strictly decreasing. - */ - - -/* static declarations */ -static int putinterval(double *poly, int degree, double *nvec, int last, double *nscale, - int nlen, double oval, int sign); - - -bool -ft_interpolate(double *data, double *ndata, double *oscale, int olen, double *nscale, int nlen, int degree) -{ - double *result, *scratch, *xdata, *ydata; - int sign, lastone, i, l; - - if ((olen < 2) || (nlen < 2)) { - fprintf(cp_err, "Error: lengths too small to interpolate.\n"); - return (FALSE); - } - if ((degree < 1) || (degree > olen)) { - fprintf(cp_err, "Error: degree is %d, can't interpolate.\n", - degree); - return (FALSE); - } - - if (oscale[1] < oscale[0]) - sign = -1; - else - sign = 1; - - scratch = (double *) tmalloc((degree + 1) * (degree + 2) * - sizeof (double)); - result = (double *) tmalloc((degree + 1) * sizeof (double)); - xdata = (double *) tmalloc((degree + 1) * sizeof (double)); - ydata = (double *) tmalloc((degree + 1) * sizeof (double)); - - /* Deal with the first degree pieces. */ - bcopy((char *) data, (char *) ydata, (degree + 1) * sizeof (double)); - bcopy((char *) oscale, (char *) xdata, (degree + 1) * sizeof (double)); - - while (!ft_polyfit(xdata, ydata, result, degree, scratch)) { - /* If it doesn't work this time, bump the interpolation - * degree down by one. - */ - - if (--degree == 0) { - fprintf(cp_err, "ft_interpolate: Internal Error.\n"); - return (FALSE); - } - - } - - /* Add this part of the curve. What we do is evaluate the polynomial - * at those points between the last one and the one that is greatest, - * without being greater than the leftmost old scale point, or least - * if the scale is decreasing at the end of the interval we are looking - * at. - */ - lastone = -1; - for (i = 0; i < degree; i++) { - lastone = putinterval(result, degree, ndata, lastone, - nscale, nlen, xdata[i], sign); - } - - /* Now plot the rest, piece by piece. l is the - * last element under consideration. - */ - for (l = degree + 1; l < olen; l++) { - - /* Shift the old stuff by one and get another value. */ - for (i = 0; i < degree; i++) { - xdata[i] = xdata[i + 1]; - ydata[i] = ydata[i + 1]; - } - ydata[i] = data[l]; - xdata[i] = oscale[l]; - - while (!ft_polyfit(xdata, ydata, result, degree, scratch)) { - if (--degree == 0) { - fprintf(cp_err, - "interpolate: Internal Error.\n"); - return (FALSE); - } - } - lastone = putinterval(result, degree, ndata, lastone, - nscale, nlen, xdata[i], sign); - } - if (lastone < nlen - 1) /* ??? */ - ndata[nlen - 1] = data[olen - 1]; - tfree(scratch); - tfree(xdata); - tfree(ydata); - tfree(result); - return (TRUE); -} - -/* Takes n = (degree+1) doubles, and fills in result with the n coefficients - * of the polynomial that will fit them. It also takes a pointer to an - * array of n ^ 2 + n doubles to use for scratch -- we want to make this - * fast and avoid doing mallocs for each call. - */ - -bool -ft_polyfit(double *xdata, double *ydata, double *result, int degree, double *scratch) -{ - register double *mat1 = scratch; - register int l, k, j, i; - register int n = degree + 1; - register double *mat2 = scratch + n * n; /* XXX These guys are hacks! */ - double d; - -/* -fprintf(cp_err, "n = %d, xdata = ( ", n); - for (i = 0; i < n; i++) - fprintf(cp_err, "%G ", xdata[i]); - fprintf(cp_err, ")\n"); - fprintf(cp_err, "ydata = ( "); - for (i = 0; i < n; i++) - fprintf(cp_err, "%G ", ydata[i]); - fprintf(cp_err, ")\n"); -*/ - - bzero((char *) result, n * sizeof(double)); - bzero((char *) mat1, n * n * sizeof (double)); - bcopy((char *) ydata, (char *) mat2, n * sizeof (double)); - - /* Fill in the matrix with x^k for 0 <= k <= degree for each point */ - l = 0; - for (i = 0; i < n; i++) { - d = 1.0; - for (j = 0; j < n; j++) { - mat1[l] = d; - d *= xdata[i]; - l += 1; - } - } - - /* Do Gauss-Jordan elimination on mat1. */ - for (i = 0; i < n; i++) { - int lindex; - double largest; - /* choose largest pivot */ - for (j=i, largest = mat1[i * n + i], lindex = i; j < n; j++) { - if (fabs(mat1[j * n + i]) > largest) { - largest = fabs(mat1[j * n + i]); - lindex = j; - } - } - if (lindex != i) { - /* swap rows i and lindex */ - for (k = 0; k < n; k++) { - d = mat1[i * n + k]; - mat1[i * n + k] = mat1[lindex * n + k]; - mat1[lindex * n + k] = d; - } - d = mat2[i]; - mat2[i] = mat2[lindex]; - mat2[lindex] = d; - } - /* Make sure we have a non-zero pivot. */ - if (mat1[i * n + i] == 0.0) { - /* this should be rotated. */ - return (FALSE); - } - for (j = i + 1; j < n; j++) { - d = mat1[j * n + i] / mat1[i * n + i]; - for (k = 0; k < n; k++) - mat1[j * n + k] -= d * mat1[i * n + k]; - mat2[j] -= d * mat2[i]; - } - } - - for (i = n - 1; i > 0; i--) - for (j = i - 1; j >= 0; j--) { - d = mat1[j * n + i] / mat1[i * n + i]; - for (k = 0; k < n; k++) - mat1[j * n + k] -= - d * mat1[i * n + k]; - mat2[j] -= d * mat2[i]; - } - - /* Now write the stuff into the result vector. */ - for (i = 0; i < n; i++) { - result[i] = mat2[i] / mat1[i * n + i]; - /* printf(cp_err, "result[%d] = %G\n", i, result[i]);*/ - } - -#define ABS_TOL 0.001 -#define REL_TOL 0.001 - - /* Let's check and make sure the coefficients are ok. If they aren't, - * just return FALSE. This is not the best way to do it. - */ - for (i = 0; i < n; i++) { - d = ft_peval(xdata[i], result, degree); - if (fabs(d - ydata[i]) > ABS_TOL) { - /* - fprintf(cp_err, - "Error: polyfit: x = %le, y = %le, int = %le\n", - xdata[i], ydata[i], d); - printmat("mat1", mat1, n, n); - printmat("mat2", mat2, n, 1); - */ - return (FALSE); - } else if (fabs(d - ydata[i]) / (fabs(d) > ABS_TOL ? fabs(d) : - ABS_TOL) > REL_TOL) { - /* - fprintf(cp_err, - "Error: polyfit: x = %le, y = %le, int = %le\n", - xdata[i], ydata[i], d); - printmat("mat1", mat1, n, n); - printmat("mat2", mat2, n, 1); - */ - return (FALSE); - } - } - - return (TRUE); -} - -/* Returns thestrchr of the last element that was calculated. oval is the - * value of the old scale at the end of the interval that is being interpolated - * from, and sign is 1 if the old scale was increasing, and -1 if it was - * decreasing. - */ - -static int -putinterval(double *poly, int degree, double *nvec, int last, double *nscale, int nlen, double oval, int sign) -{ - int end, i; - - /* See how far we have to go. */ - for (end = last + 1; end < nlen; end++) - if (nscale[end] * sign > oval * sign) - break; - end--; - - for (i = last + 1; i <= end; i++) - nvec[i] = ft_peval(nscale[i], poly, degree); - return (end); -} - - -double -ft_peval(double x, double *coeffs, int degree) -{ - double y; - int i; - - if (!coeffs) - return 0.0; /* XXX Should not happen */ - - y = coeffs[degree]; /* there are (degree+1) coeffs */ - - for (i = degree - 1; i >= 0; i--) { - y *= x; - y += coeffs[i]; - } - - return y; -} - void lincopy(struct dvec *ov, double *newscale, int newlen, struct dvec *oldscale) { @@ -311,12 +46,3 @@ lincopy(struct dvec *ov, double *newscale, int newlen, struct dvec *oldscale) return; } -void -ft_polyderiv(double *coeffs, int degree) -{ - int i; - - for (i = 0; i < degree; i++) { - coeffs[i] = (i + 1) * coeffs[i + 1]; - } -} diff --git a/src/maths/Makefile.am b/src/maths/Makefile.am index 793998744..e9f72164b 100644 --- a/src/maths/Makefile.am +++ b/src/maths/Makefile.am @@ -1,4 +1,4 @@ # Process this file with automake -SUBDIRS = cmaths ni sparse +SUBDIRS = cmaths ni sparse poly MAINTAINERCLEANFILES = Makefile.in diff --git a/src/maths/poly/.cvsignore b/src/maths/poly/.cvsignore new file mode 100644 index 000000000..e440fafda --- /dev/null +++ b/src/maths/poly/.cvsignore @@ -0,0 +1,3 @@ +Makefile.in +Makefile +.deps diff --git a/src/maths/poly/Makefile.am b/src/maths/poly/Makefile.am new file mode 100644 index 000000000..bdff80280 --- /dev/null +++ b/src/maths/poly/Makefile.am @@ -0,0 +1,18 @@ +## Process this file with automake to produce Makefile.in + +noinst_LIBRARIES = libpoly.a + +libpoly_a_SOURCES = \ + interpolate.c \ + interpolate.h \ + polyfit.c \ + polyfit.h \ + polyderiv.c \ + polyderiv.h \ + polyeval.c \ + polyeval.h + + +INCLUDES = -I$(top_srcdir)/src/include + +MAINTAINERCLEANFILES = Makefile.in diff --git a/src/maths/poly/interpolate.c b/src/maths/poly/interpolate.c new file mode 100644 index 000000000..b9482125f --- /dev/null +++ b/src/maths/poly/interpolate.c @@ -0,0 +1,122 @@ +#include +#include +#include + +#include "interpolate.h" +#include "polyeval.h" +#include "polyfit.h" + +/* Returns thestrchr of the last element that was calculated. oval is + * the value of the old scale at the end of the interval that is being + * interpolated from, and sign is 1 if the old scale was increasing, + * and -1 if it was decreasing. */ +static int +putinterval(double *poly, int degree, double *nvec, + int last, double *nscale, int nlen, double oval, int sign) +{ + int end, i; + + /* See how far we have to go. */ + for (end = last + 1; end < nlen; end++) + if (nscale[end] * sign > oval * sign) + break; + end--; + + for (i = last + 1; i <= end; i++) + nvec[i] = ft_peval(nscale[i], poly, degree); + return (end); +} + + +/* Interpolate data from oscale to nscale. data is assumed to be olen long, + * ndata will be nlen long. Returns FALSE if the scales are too strange + * to deal with. Note that we are guaranteed that either both scales are + * strictly increasing or both are strictly decreasing. + */ +bool +ft_interpolate(double *data, double *ndata, double *oscale, int olen, + double *nscale, int nlen, int degree) +{ + double *result, *scratch, *xdata, *ydata; + int sign, lastone, i, l; + + if ((olen < 2) || (nlen < 2)) { + fprintf(cp_err, "Error: lengths too small to interpolate.\n"); + return (FALSE); + } + if ((degree < 1) || (degree > olen)) { + fprintf(cp_err, "Error: degree is %d, can't interpolate.\n", + degree); + return (FALSE); + } + + if (oscale[1] < oscale[0]) + sign = -1; + else + sign = 1; + + scratch = (double *) tmalloc((degree + 1) * (degree + 2) * + sizeof (double)); + result = (double *) tmalloc((degree + 1) * sizeof (double)); + xdata = (double *) tmalloc((degree + 1) * sizeof (double)); + ydata = (double *) tmalloc((degree + 1) * sizeof (double)); + + /* Deal with the first degree pieces. */ + bcopy((char *) data, (char *) ydata, (degree + 1) * sizeof (double)); + bcopy((char *) oscale, (char *) xdata, (degree + 1) * sizeof (double)); + + while (!ft_polyfit(xdata, ydata, result, degree, scratch)) { + /* If it doesn't work this time, bump the interpolation + * degree down by one. + */ + + if (--degree == 0) { + fprintf(cp_err, "ft_interpolate: Internal Error.\n"); + return (FALSE); + } + + } + + /* Add this part of the curve. What we do is evaluate the polynomial + * at those points between the last one and the one that is greatest, + * without being greater than the leftmost old scale point, or least + * if the scale is decreasing at the end of the interval we are looking + * at. + */ + lastone = -1; + for (i = 0; i < degree; i++) { + lastone = putinterval(result, degree, ndata, lastone, + nscale, nlen, xdata[i], sign); + } + + /* Now plot the rest, piece by piece. l is the + * last element under consideration. + */ + for (l = degree + 1; l < olen; l++) { + + /* Shift the old stuff by one and get another value. */ + for (i = 0; i < degree; i++) { + xdata[i] = xdata[i + 1]; + ydata[i] = ydata[i + 1]; + } + ydata[i] = data[l]; + xdata[i] = oscale[l]; + + while (!ft_polyfit(xdata, ydata, result, degree, scratch)) { + if (--degree == 0) { + fprintf(cp_err, + "interpolate: Internal Error.\n"); + return (FALSE); + } + } + lastone = putinterval(result, degree, ndata, lastone, + nscale, nlen, xdata[i], sign); + } + if (lastone < nlen - 1) /* ??? */ + ndata[nlen - 1] = data[olen - 1]; + tfree(scratch); + tfree(xdata); + tfree(ydata); + tfree(result); + return (TRUE); +} diff --git a/src/maths/poly/interpolate.h b/src/maths/poly/interpolate.h new file mode 100644 index 000000000..f39fd88ff --- /dev/null +++ b/src/maths/poly/interpolate.h @@ -0,0 +1,8 @@ +#ifndef _INTERPOLATE_H +#define _INTERPOLATE_H + +#include + +bool ft_interpolate(double *data, double *ndata, double *oscale, int olen, double *nscale, int nlen, int degree); + +#endif diff --git a/src/maths/poly/poly.h b/src/maths/poly/poly.h new file mode 100644 index 000000000..54235a1fa --- /dev/null +++ b/src/maths/poly/poly.h @@ -0,0 +1,12 @@ +/* The public interface to the polygon library. */ + +#ifndef _POLY_H +#define _POLY_H + +#include "interpolate.h" +#include "polyderiv.h" +#include "polyeval.h" +#include "polyfit.h" + +#endif + diff --git a/src/maths/poly/polyderiv.c b/src/maths/poly/polyderiv.c new file mode 100644 index 000000000..80666a38e --- /dev/null +++ b/src/maths/poly/polyderiv.c @@ -0,0 +1,11 @@ +#include "polyderiv.h" + +void +ft_polyderiv(double *coeffs, int degree) +{ + int i; + + for (i = 0; i < degree; i++) { + coeffs[i] = (i + 1) * coeffs[i + 1]; + } +} diff --git a/src/maths/poly/polyderiv.h b/src/maths/poly/polyderiv.h new file mode 100644 index 000000000..a0a275edf --- /dev/null +++ b/src/maths/poly/polyderiv.h @@ -0,0 +1,6 @@ +#ifndef _POLYDERIV_H +#define _POLYDERIV_H + +void ft_polyderiv(double *coeffs, int degree); + +#endif diff --git a/src/maths/poly/polyeval.c b/src/maths/poly/polyeval.c new file mode 100644 index 000000000..0860cbe03 --- /dev/null +++ b/src/maths/poly/polyeval.c @@ -0,0 +1,20 @@ +#include "polyeval.h" + +double +ft_peval(double x, double *coeffs, int degree) +{ + double y; + int i; + + if (!coeffs) + return 0.0; /* XXX Should not happen */ + + y = coeffs[degree]; /* there are (degree+1) coeffs */ + + for (i = degree - 1; i >= 0; i--) { + y *= x; + y += coeffs[i]; + } + + return y; +} diff --git a/src/maths/poly/polyeval.h b/src/maths/poly/polyeval.h new file mode 100644 index 000000000..1ed9c8074 --- /dev/null +++ b/src/maths/poly/polyeval.h @@ -0,0 +1,6 @@ +#ifndef _POLYEVAL_H +#define _POLYEVAL_H + +double ft_peval(double x, double *coeffs, int degree); + +#endif diff --git a/src/maths/poly/polyfit.c b/src/maths/poly/polyfit.c new file mode 100644 index 000000000..44ba47cfa --- /dev/null +++ b/src/maths/poly/polyfit.c @@ -0,0 +1,116 @@ +#include + +#include "polyfit.h" +#include "polyeval.h" + +/* Takes n = (degree+1) doubles, and fills in result with the n + * coefficients of the polynomial that will fit them. It also takes a + * pointer to an array of n ^ 2 + n doubles to use for scratch -- we + * want to make this fast and avoid doing mallocs for each call. */ +bool +ft_polyfit(double *xdata, double *ydata, double *result, + int degree, double *scratch) +{ + double *mat1 = scratch; + int l, k, j, i; + int n = degree + 1; + double *mat2 = scratch + n * n; /* XXX These guys are hacks! */ + double d; + + memset((char *) result, 0, n * sizeof(double)); + memset((char *) mat1, 0, n * n * sizeof (double)); + memcpy((char *) ydata, (char *) mat2, n * sizeof (double)); + + /* Fill in the matrix with x^k for 0 <= k <= degree for each point */ + l = 0; + for (i = 0; i < n; i++) { + d = 1.0; + for (j = 0; j < n; j++) { + mat1[l] = d; + d *= xdata[i]; + l += 1; + } + } + + /* Do Gauss-Jordan elimination on mat1. */ + for (i = 0; i < n; i++) { + int lindex; + double largest; + /* choose largest pivot */ + for (j=i, largest = mat1[i * n + i], lindex = i; j < n; j++) { + if (fabs(mat1[j * n + i]) > largest) { + largest = fabs(mat1[j * n + i]); + lindex = j; + } + } + if (lindex != i) { + /* swap rows i and lindex */ + for (k = 0; k < n; k++) { + d = mat1[i * n + k]; + mat1[i * n + k] = mat1[lindex * n + k]; + mat1[lindex * n + k] = d; + } + d = mat2[i]; + mat2[i] = mat2[lindex]; + mat2[lindex] = d; + } + /* Make sure we have a non-zero pivot. */ + if (mat1[i * n + i] == 0.0) { + /* this should be rotated. */ + return (FALSE); + } + for (j = i + 1; j < n; j++) { + d = mat1[j * n + i] / mat1[i * n + i]; + for (k = 0; k < n; k++) + mat1[j * n + k] -= d * mat1[i * n + k]; + mat2[j] -= d * mat2[i]; + } + } + + for (i = n - 1; i > 0; i--) + for (j = i - 1; j >= 0; j--) { + d = mat1[j * n + i] / mat1[i * n + i]; + for (k = 0; k < n; k++) + mat1[j * n + k] -= + d * mat1[i * n + k]; + mat2[j] -= d * mat2[i]; + } + + /* Now write the stuff into the result vector. */ + for (i = 0; i < n; i++) { + result[i] = mat2[i] / mat1[i * n + i]; + /* printf(cp_err, "result[%d] = %G\n", i, result[i]);*/ + } + +#define ABS_TOL 0.001 +#define REL_TOL 0.001 + + /* Let's check and make sure the coefficients are ok. If they aren't, + * just return FALSE. This is not the best way to do it. + */ + for (i = 0; i < n; i++) { + d = ft_peval(xdata[i], result, degree); + if (fabs(d - ydata[i]) > ABS_TOL) { + /* + fprintf(cp_err, + "Error: polyfit: x = %le, y = %le, int = %le\n", + xdata[i], ydata[i], d); + printmat("mat1", mat1, n, n); + printmat("mat2", mat2, n, 1); + */ + return (FALSE); + } else if (fabs(d - ydata[i]) / (fabs(d) > ABS_TOL ? fabs(d) : + ABS_TOL) > REL_TOL) { + /* + fprintf(cp_err, + "Error: polyfit: x = %le, y = %le, int = %le\n", + xdata[i], ydata[i], d); + printmat("mat1", mat1, n, n); + printmat("mat2", mat2, n, 1); + */ + return (FALSE); + } + } + + return (TRUE); +} diff --git a/src/maths/poly/polyfit.h b/src/maths/poly/polyfit.h new file mode 100644 index 000000000..e27cf3c49 --- /dev/null +++ b/src/maths/poly/polyfit.h @@ -0,0 +1,9 @@ +#ifndef _POLYFIT_H +#define _POLYFIT_H + +#include + +bool ft_polyfit(double *xdata, double *ydata, double *result, + int degree, double *scratch); + +#endif