From ff2138b1867e9bbc54d0d3b338eda08c798bcd14 Mon Sep 17 00:00:00 2001
From: rlar <rlar>
Date: Sun, 15 Jan 2012 17:19:35 +0000
Subject: [PATCH] hisim2, bug fix, a numerical problem in the hisim2 and
 hisimhv models

First seen when compiled with gcc-4.6.2 -g -O1

The macro `Fn_SZ' which boils down to
  1/2 * (x + sqrt(x*x + c*c))
was used in a context where a negative result
  blew up the following computations.
(used to compute `Egidl', which is used to
 compute exp(-1 / (Egidl + small_constant)))

For large negative values of x the computation
boils down to
  1/2 ( x + almost(|x|) )
where the summands almost cancel each other,
sometimes yielding a small negative result.

small_constant was too small to avoid a big
  result for -1/(Egidl + small_constant)
yielding an `inf' during the exp() computation,
which was later multiplied with a `0'
  yielding a NaN,
which was carried forward
  during the rest of computations.

Because the error of the cancellation scales
with the magnitude of x, no larger `small_constant'
could have avoided the problem.

Presumably the problem was amplified
  by a mixture of precisions (double versus extended float)
  of intermediate values.
(the program wasn't compiled for sse)

( x was -2.812500e+06,
  c was  1.000000e-02,
  Fn_SZ result was -1.853095e-11

  thus the cancellation remainder
  was of relative size
    6.6e-18
  which is approximately
    2^-57
  and thus more accurate
   as a `double float' could have delivered
)
---
 src/spicelib/devices/hisim2/hsm2eval.c | 2 ++
 src/spicelib/devices/hisim2/hsm2temp.c | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/src/spicelib/devices/hisim2/hsm2eval.c b/src/spicelib/devices/hisim2/hsm2eval.c
index e10fd5f66..00caef669 100644
--- a/src/spicelib/devices/hisim2/hsm2eval.c
+++ b/src/spicelib/devices/hisim2/hsm2eval.c
@@ -205,6 +205,8 @@ double TMF1 , TMF2 , TMF3 , TMF4 ;
     TMF2 = sqrt ( ( x ) *  ( x ) + 4.0 * ( delta ) * ( delta) ) ; \
     dx = 0.5 * ( 1.0 + ( x ) / TMF2 ) ; \
     y = 0.5 * ( ( x ) + TMF2 ) ; \
+    if (y < 0) \
+        y = 0 ; \
   }
 #if 0
 /*---------------------------------------------------*
diff --git a/src/spicelib/devices/hisim2/hsm2temp.c b/src/spicelib/devices/hisim2/hsm2temp.c
index 4d851c22f..b71071fd9 100644
--- a/src/spicelib/devices/hisim2/hsm2temp.c
+++ b/src/spicelib/devices/hisim2/hsm2temp.c
@@ -41,6 +41,8 @@
 #define Fn_SZtemp( y , x , delta ) { \
     T1 = sqrt ( ( x ) *  ( x ) + 4.0 * ( delta ) * ( delta) ) ; \
     y = 0.5 * ( ( x ) + T1 ) ; \
+    if (y < 0) \
+        y = 0 ; \
   }
 
 #define Fn_SUtemp( y , x , xmax , delta ) { \