Date: Mon, 1 Jul 1996 08:18:06 +0200 (IST)
From: Eli Zaretskii <eliz AT is DOT elta DOT co DOT il>
To: "Rafael R. Sevilla" <rsevilla AT upd DOT edu DOT ph>
Cc: DJGPP Mail Server <djgpp AT delorie DOT com>
Subject: Re: Floating point exceptions
In-Reply-To: <Pine.SOL.3.91.960625233137.905A-100000@sauron>
Message-Id: <Pine.SUN.3.91.960701075043.16855B-100000@is>
Mime-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII


On Sun, 30 Jun 1996, Rafael R. Sevilla wrote:

> There's a
> portion of code in the ray tracer that inverts a 3x3 matrix, and one of my
> scenes happens to create a matrix that has a zero determinant. The Sparc
> eats it, no sweat, and presses on without telling me anything! Since I use
> curses with it, I have a signal handler that takes care of such things,
> lest I mess up the screen, but I don't get a word from the Sparc. Can
> anyone tell me anything about this rather odd behavior? Must I explicitly 
> check for division by zero everywhere just to ensure that my code acts 
> the same way across all platforms?

IMHO, any code that causes division by zero as part of its normal
operation either has a bug or uses inappropriate algorithms.  If I
understand correctly, your code inverts matrices by computing their
determinant.  This is the worst possible method!  It is very unstable
(i.e. it magnifies the inaccuracies in the matrix elements by a large
factor), and I've seen it produce terribly wrong results even for
perfectly ``healthy'' matrices; with badly-conditioned matrices, it almost
always will divide by zero.  (It seems to be a rule that every method
which you've been taught in high school during your math classes turns out
to be numerically the most unstable.)  If you are serious about inverting
matrices (like if you want to make it part of a production code or a
public library), get any of the textbooks on matrix software and use any
of the algorithms described there.  And don't reinvent the wheel, just
copy code from one of the existing libraries floating around.  I recommend
an algorithm called SVD (Singular Value Decomposition).  This is *the*
most stable method ever, it will NEVER divide by zero or fail to produce a
meaningful result, even if your matrix is truely singular, although it 
will bloat your program by some 10KB (TANSTAAFL). 

Catching FP exceptions is a means for writing code that doesn't bomb even 
if there is a bug in it, or if a user feeds it with an input that doesn't 
make sense.  It should not be used (IMHO) for ``correcting'' results 
during normal operation with perfectly legal inputs.

> Is there something different about the way Sparcs and Pentiums handle
> floating point exceptions? I've been developing a ray tracer portable
> between djgpp, Solaris 5.3, and Linux, and the Pentium versions (djgpp and
> Linux) seem to not work the same way as the Sun version. In particular the
> Sparc seems to forget that it's dividing by zero somewhere!

I don't know about Linux, but for DJGPP you can try setting the flags of 
the FPU so that division by zero doesn't cause an exception.  Check out 
the description of `_control87' library function in the libc docs.

It is possible that on Sparc, the FPU is set this way by default.  In 
that case, you should silently get a NaN when you divide by zero.  Or 
maybe the OS catches the FP exception and does something about it 
(because otherwise the CPU cannot continue, like x86 sometimes does in 
case of integer zero-divide).