diff --git a/ChangeLog b/ChangeLog index 88a4a6573..f5b798e1a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +2008-05-11 Dietmar Warning + * src/misc/missing_math.c,h, src/math/misc/*.*: move math function to one place + * src/include/ngspice.h, src/main.c, src/Makefile.am, src/maths/Makefile.am, + src/maths/misc/Makefile.am, src/misc/Makefile.am, configure.in: organization of libmathmisc.a + * src/maths/misc/erfc.c: better erfc for lossy transmission line + 2008-05-10 Holger Vogt * src/frontend/resource.c: Memory information is now stemming from the /proc file system (LINUX) or using GlobalMemoryStatusEx and diff --git a/configure.in b/configure.in index 0310f965c..4da96f19c 100644 --- a/configure.in +++ b/configure.in @@ -630,10 +630,6 @@ if test "$enable_cider" = "yes"; then $CIDERDIR/input/libciderinput.a \ $CIDERDIR/support/libcidersuprt.a" - CIDERMATH="misc" - - CIDERMATHDIR="maths/misc/libmathmisc.a" - NUMDEV=" spicelib/devices/nbjt/libnbjt.a \ spicelib/devices/nbjt2/libnbjt2.a \ spicelib/devices/numd/libnumd.a \ @@ -650,16 +646,12 @@ CIDERSCRIPTS="devload devaxis ciderinit" else CIDERLIB="" CIDERSIM="" - CIDERMATH="" - CIDERMATHDIR="" NUMDEV="" NUMDEVDIR="" CIDERSCRIPTS="" fi AC_SUBST(CIDERDIR) AC_SUBST(CIDERSIM) -AC_SUBST(CIDERMATH) -AC_SUBST(CIDERMATHDIR) AC_SUBST(NUMDEV) AC_SUBST(NUMDEVDIR) AC_SUBST(CIDERSCRIPTS) diff --git a/src/Makefile.am b/src/Makefile.am index e8f088f87..e54586c65 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -109,7 +109,7 @@ ngspice_LDADD = \ @CIDERSIM@ \ maths/deriv/libderiv.a \ maths/cmaths/libcmaths.a \ - @CIDERMATHDIR@ \ + maths/misc/libmathmisc.a \ maths/poly/libpoly.a \ maths/ni/libni.a \ maths/sparse/libsparse.a \ @@ -142,6 +142,7 @@ ngnutmeg_LDADD = \ @NUMPARAMLIB@ \ frontend/help/libhlp.a \ maths/cmaths/libcmaths.a \ + maths/misc/libmathmisc.a \ maths/poly/libpoly.a \ misc/libmisc.a \ spicelib/parser/libinp.a @@ -173,6 +174,7 @@ ngsconvert_SOURCES = ngsconvert.c ngsconvert_LDADD = \ frontend/libfte.a \ frontend/parser/libparser.a \ + maths/misc/libmathmisc.a \ misc/libmisc.a ## proc2mod: diff --git a/src/include/ngspice.h b/src/include/ngspice.h index c08c7c181..55da9ff95 100644 --- a/src/include/ngspice.h +++ b/src/include/ngspice.h @@ -31,6 +31,7 @@ #include #include +#include "missing_math.h" #ifdef STDC_HEADERS # include diff --git a/src/maths/Makefile.am b/src/maths/Makefile.am index cfda08ed0..86b4017fe 100644 --- a/src/maths/Makefile.am +++ b/src/maths/Makefile.am @@ -1,6 +1,6 @@ ## Process this file with automake -SUBDIRS = cmaths ni sparse poly deriv @CIDERMATH@ +SUBDIRS = cmaths ni sparse poly deriv misc DIST_SUBDIRS = cmaths ni sparse poly deriv misc MAINTAINERCLEANFILES = Makefile.in diff --git a/src/maths/misc/Makefile.am b/src/maths/misc/Makefile.am index a94a9b970..f5815c594 100644 --- a/src/maths/misc/Makefile.am +++ b/src/maths/misc/Makefile.am @@ -8,6 +8,10 @@ libmathmisc_a_SOURCES = \ bernoull.h \ bernoull.c \ erfc.c \ + equality.c \ + isnan.c \ + logb.c \ + scalb.c \ norm.h \ norm.c diff --git a/src/maths/misc/accuracy.c b/src/maths/misc/accuracy.c index cd175ace9..a408c0f6f 100644 --- a/src/maths/misc/accuracy.c +++ b/src/maths/misc/accuracy.c @@ -69,6 +69,9 @@ evalAccLimits(void) 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 */ @@ -149,7 +152,6 @@ evalAccLimits(void) muLim *= 2.0; MutLim = muLim; - } /* diff --git a/src/maths/misc/equality.c b/src/maths/misc/equality.c new file mode 100644 index 000000000..9110bed3a --- /dev/null +++ b/src/maths/misc/equality.c @@ -0,0 +1,62 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +**********/ + +#include +#include "ngspice.h" + +#ifdef _MSC_VER +typedef __int64 long64; +#else +typedef long long long64; +#endif + +#define Abs(x) ((x) < 0 ? -(x) : (x)) + +/* From Bruce Dawson, Comparing floating point numbers, + http://www.cygnus-software.com/papers/comparingfloats/Comparing%20floating%20point%20numbers.htm + Original this function is named AlmostEqual2sComplement but we leave it to AlmostEqualUlps + and can leave the code (measure.c, dctran.c) unchanged. The transformation to the 2's complement + prevent problems around 0.0. + One Ulp is equivalent to a maxRelativeError of between 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000. + Practical: 3 < maxUlps < some hundred's (or thousand's) - depending on numerical requirements. +*/ +bool AlmostEqualUlps(double A, double B, int maxUlps) +{ + long64 aInt, bInt, intDiff; + + if (A == B) + return TRUE; + + /* If not - the entire method can not work */ + assert(sizeof(double) == sizeof(long64)); + + /* Make sure maxUlps is non-negative and small enough that the */ + /* default NAN won't compare as equal to anything. */ + assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024); + aInt = *(long64*)&A; + /* Make aInt lexicographically ordered as a twos-complement int */ + if (aInt < 0) +#ifdef _MSC_VER + aInt = 0x8000000000000000 - aInt; +#else + aInt = 0x8000000000000000LL - aInt; +#endif + bInt = *(long64*)&B; + /* Make bInt lexicographically ordered as a twos-complement int */ + if (bInt < 0) +#ifdef _MSC_VER + bInt = 0x8000000000000000 - bInt; +#else + bInt = 0x8000000000000000LL - bInt; +#endif +#ifdef _MSC_VER + intDiff = Abs(aInt - bInt); +#else + intDiff = llabs(aInt - bInt); +#endif +/* printf("A:%e B:%e aInt:%d bInt:%d diff:%d\n", A, B, aInt, bInt, intDiff); */ + if (intDiff <= maxUlps) + return TRUE; + return FALSE; +} diff --git a/src/maths/misc/isnan.c b/src/maths/misc/isnan.c new file mode 100644 index 000000000..87e0ad110 --- /dev/null +++ b/src/maths/misc/isnan.c @@ -0,0 +1,53 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +**********/ + +#include "ngspice.h" + +#ifndef HAVE_ISNAN + +/* isnan (originally) for SOI devices in MINGW32 hvogt (dev.c) */ +union ieee754_double +{ + double d; + + /* This is the IEEE 754 double-precision format. */ + struct + { + /* Together these comprise the mantissa. */ + unsigned int mantissa1:32; + unsigned int mantissa0:20; + unsigned int exponent:11; + unsigned int negative:1; + } ieee; + struct + { + /* Together these comprise the mantissa. */ + unsigned int mantissa1:32; + unsigned int mantissa0:19; + unsigned int quiet_nan:1; + unsigned int exponent:11; + unsigned int negative:1; + } ieee_nan; +}; + +int +isnan(double value) +{ + union ieee754_double u; + + u.d = value; + + /* IEEE 754 NaN's have the maximum possible + exponent and a nonzero mantissa. */ + return ((u.ieee.exponent & 0x7ff) == 0x7ff && + (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)); + +} + +/* + * end isnan.c + */ +#else /* HAVE_ISNAN */ +int Dummy_Symbol_4; +#endif /* HAVE_ISNAN */ diff --git a/src/maths/misc/logb.c b/src/maths/misc/logb.c new file mode 100644 index 000000000..dd4b46cc6 --- /dev/null +++ b/src/maths/misc/logb.c @@ -0,0 +1,36 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +**********/ + +#include "ngspice.h" + +#ifndef HAVE_LOGB + +double +logb(double x) +{ + double y = 0.0; + + if (x != 0.0) + { + if (x < 0.0) + x = - x; + while (x > 2.0) + { + y += 1.0; + x /= 2.0; + } + while (x < 1.0) + { + y -= 1.0; + x *= 2.0; + } + } + else + y = 0.0; + + return y; +} +#else +int Dummy_Symbol_3; +#endif diff --git a/src/maths/misc/scalb.c b/src/maths/misc/scalb.c new file mode 100644 index 000000000..18c127654 --- /dev/null +++ b/src/maths/misc/scalb.c @@ -0,0 +1,34 @@ +/********** +Copyright 1991 Regents of the University of California. All rights reserved. +**********/ + +#include "ngspice.h" + +#ifndef HAVE_SCALB +# ifdef HAVE_SCALBN +# define scalb scalbn +#else /* Chris Inbody */ + +double +scalb(double x, int n) +{ + double y, z = 1.0, k = 2.0; + + if (n < 0) { + n = -n; + k = 0.5; + } + + if (x != 0.0) + for (y = 1.0; n; n >>= 1) { + y *= k; + if (n & 1) + z *= y; + } + + return x * z; +} +# endif /* HAVE_SCALBN */ +#else /* HAVE_SCALB */ +int Dummy_Symbol_1; +#endif /* HAVE_SCALB */ diff --git a/src/misc/Makefile.am b/src/misc/Makefile.am index d958c5e00..831d3d0c3 100644 --- a/src/misc/Makefile.am +++ b/src/misc/Makefile.am @@ -14,8 +14,6 @@ libmisc_a_SOURCES = \ dup2.h \ ivars.c \ ivars.h \ - missing_math.c \ - missing_math.h \ mktemp.c \ mktemp.h \ printnum.c \ diff --git a/src/misc/missing_math.c b/src/misc/missing_math.c deleted file mode 100644 index 4ec08b51f..000000000 --- a/src/misc/missing_math.c +++ /dev/null @@ -1,203 +0,0 @@ -/********** -Copyright 1991 Regents of the University of California. All rights reserved. -$Id$ -**********/ - -/* - * Missing math functions - */ -#include -#include "config.h" -#include "ngspice.h" -#include "missing_math.h" - -#ifdef _MSC_VER -typedef __int64 long64; -#else -typedef long long long64; -#endif - -#define Abs(x) ((x) < 0 ? -(x) : (x)) - -/* From Bruce Dawson, Comparing floating point numbers, - http://www.cygnus-software.com/papers/comparingfloats/Comparing%20floating%20point%20numbers.htm - Original this function is named AlmostEqual2sComplement but we leave it to AlmostEqualUlps - and can leave the code (measure.c, dctran.c) unchanged. The transformation to the 2's complement - prevent problems around 0.0. - One Ulp is equivalent to a maxRelativeError of between 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000. - Practical: 3 < maxUlps < some hundred's (or thousand's) - depending on numerical requirements. -*/ -bool AlmostEqualUlps(double A, double B, int maxUlps) -{ - long64 aInt, bInt, intDiff; - - if (A == B) - return TRUE; - - /* If not - the entire method can not work */ - assert(sizeof(double) == sizeof(long64)); - - /* Make sure maxUlps is non-negative and small enough that the */ - /* default NAN won't compare as equal to anything. */ - assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024); - aInt = *(long64*)&A; - /* Make aInt lexicographically ordered as a twos-complement int */ - if (aInt < 0) -#ifdef _MSC_VER - aInt = 0x8000000000000000 - aInt; -#else - aInt = 0x8000000000000000LL - aInt; -#endif - bInt = *(long64*)&B; - /* Make bInt lexicographically ordered as a twos-complement int */ - if (bInt < 0) -#ifdef _MSC_VER - bInt = 0x8000000000000000 - bInt; -#else - bInt = 0x8000000000000000LL - bInt; -#endif -#ifdef _MSC_VER - intDiff = Abs(aInt - bInt); -#else - intDiff = llabs(aInt - bInt); -#endif -/* printf("A:%e B:%e aInt:%d bInt:%d diff:%d\n", A, B, aInt, bInt, intDiff); */ - if (intDiff <= maxUlps) - return TRUE; - return FALSE; -} - -#ifndef HAVE_LOGB - -double -logb(double x) -{ - double y = 0.0; - - if (x != 0.0) - { - if (x < 0.0) - x = - x; - while (x > 2.0) - { - y += 1.0; - x /= 2.0; - } - while (x < 1.0) - { - y -= 1.0; - x *= 2.0; - } - } - else - y = 0.0; - - return y; -} -#endif - - - - -#ifndef HAVE_SCALB -# ifdef HAVE_SCALBN -# define scalb scalbn -#else /* Chris Inbody */ - -double -scalb(double x, int n) -{ - double y, z = 1.0, k = 2.0; - - if (n < 0) { - n = -n; - k = 0.5; - } - - if (x != 0.0) - for (y = 1.0; n; n >>= 1) { - y *= k; - if (n & 1) - z *= y; - } - - return x * z; -} -# endif /* HAVE_SCALBN */ -#endif /* HAVE_SCALB */ - - -#ifndef HAVE_ERFC -/* From C. Hastings, Jr., Approximations for digital computers, - Princeton Univ. Press, 1955. - Approximation accurate to within 1.5E-7 - (making some assumptions about your machine's floating point mechanism) -*/ - -double -erfc(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); -} -#endif - - - -#ifndef HAVE_ISNAN -/* isnan (originally) for SOI devices in MINGW32 hvogt (dev.c) */ - -union ieee754_double -{ - double d; - - /* This is the IEEE 754 double-precision format. */ - struct - { - /* Together these comprise the mantissa. */ - unsigned int mantissa1:32; - unsigned int mantissa0:20; - unsigned int exponent:11; - unsigned int negative:1; - } ieee; - struct - { - /* Together these comprise the mantissa. */ - unsigned int mantissa1:32; - unsigned int mantissa0:19; - unsigned int quiet_nan:1; - unsigned int exponent:11; - unsigned int negative:1; - } ieee_nan; -}; - -int -isnan(double value) - -{ - union ieee754_double u; - - u.d = value; - - /* IEEE 754 NaN's have the maximum possible - exponent and a nonzero mantissa. */ - return ((u.ieee.exponent & 0x7ff) == 0x7ff && - (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)); - -} - -/* - * end isnan.c - */ -#endif /* HAVE_ISNAN */ - diff --git a/src/misc/missing_math.h b/src/misc/missing_math.h deleted file mode 100644 index 40358a72a..000000000 --- a/src/misc/missing_math.h +++ /dev/null @@ -1,29 +0,0 @@ -/************* - * Header file for missing_math.c - * 1999 E. Rouat - ************/ - -#ifndef MISSING_MATH_H_INCLUDED -#define MISSING_MATH_H_INCLUDED - -bool AlmostEqualUlps(double, double, int); - -#ifndef HAVE_ERFC -double erfc(double); -#endif - -#ifndef HAVE_LOGB -double logb(double); -#endif - -#ifndef HAVE_SCALB -# ifndef HAVE_SCALBN -double scalb(double, int); -# endif -#endif - -#ifndef HAVE_ISNAN -int isnan(double); -#endif - -#endif /* MISSING_MATH_H_INCLUDED */