182 lines
5.0 KiB
C
182 lines
5.0 KiB
C
/**********
|
|
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, dif; /* Introduced to avoid numerical trap if
|
|
using non IEEE754 FPU */
|
|
|
|
#ifndef CIDER
|
|
double Acc, BMin, BMax, ExpLim, MuLim, MutLim;
|
|
#endif
|
|
|
|
/* 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);
|
|
|
|
dif = fabs(xhold - xh);
|
|
if (dif <= DBL_EPSILON) 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 );
|
|
}
|
|
}
|
|
|
|
*/
|