Re: [Bug-gnucap] insufficient limits checking in tanh model

[Al Davis]

On Saturday 05 April 2003 02:10 pm, address@hidden wrote:
> I've been looking at gnucap-0.33 a little and it seems that
> the bm_tanh.cc file has a potential problem on some systems. 
> In this function:
>
> void EVAL_BM_TANH::tr_eval(ELEMENT* d)const
> {
>   double x = ioffset(d->_y0.x);
>   double aa = x * _gain/_limit;
>   double cosine = cosh(aa);
>   double f1 = _gain / (cosine*cosine);
>   double f0 = _limit * tanh(aa);
>   d->_y0 = FPOLY1(x, f0, f1);
>   tr_final_adjust(&(d->_y0), d->f_is_value());
> }
>
> its possible for 'aa' to be large enough to overflow cosh and
> product 'inf' as the result.  Then you're relying on _gain /
> (inf * inf) to give 0.  In particular I've noticed that on an
> alpha system (alpha processors do not implement full IEEE 754
> math in hardware) I get floating point exceptions at that
> line of the code.  If I compile with the -mieee flag to gcc,
> traps are generated in the code to allow the software
> floating point completion code to kick in and the program
> runs without generating the SIGFPE's that I was getting
> before.  However, not all systems implement this sort of
> floating point math so you may find problems on some systems.

Can you post a circuit that demonstrates the problem?

It doesn't happen on any system I have access to, but I 
understand why it will on some others.

What does it do?

I expect it will print "floating point error" and continue, 
eventually giving you the correct answer.  Does the command 
"option noerror" suppress the message?  If that is what you 
get, and you can suppress it, that means the sigfpe trap is 
working correctly.

Aside from that ..

There should be a trap there.  Not all systems have sech, so I 
can't eliminate the problem that way.  Besides, it might be 
implemented as 1/cosh.  Also, tanh might be implemented as 
sinh/cosh, making it also prone to the problem.  All of these 
could be implemented using exp(x) and exp(-x), which is also 
prone to the problem.

I think the only reliable way to handle it is to trap large 
"aa". 

Try this:

const double LOGBIGBIG = log(BIGBIG);
/*--------------------------------------------------------------------------*/
void EVAL_BM_TANH::tr_eval(ELEMENT* d)const
{
  double x = ioffset(d->_y0.x);
  double aa = x * _gain/_limit;
  double f1, f0;
  {if (aa > LOGBIGBIG) {
    f1 = 0;
    f0 = _limit;
  }else if (aa < -LOGBIGBIG) {
    f1 = 0;
    f0 = -_limit;
  }else{
    double cosine = cosh(aa);
    f1 = _gain / (cosine*cosine);
    f0 = _limit * tanh(aa);
  }}
  d->_y0 = FPOLY1(x, f0, f1);
  tr_final_adjust(&(d->_y0), d->f_is_value());
}
/*--------------------------------------------------------------------------*/