Cider simulator (math support routines) Import.

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

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 );
+    }
+}
+
+*/

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 */

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;
+}

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 */

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 );
+    }
+}
+

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 );
+}

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 */

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 <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <fpu_control.h>
+
+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);
+
+}

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 <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <fpu_control.h>
+
+
+
+
+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);
+}