diff --git a/ChangeLog b/ChangeLog index be8449d98..dcb051023 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,16 @@ +2008-01-02 Dietmar Warning + * src/frontend/outitf.c: Fixed rawfile ascii generation. Same like below. + * src/frontend/inp.c: don't need local buffer w/o getcwd + * src/conf.c: belong spice3 manual ascii is default anyway + * src/misc/missing_math.*, src/include/missig_math.h, /src/frontend/measure.c, + src/spicelib/analysis/dctran.c: Using a real double compare to equal. + 2008-01-02 Paolo Nenzi * src/frontend/rawfile.c: Fixed rawfile ascii generation. The prevoius patch produced incorrect string like v(v(1)) for v(1) in the output file. 2007-12-31 Holger Vogt - * src/frontend/com_chdir.c: fix for the crashing of ngspice under Windows when +don't need buffer w/o getcwd fix for the crashing of ngspice under Windows when started from windows explorer. * src/frontend/inp.c: ngspice crashed when executing a file consisting of a simple control section. Fixed. diff --git a/src/frontend/measure.c b/src/frontend/measure.c index 79f3dc238..190b79420 100644 --- a/src/frontend/measure.c +++ b/src/frontend/measure.c @@ -74,7 +74,7 @@ get_volt_time( struct dvec *time, struct dvec *values, double value, char polari *failed = TRUE; } } - if ( AlmostEqualUlps( comp_time, 0, 3 ) ) *failed = TRUE; + if ( AlmostEqualUlps( comp_time, 0, 100 ) ) *failed = TRUE; return comp_time; } @@ -139,7 +139,7 @@ measure2( char *meas_type, char *vec_name, char vec_type, double from, double to *result_time = time->v_realdata[i]; } else { *result = ( strcmp( meas_type, "max" ) == 0 ) ? max( *result, vec->v_realdata[i] ) : min( *result, vec->v_realdata[i] ); - if ( !AlmostEqualUlps( prev_result, *result, 3 ) ) *result_time = time->v_realdata[i]; + if ( !AlmostEqualUlps( prev_result, *result, 100 ) ) *result_time = time->v_realdata[i]; } } } @@ -164,23 +164,23 @@ measure2( char *meas_type, char *vec_name, char vec_type, double from, double to // see if y-value constant for ( i = 0; i < xy_size-1; i++ ) - if ( !AlmostEqualUlps( *(y+i), *(y+i+1), 3 ) ) constant_y = FALSE; + if ( !AlmostEqualUlps( *(y+i), *(y+i+1), 100 ) ) constant_y = FALSE; // Compute Integral (area under curve) i = 0; while ( i < xy_size-1 ) { // Simpson's 3/8 Rule - if ( AlmostEqualUlps( *(width+i), *(width+i+1), 3 ) && AlmostEqualUlps( *(width+i), *(width+i+2), 3 ) ) { + if ( AlmostEqualUlps( *(width+i), *(width+i+1), 100 ) && AlmostEqualUlps( *(width+i), *(width+i+2), 100 ) ) { sum1 += 3*(*(width+i))*(*(y+i) + 3*(*(y+i+1) + *(y+i+2)) + *(y+i+3))/8; i += 3; } // Simpson's 1/3 Rule - else if ( AlmostEqualUlps( *(width+i), *(width+i+1), 3 ) ) { + else if ( AlmostEqualUlps( *(width+i), *(width+i+1), 100 ) ) { sum2 += *(width+i)*(*(y+i) + 4*(*(y+i+1)) + *(y+i+2))/3; i += 2; } // Trapezoidal Rule - else if ( !AlmostEqualUlps( *(width+i), *(width+i+1), 3 ) ) { + else if ( !AlmostEqualUlps( *(width+i), *(width+i+1), 100 ) ) { sum3 += *(width+i)*(*(y+i) + *(y+i+1))/2; i++; } diff --git a/src/include/missing_math.h b/src/include/missing_math.h index fb7009ba8..52e72ceb9 100644 --- a/src/include/missing_math.h +++ b/src/include/missing_math.h @@ -8,7 +8,7 @@ Copyright 1999 Emmanuel Rouat #ifndef MISSING_MATH_H_INCLUDED #define MISSING_MATH_H_INCLUDED -bool AlmostEqualUlps(float, float, int); +bool AlmostEqualUlps(double, double, int); #ifndef HAVE_ERFC extern double erfc(double); diff --git a/src/misc/missing_math.c b/src/misc/missing_math.c index 421bd67b4..4ec08b51f 100644 --- a/src/misc/missing_math.c +++ b/src/misc/missing_math.c @@ -11,22 +11,59 @@ $Id$ #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)) -/* Initial AlmostEqualULPs version - fast and simple, but */ -/* some limitations. */ -bool AlmostEqualUlps(float A, float B, int maxUlps) +/* 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) { - int intDiff; - assert(sizeof(float) == sizeof(int)); + long64 aInt, bInt, intDiff; if (A == B) return TRUE; - intDiff = abs(*(int*)&A - *(int*)&B); - + /* 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/misc/missing_math.h b/src/misc/missing_math.h index 6b479d37b..40358a72a 100644 --- a/src/misc/missing_math.h +++ b/src/misc/missing_math.h @@ -6,7 +6,7 @@ #ifndef MISSING_MATH_H_INCLUDED #define MISSING_MATH_H_INCLUDED -bool AlmostEqualUlps(float, float, int); +bool AlmostEqualUlps(double, double, int); #ifndef HAVE_ERFC double erfc(double); diff --git a/src/spicelib/analysis/dctran.c b/src/spicelib/analysis/dctran.c index 6972f6676..bb4e32715 100644 --- a/src/spicelib/analysis/dctran.c +++ b/src/spicelib/analysis/dctran.c @@ -459,7 +459,7 @@ nextTime: ckt->CKTstat->STAToldIter = ckt->CKTstat->STATnumIter; if(check_autostop("tran") || fabs(ckt->CKTtime - ckt->CKTfinalTime) < ckt->CKTminBreak || - AlmostEqualUlps( ckt->CKTtime, ckt->CKTfinalTime, 3 ) ) { + AlmostEqualUlps( ckt->CKTtime, ckt->CKTfinalTime, 100 ) ) { #ifdef STEPDEBUG printf(" done: time is %g, final time is %g, and tol is %g\n", ckt->CKTtime,ckt->CKTfinalTime,ckt->CKTminBreak); @@ -526,7 +526,7 @@ resume: #ifdef XSPICE /* gtri - begin - wbk - Cut integration order if first timepoint after breakpoint */ //if(ckt->CKTtime == g_mif_info.breakpoint.last) - if ( AlmostEqualUlps( ckt->CKTtime, g_mif_info.breakpoint.last, 3 ) ) + if ( AlmostEqualUlps( ckt->CKTtime, g_mif_info.breakpoint.last, 100 ) ) ckt->CKTorder = 1; /* gtri - end - wbk - Cut integration order if first timepoint after breakpoint */ @@ -535,7 +535,7 @@ resume: /* are we at a breakpoint, or indistinguishably close? */ //if ((ckt->CKTtime == *(ckt->CKTbreaks)) || (*(ckt->CKTbreaks) - - if ( AlmostEqualUlps( ckt->CKTtime, *(ckt->CKTbreaks), 3 ) || (*(ckt->CKTbreaks) - + if ( AlmostEqualUlps( ckt->CKTtime, *(ckt->CKTbreaks), 100 ) || (*(ckt->CKTbreaks) - (ckt->CKTtime) <= ckt->CKTdelmin)) { /* first timepoint after a breakpoint - cut integration order */ /* and limit timestep to .1 times minimum of time to next breakpoint, @@ -589,7 +589,7 @@ resume: #ifdef STEPDEBUG printf(" brk_pt: %g ckt_time: %g ckt_min_break: %g\n",*(ckt->CKTbreaks), ckt->CKTtime, ckt->CKTminBreak); #endif - if(AlmostEqualUlps(*(ckt->CKTbreaks),ckt->CKTtime,3) || *(ckt->CKTbreaks) <= (ckt->CKTtime + ckt->CKTminBreak)) { + if(AlmostEqualUlps(*(ckt->CKTbreaks),ckt->CKTtime, 100) || *(ckt->CKTbreaks) <= (ckt->CKTtime + ckt->CKTminBreak)) { #ifdef STEPDEBUG printf("throwing out permanent breakpoint times <= current time (brk pt: %g)\n",*(ckt->CKTbreaks)); printf(" ckt_time: %g ckt_min_break: %g\n",ckt->CKTtime, ckt->CKTminBreak);