Using a real double compare to equal.

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 <p.nenzi@ieee.org>
  * 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.

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

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

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

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

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