From 63d16b9f2eb9a8c006c32709db7b882b9bb667b2 Mon Sep 17 00:00:00 2001 From: pnenzi Date: Mon, 11 Aug 2003 19:38:20 +0000 Subject: [PATCH] Cider simulator (math support routines) Import. --- src/maths/misc/Makefile.am | 14 +++ src/maths/misc/accuracy.c | 177 +++++++++++++++++++++++++++++++++ src/maths/misc/accuracy.h | 14 +++ src/maths/misc/bernoull.c | 87 ++++++++++++++++ src/maths/misc/bernoull.h | 15 +++ src/maths/misc/erfc.c | 47 +++++++++ src/maths/misc/norm.c | 78 +++++++++++++++ src/maths/misc/norm.h | 17 ++++ src/maths/misc/test_accuracy.c | 64 ++++++++++++ src/maths/misc/test_erfc.c | 104 +++++++++++++++++++ 10 files changed, 617 insertions(+) create mode 100644 src/maths/misc/Makefile.am create mode 100644 src/maths/misc/accuracy.c create mode 100644 src/maths/misc/accuracy.h create mode 100644 src/maths/misc/bernoull.c create mode 100644 src/maths/misc/bernoull.h create mode 100644 src/maths/misc/erfc.c create mode 100644 src/maths/misc/norm.c create mode 100644 src/maths/misc/norm.h create mode 100644 src/maths/misc/test_accuracy.c create mode 100644 src/maths/misc/test_erfc.c diff --git a/src/maths/misc/Makefile.am b/src/maths/misc/Makefile.am new file mode 100644 index 000000000..98b6b7d46 --- /dev/null +++ b/src/maths/misc/Makefile.am @@ -0,0 +1,14 @@ +## Process this file with automake to produce Makefile.in + +noinst_LIBRARIES = libmathsmisc.a + +libmathsmisc_a_SOURCES = \ + accuracy.c \ + accuracy.h \ + bernoull.c \ + erfc.c \ + norm.c + +INCLUDES = -I$(top_srcdir)/src/include + +MAINTAINERCLEANFILES = Makefile.in diff --git a/src/maths/misc/accuracy.c b/src/maths/misc/accuracy.c new file mode 100644 index 000000000..73f70389d --- /dev/null +++ b/src/maths/misc/accuracy.c @@ -0,0 +1,177 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group +Modified: 2001 Paolo Nenzi +**********/ + +#include "ngspice.h" + + +/* XXX Globals are hiding here. + * All globals have been moved to main.c were all true globals + * should reside. May be in the future that the machine dependent + * global symbols will be gathered into a structure. + * + * Paolo Nenzi 2002 + */ + +/* + * void + * evalAccLimits(void) + * + * This function computes the floating point accuracy on your machine. + * This part of the code is critical and will affect the result of + * CIDER simulator. CIDER uses directly some constants calculated here + * while SPICE does not (yet), but remember that spice runs on the same + * machine that runs CIDER and uses the same floating point implementation. + * + * A note to x86 Linux users: + * + * Intel processors (from i386 up to the latest Pentiums) do FP calculations + * in two ways: + * 53-bit mantissa mode (64 bits overall) + * 64-bit mantissa mode (80 bits overall) + * + * The 53-bit mantissa mode conforms to the IEEE 754 standard (double + * precision), the other (64-bit mantissa mode) does not conform to + * IEEE 754, butlet programmers to use the higher precision "long double" + * data type. + * + * Now the problem: the x86 FPU can be in on mode only, therefore is + * not possible to have IEEE conformant double and "long double" data + * type at the same time. You have to choose which one you prefer. + * + * Linux developers decided that was better to give "long double" to + * programmers that to provide a standard "double" implementation and, + * by default, Linux set the FPU in 64 bit mantissa mode. FreeBSD, on + * the other side, adopt the opposite solution : the FPU is in 53 bit + * mantissa mode. + * + * Since no one but the programmers really knows what a program requires, + * the final decision on the FPU mode of operation is left to her/him. + * It is possible to alter the FPU mode of operation using instruction + * produced by the operating system. + * + * Thanks to AMAKAWA Shuhei for the information on x86 FPU Linux and + * freeBSD on which the above text was derived. + * Paolo Nenzi 2002. + */ + +void +evalAccLimits(void) +{ + double acc = 1.0; + double xl = 0.0; + double xu = 1.0; + double xh, x1, x2, expLim; + double muLim, temp1, temp2, temp3, temp4; + + double xhold; /* Introduced to avoid numerical trap if + using non IEEE754 FPU */ + + +/* First we compute accuracy */ + + for( ; (acc + 1.0) > 1.0 ; ) { + acc *= 0.5; + } + acc *= 2.0; + Acc = acc; + +/* + * This loop has been modified to include a variable to track + * xh change. If it does not change, we exit the loop. This is + * an ugly countermeasure to avoid infinite cycling when in + * x86 64-bit mantissa mode. + * + * Paolo Nenzi 2002 + */ + + xh = 0.5 * (xl + xu); + for( ; (xu-xl > (2.0 * acc * (xu + xl))); ) { + x1 = 1.0 / ( 1.0 + (0.5 * xh) ); + x2 = xh / ( exp(xh) - 1.0 ); + if( (x1 - x2) <= (acc * (x1 + x2))) { + xl = xh; + xhold = xh; + } else { + xu = xh; + xhold = xh; + } + xh = 0.5 * (xl + xu); + if (xhold == xh) break; + } + BMin = xh; + BMax = -log( acc ); + + +/* + * This loop calculate the maximum exponent x for + * which the function exp(-x) returns a value greater + * than 0. The result is used to prevent underflow + * on large negative arguments to exponential. + * AFAIK: used only in Bernoulli function. + */ + + expLim = 80.0; + for( ; exp( -expLim ) > 0.0; ) { + expLim += 1.0; + } + expLim -= 1.0; + ExpLim = expLim; + +/* + * What this loop does ??? + */ + + muLim = 1.0; + temp4 = 0.0; + for( ; 1.0 - temp4 > acc; ){ + muLim *= 0.5; + temp1 = muLim; + temp2 = pow( temp1 , 0.333 ) ; + temp3 = 1.0 / (1.0 + temp1 * temp2 ); + temp4 = pow( temp3 , 0.37/1.333 ); + } + muLim *= 2.0; + MuLim = muLim; + + muLim = 1.0; + temp3 = 0.0; + for( ; 1.0 - temp3 > acc; ){ + muLim *= 0.5; + temp1 = muLim; + temp2 = 1.0 / (1.0 + temp1 * temp1 ); + temp3 = sqrt( temp2 ); + } + muLim *= 2.0; + MutLim = muLim; + + +} + +/* + * Other misterious code. + * Berkeley's people love to leave spare code for info archeologists. + * +main () +{ + double x; + double bx, dbx, bMx, dbMx; + + evalAccLimits(); + for( x = 0.0; x <= 100.0 ; x += 1.0 ) { + bernoulliNew( x, &bx, &dbx, &bMx, &dbMx, 1); + printf( "\nbernoulliNew: x = %e bx = %e dbx = %e bMx = %e dbMx = %e ", x, bx, dbx, bMx, dbMx ); + bernoulli( x, &bx, &dbx, &bMx, &dbMx); + printf( "\nbernoulli: x = %e bx = %e dbx = %e bMx = %e dbMx = %e ", x, bx, dbx, bMx, dbMx ); + } + for( x = 0.0; x >= -100.0 ; x -= 1.0 ) { + bernoulliNew( x, &bx, &dbx, &bMx, &dbMx, 1); + printf( "\nbernoulliNew: x = %e bx = %e dbx = %e bMx = %e dbMx = %e ", x, bx, dbx, bMx, dbMx ); + bernoulli( x, &bx, &dbx, &bMx, &dbMx); + printf( "\nbernoulli: x = %e bx = %e dbx = %e bMx = %e dbMx = %e ", x, bx, dbx, bMx, dbMx ); + } +} + +*/ diff --git a/src/maths/misc/accuracy.h b/src/maths/misc/accuracy.h new file mode 100644 index 000000000..42b83bfbd --- /dev/null +++ b/src/maths/misc/accuracy.h @@ -0,0 +1,14 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +Authors: 1987 Karti Mayaram, 1991 David Gates +**********/ +/* + * Definitions of Globals for Machine Accuracy Limits + */ +#ifndef ACC_H +#define ACC_H +extern double BMin; /* lower limit for B(x) */ +extern double BMax; /* upper limit for B(x) */ +extern double ExpLim; /* limit for exponential */ +extern double Accuracy; /* accuracy of the machine */ +#endif /* ACC_H */ diff --git a/src/maths/misc/bernoull.c b/src/maths/misc/bernoull.c new file mode 100644 index 000000000..4cf74b475 --- /dev/null +++ b/src/maths/misc/bernoull.c @@ -0,0 +1,87 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group +**********/ + +#include "ngspice.h" +#include "accuracy.h" + +/* + * this function computes the bernoulli function + * f(x) = x / ( exp (x) - 1 ) and its derivatives. the function + * f(-x) alongwith its derivative is also computed. + */ + +/* + * input delta-psi + * outputs f(x), df(x)/dx, f(-x), and df(-x)/dx + */ + +void bernoulli (double x, double *pfx, double *pDfxDx, double *pfMx, + double *pDfMxDx, BOOLEAN derivAlso) +{ + double fx, fMx; + double dFxDx = 0.0; + double dFMxDx = 0.0; + double expX, temp; + + if( x <= -BMax ) { + fx = -x; + if( x <= -ExpLim ) { + fMx = 0.0; + if( derivAlso ) { + dFxDx = -1.0; + dFMxDx = 0.0; + } + } + else { + expX = exp( x ); + fMx = -x * expX; + if( derivAlso ) { + dFxDx = fMx - 1.0; + dFMxDx = -expX * ( 1.0 + x ); + } + } + } + else if ( ABS( x) <= BMin ) { + fx = 1.0 / (1.0 + 0.5 * x ); + fMx = 1.0 / (1.0 - 0.5 * x ); + if( derivAlso ) { + temp = 1.0 + x; + dFxDx = -(0.5 + x / 3.0) / temp; + dFMxDx = (0.5 + 2 * x /3.0 )/ temp; + } + } + else if ( x >= BMax ) { + fMx = x; + if( x >= ExpLim ) { + fx = 0.0; + if( derivAlso ) { + dFxDx = 0.0; + dFMxDx = 1.0; + } + } + else { + expX = exp( -x ); + fx = x * expX; + if( derivAlso ) { + dFxDx = expX * ( 1.0 - x ); + dFMxDx = 1.0 - fx; + } + } + } + else { + expX = exp( x ); + temp = 1.0 / ( expX - 1.0 ); + fx = x * temp; + fMx = expX * fx; + if( derivAlso ) { + dFxDx = temp * ( 1.0 - fMx ); + dFMxDx = temp * ( expX - fMx ); + } + } + *pfx = fx; + *pfMx = fMx; + *pDfxDx = dFxDx; + *pDfMxDx = dFMxDx; +} diff --git a/src/maths/misc/bernoull.h b/src/maths/misc/bernoull.h new file mode 100644 index 000000000..79d5823f3 --- /dev/null +++ b/src/maths/misc/bernoull.h @@ -0,0 +1,15 @@ +/********** + (C) Paolo Nenzi 2001 + **********/ + +/* + * Bernoulli function + */ + +#ifndef BERNOULL_H +#define BERNOULL_H + +extern void bernoulli (double, double *, double *, double *, + double *, BOOLEAN); + +#endif /* BERNOULL_H */ diff --git a/src/maths/misc/erfc.c b/src/maths/misc/erfc.c new file mode 100644 index 000000000..8dcb59f4d --- /dev/null +++ b/src/maths/misc/erfc.c @@ -0,0 +1,47 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group +**********/ + +#include "ngspice.h" +#include "numglobs.h" +#include "numconst.h" + +/* erfc computes the erfc(x) the code is from sedan's derfc.f */ + +double erfc ( double x) +{ + double sqrtPi, n, temp1, xSq, sum1, sum2; + sqrtPi = sqrt( PI ); + x = ABS( x ); + n = 1.0; + xSq = 2.0 * x * x; + sum1 = 0.0; + + if ( x > 3.23 ) { + /* asymptotic expansion */ + temp1 = exp( - x * x ) / ( sqrtPi * x ); + sum2 = temp1; + + while ( sum1 != sum2 ) { + sum1 = sum2; + temp1 = -1.0 * ( temp1 / xSq ); + sum2 += temp1; + n += 2.0; + } + return( sum2 ); + } + else { + /* series expansion for small x */ + temp1 = ( 2.0 / sqrtPi ) * exp( - x * x ) * x; + sum2 = temp1; + while ( sum1 != sum2 ) { + n += 2.0; + sum1 = sum2; + temp1 *= xSq / n; + sum2 += temp1; + } + return( 1.0 - sum2 ); + } +} + diff --git a/src/maths/misc/norm.c b/src/maths/misc/norm.c new file mode 100644 index 000000000..6a9576f67 --- /dev/null +++ b/src/maths/misc/norm.c @@ -0,0 +1,78 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group +**********/ + +#include "ngspice.h" + +/* functions to compute max and one norms of a given vector of doubles */ + +double +maxNorm(double *vector, int size) +{ + double norm = 0.0; + double candidate; + int index, nIndex; + + nIndex = 1; + for( index = 1; index <= size; index++ ) { + candidate = fabs(vector[ index ]); + if( candidate > norm ) { + norm = candidate; + nIndex = index; + } + } + /* + printf("\n maxNorm: index = %d", nIndex); + */ + return( norm ); +} + + +double +oneNorm(double *vector, int size) +{ + double norm = 0.0; + double value; + int index; + + for( index = 1; index <= size; index++ ) { + value = vector[ index ]; + if( value < 0.0 ) + norm -= value; + else + norm += value; + } + return( norm ); +} + +double +l2Norm(double *vector, int size) +{ + double norm = 0.0; + double value; + int index; + + for( index = 1; index <= size; index++ ) { + value = vector[ index ]; + norm += value * value; + } + norm = sqrt( norm ); + return( norm ); +} + +/* + * dot(): + * computes dot product of two vectors + */ +double +dot(double *vector1, double *vector2, int size) +{ + register double dot = 0.0; + register int index; + + for( index = 1; index <= size; index++ ) { + dot += vector1[ index ] * vector2[ index ]; + } + return( dot ); +} diff --git a/src/maths/misc/norm.h b/src/maths/misc/norm.h new file mode 100644 index 000000000..fdc496624 --- /dev/null +++ b/src/maths/misc/norm.h @@ -0,0 +1,17 @@ +/********** + (C) Paolo Nenzi 2001 + **********/ + +/* + * Bernoulli function + */ + +#ifndef NORM_H +#define NORM_H + +extern double maxNorm(double *, int); +extern double oneNorm(double *, int); +extern double l2Norm(double *, int); +extern double dot(double *, double *, int); + +#endif /* NORM_H */ diff --git a/src/maths/misc/test_accuracy.c b/src/maths/misc/test_accuracy.c new file mode 100644 index 000000000..659dc9203 --- /dev/null +++ b/src/maths/misc/test_accuracy.c @@ -0,0 +1,64 @@ +/* Paolo Nenzi 2002 - This program tests some machine + * dependent variables. + */ + +/* Nota: + * + * Compilare due volte nel seguente modo: + * + * gcc test_accuracy.c -o test_64accuracy -lm + * gcc -DIEEEDOUBLE test_accuracy.c -o test_53accuracy -lm + * + */ + +#include +#include +#include +#include + +int main (void) +{ + fpu_control_t prec; + double acc=1.0; + double xl = 0.0; + double xu = 1.0; + double xh, x1, x2; + double xhold =0.0; + + +#ifdef IEEEDOUBLE + _FPU_GETCW(prec); + prec &= ~_FPU_EXTENDED; + prec |= _FPU_DOUBLE; + _FPU_SETCW(prec); +#endif + + for( ; (acc + 1.0) > 1.0 ; ) { + acc *= 0.5; + } + acc *= 2.0; + printf("Accuracy: %e\n", acc); + printf("------------------------------------------------------------------\n"); + + xh = 0.5 * (xl + xu); + + for( ; (xu-xl > (2.0 * acc * (xu + xl))); ) { + + x1 = 1.0 / ( 1.0 + (0.5 * xh) ); + x2 = xh / ( exp(xh) - 1.0 ); + + if( (x1 - x2) <= (acc * (x1 + x2))) { + xl = xh; + xhold = xh; + } else { + xu = xh; + xhold = xh; + } + xh = 0.5 * (xl + xu); +// if (xhold == xh) break; + } +printf("xu-xl: %e \t cond: %e \t xh: %e\n", (xu-xl), (2.0 * acc * (xu + xl)), xh); + +exit(1); + +} diff --git a/src/maths/misc/test_erfc.c b/src/maths/misc/test_erfc.c new file mode 100644 index 000000000..dda072e84 --- /dev/null +++ b/src/maths/misc/test_erfc.c @@ -0,0 +1,104 @@ +/* Paolo Nenzi 2002 - This program tests function + * implementations. + */ + + +/* + * The situation on erfc functions in spice/cider: + * + * First we have the ierfc in spice, a sort of interpolation, which is + * fast to compute but gives is not so "good" + * Then we have derfc from cider, which is accurate but slow + * then we have glibc implementation. + * + * Proposal: + * + * Use Linux implementation as defualt and then test cider one. + */ + + +#include +#include +#include +#include + + + + +double +derfc ( double x) +{ + double sqrtPi, n, temp1, xSq, sum1, sum2; + sqrtPi = sqrt( M_PI ); + x = fabs( x ); + n = 1.0; + xSq = 2.0 * x * x; + sum1 = 0.0; + + if ( x > 3.23 ) { + /* asymptotic expansion */ + temp1 = exp( - x * x ) / ( sqrtPi * x ); + sum2 = temp1; + + while ( sum1 != sum2 ) { + sum1 = sum2; + temp1 = -1.0 * ( temp1 / xSq ); + sum2 += temp1; + n += 2.0; + } + return( sum2 ); + } + else { + /* series expansion for small x */ + temp1 = ( 2.0 / sqrtPi ) * exp( - x * x ) * x; + sum2 = temp1; + while ( sum1 != sum2 ) { + n += 2.0; + sum1 = sum2; + temp1 *= xSq / n; + sum2 += temp1; + } + return( 1.0 - sum2 ); + } +} + +double +ierfc(double x) + +{ + double t, z; + + t = 1/(1 + 0.3275911*x); + z = 1.061405429; + z = -1.453152027 + t * z; + z = 1.421413741 + t * z; + z = -0.284496736 + t * z; + z = 0.254829592 + t * z; + z = exp(-x*x) * t * z; + + return(z); +} + + + +int main (void) +{ + fpu_control_t prec; + double x = -100.0; + double y1= 0.0, y2 = 0.0; + +// _FPU_GETCW(prec); +// prec &= ~_FPU_EXTENDED; +// prec |= _FPU_DOUBLE; +// _FPU_SETCW(prec); + + +for (;(x <= 100.0);) +{ + y1 = ierfc(x); + y2 = derfc(x); + printf("A: %e \t s: %e \t c: %e \t (c-lin): %e\n", x, y1, y2, y2-erfc(x) ); + x = x + 0.1; + } + exit(1); +}