[fpc-pascal] Division by Zero: EDivByZero and EZeroDivide

[Bart]

Hi,

Is there a reason why in fpc both floating point and integer division
by zero raise an EDivByZero exception?

See: http://docwiki.embarcadero.com/VCL/en/SysUtils.EZeroDivide

SysUtils.EZeroDivide exception is raised when an application attempts
to divide a floating-point value by zero.
Note: Integer divide-by-zero errors raise the SysUtils.EDivByZero exception.

program zerodiv;

{$IFDEF FPC}
{$mode objfpc}{$H+}
{$ENDIF}

uses Sysutils;

var r,x: extended;
i,j: integer;

begin
i := 0;
r := 0.0;
try
write('Integer division by zero: 1 div 0 -> ');
j := 1 div i;
writeln(j);
except
on e: exception do writeln(e.classname);
end;
try
write('Float division by zero : 1.0 / 0.0 -> ');
x := 1.0 / r;
writeln(x);
except
on e: exception do writeln(e.classname);
end;
end.

Output with Delphi (3.0:
Integer division by zero: 1 div 0 -> EDivByZero
Float division by zero : 1.0 / 0.0 -> EZeroDivide

Output with Fpc 2.2.4 (tested on win32 and Linux i386)
Integer division by zero: 1 div 0 -> EDivByZero
Float division by zero : 1.0 / 0.0 -> EDivByZero

Is this by design, or should I report it as a bug?
(B.t.w. setting mode to delphi makes no difference)

Bart

[Tom Verhoeff]

Post by Bart
Is there a reason why in fpc both floating point and integer division
by zero raise an EDivByZero exception?
See: http://docwiki.embarcadero.com/VCL/en/SysUtils.EZeroDivide
SysUtils.EZeroDivide exception is raised when an application attempts
to divide a floating-point value by zero.
Note: Integer divide-by-zero errors raise the SysUtils.EDivByZero exception.

It would be nicer if one had the ability to make floating-point division
by zero return an IEEE 754 plus/minus infinity, without raising an
exception. This is e.g. useful when one needs to evaluate rational
functions. Without such an infinity, you need to make a nasty case
analysis, which furthermore depends on the rational function.

A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk. The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0. When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities. Also see

<http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf>

for other examples and further motivation.

Unfortunately, many compilers still do not fully support the IEEE 754
floating-point standard.

Best regards,

Tom

[Jorge Aldo G. de F. Junior]

Wouldnt a NaN (Not a number) be more "matematically correct" result (I saw
that on an old book about i387)

Matematically division by zero is an "mathematical impossibility", so NaN
would be more logical

Is there a way to deal with NaN's and other i387 (and similars) conventions
directly in Pascal ?

Something like :

If A / B = NaN then something();

Post by Bart

Post by Bart
Is there a reason why in fpc both floating point and integer division
by zero raise an EDivByZero exception?
See: http://docwiki.embarcadero.com/VCL/en/SysUtils.EZeroDivide
SysUtils.EZeroDivide exception is raised when an application attempts
to divide a floating-point value by zero.
Note: Integer divide-by-zero errors raise the SysUtils.EDivByZero

exception.
It would be nicer if one had the ability to make floating-point division
by zero return an IEEE 754 plus/minus infinity, without raising an
exception. This is e.g. useful when one needs to evaluate rational
functions. Without such an infinity, you need to make a nasty case
analysis, which furthermore depends on the rational function.
A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk. The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0. When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities. Also see
<http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf<http://www.cs.berkeley.edu/%7Ewkahan/ieee754status/why-ieee.pdf>
for other examples and further motivation.
Unfortunately, many compilers still do not fully support the IEEE 754
floating-point standard.
Best regards,
Tom
--
PHONE: [phone] | Technische Universiteit Eindhoven
FAX: [phone] | PO Box 513, NL-5600 MB Eindhoven
http://www.win.tue.nl/~wstomv/ <http://www.win.tue.nl/%7Ewstomv/> | The
Netherlands
_______________________________________________
http://lists.freepascal.org/mailman/listinfo/fpc-pascal

[Jürgen Hestermann]

Post by Jorge Aldo G. de F. Junior
Is there a way to deal with NaN's and other i387 (and similars)
conventions directly in Pascal ?

Yes, if not yet possible (I didn't have a closer look) I would like to
have that too. But I think it's heavily depending on the processor...

(Details about NaN's and infinity at http://en.wikipedia.org/wiki/NaN)

[Aleksa Todorovic]

On Sun, Oct 18, 2009 at 18:49, Jorge Aldo G. de F. Junior

Post by Jorge Aldo G. de F. Junior
Wouldnt a NaN (Not a number) be more "matematically correct" result (I saw
that on an old book about i387)
Matematically division by zero is an "mathematical impossibility", so NaN
would be more logical
Is there a way to deal with NaN's and other i387 (and similars) conventions
directly in Pascal ?

Take a look at SysUtils, functions IsNaN, IsInfinite.

Post by Jorge Aldo G. de F. Junior

Post by Tom Verhoeff

Post by Bart
Is there a reason why in fpc both floating point and integer division
by zero raise an EDivByZero exception?
See: http://docwiki.embarcadero.com/VCL/en/SysUtils.EZeroDivide
SysUtils.EZeroDivide exception is raised when an application attempts
to divide a floating-point value by zero.
Note:  Integer divide-by-zero errors raise the SysUtils.EDivByZero
exception.

Maybe because is more consistent to have one type for both exceptions
- when you see it, you know what it means. How could anybody make
distinction between EDivByZero and EZeroDivide without looking at
documentation?

Post by Jorge Aldo G. de F. Junior

Post by Tom Verhoeff
It would be nicer if one had the ability to make floating-point division
by zero return an IEEE 754 plus/minus infinity, without raising an
exception.  This is e.g. useful when one needs to evaluate rational
functions.  Without such an infinity, you need to make a nasty case
analysis, which furthermore depends on the rational function.

There is function SetExceptionMask in SysUtils which you can use to
control which exceptions will be thrown at runtime. I haven't used it,
but looks like a good place to start...

[Tom Verhoeff]

Post by Jorge Aldo G. de F. Junior
Wouldnt a NaN (Not a number) be more "matematically correct" result (I saw
that on an old book about i387)

No, read Kahan's article that I pointed to.

Post by Jorge Aldo G. de F. Junior
Matematically division by zero is an "mathematical impossibility", so NaN
would be more logical

There are situations (that occur e.g. in numerical approximations of
mathematical models in physics) where division by zero is a well-defined
infinity. For that reason, the IEEE 754 floating-point standard has
distinct plus zero, minus zero, plus infinity, and minus infinity.
There are precise (and useful) rules to calculate with these special
numbers. They are not NaNs.

Many processors fully support the IEEE 754 standard, but many compilers
don't. Partial support of the standard in the compiler makes the standard
much less useful.

In our department, especially the numericists disliked (the old)
Pascal because of the fixed array sizes and the lack of good support
for floating-point arithmetic (i.e. IEEE 754 since 1985). I guess few
(Free)Pascal users nowadays come from the numerical/scientific computing
community.

Another (longish) "must read" about floating-point:

David Goldberg.
"What Every Computer Scientist Should Know about Floating-Point Arithmetic".
<http://www.validlab.com/goldberg/paper.ps>

Note that IEEE Std 754 was revised in 2008, see e.g.

<http://en.wikipedia.org/wiki/IEEE_754-2008>

Best regards,

Tom

[Jonas Maebe]

Post by Tom Verhoeff
It would be nicer if one had the ability to make floating-point division
by zero return an IEEE 754 plus/minus infinity, without raising an
exception.

http://www.freepascal.org/docs-html/rtl/math/setexceptionmask.html

Fully cross-platform, even.


Jonas

[Bart]

Hi all,

Interesting thoughts.

However, if we raise an exception, why do we raise EDivByZero instead
of EZeroDivide?

Note: EDivByZero = class(EIntError), so IMHO it makes no sense to
raise it with a floating point error.

Bart

[Tom Verhoeff]

Post by Jonas Maebe

Post by Tom Verhoeff
It would be nicer if one had the ability to make floating-point division
by zero return an IEEE 754 plus/minus infinity, without raising an
exception.

http://www.freepascal.org/docs-html/rtl/math/setexceptionmask.html
Fully cross-platform, even.

How new is that?

It works 'partly'. Consider the program

program DivideByZero;
{ to see whether division by zero can result in Infinite values }

uses Math;

var one, zero, result: Double;

begin
SetExceptionMask ( [ exZeroDivide{, exPrecision} ] );
one := 1.0;
zero := 0.0;
result := one / zero;
writeln ( result ); { Expected output: +Inf }
result := one / result;
writeln ( result ); { Expected output: +0.0 }
end.

Outputs:
+Inf
An unhandled exception occurred at $000060B0 :
EDivByZero : Division by zero
$000060B0
$00005850
$0001889B

In order to allow 1 / +Inf, you also need to mask exPrecision
(i.e. uncomment it from the set passed to SetExceptionMask).
I am surprised that division by an Infinity raises EDivByZero: Division by zero

Best regards,

Tom

[Jonas Maebe]

Post by Tom Verhoeff

Post by Jonas Maebe
http://www.freepascal.org/docs-html/rtl/math/setexceptionmask.html
Fully cross-platform, even.

How new is that?

Not very. I think it exists at least since 2.0.0.

Post by Tom Verhoeff
It works 'partly'. Consider the program
program DivideByZero;
{ to see whether division by zero can result in Infinite values }
uses Math;
var one, zero, result: Double;
begin
SetExceptionMask ( [ exZeroDivide{, exPrecision} ] );
one := 1.0;
zero := 0.0;
result := one / zero;
writeln ( result ); { Expected output: +Inf }
result := one / result;
writeln ( result ); { Expected output: +0.0 }
end.
+Inf
EDivByZero : Division by zero
$000060B0
$00005850
$0001889B
In order to allow 1 / +Inf, you also need to mask exPrecision
(i.e. uncomment it from the set passed to SetExceptionMask).
Division by zero

It probably depends on your platform. In case it's Mac OS X on x86:
floating point exception reporting via Unix signals isn't exactly its
forte. The Darwin signal handler contains code to decode the
instructions where a floating point exception occurs, because the
system generates a "generic" floating point exception on integer div-
by-zero. Maybe the above operation also causes a generic floating
point exception which is then wrongly classified by us as div-by-zero.

In general, relying on these exceptions to return a particular value
is a bad idea though. Using Unix signals for this purpose is not
supported at all on Mac OS X (we should actually use Mach exceptions
instead), because the state you get in the Unix signal handler is
transient (if another exception happens at the same time in another
thread, you may lose one exception, you may get a state that contains
a mix of both exceptions, etc). Moreover, both the format of the
context parameter of Unix signals and the interface for Mach
exceptions are still considered SPI (system private interface) by
Apple, so at this point there is no way to implement either in a way
that's guaranteed to be forward compatible.


Jonas

[Jonas Maebe]

In case it's Mac OS X on x86: floating point exception reporting via
Unix signals isn't exactly its forte.

And in case Mac OS X on ppc: the same caveats apply as for x86, except
that we don't have to any opcode decoding since the PPC doesn't raise
a hardware exception for integer div-by-zero. The system does appear
to raise the wrong exception for your test program though, but there's
not much we can do about that. This is basically our signal handler:

case sig of
SIGFPE :
begin
Case Info^.si_code Of
FPE_FLTDIV,
FPE_INTDIV : Res:=200; { floating point divide by zero }
FPE_FLTOVF : Res:=205; { floating point overflow }
FPE_FLTUND : Res:=206; { floating point underflow }
FPE_FLTRES, { floating point inexact result }
FPE_FLTINV : Res:=207; { invalid floating point operation }
Else
Res:=207; {coprocessor error}
end;

If the system say its FPE_FLTDIV (which according to /usr/include/sys/
signal.h means "[XSI] floating point divide by zero"), we report that...


Jonas

[Bart]

Jonas,

Could you comment on my opinion that in case of a floating point
divide by zero error an EZeroDivide should be raised (or no exception
at all depending on the FPU exception mask) and only in case of
integer divide by zero (div operation) an EDivByZero should be raised?

This was my original question and no one seems to have an opinion on that.

Bart

[Jonas Maebe]

Post by Bart
Could you comment on my opinion that in case of a floating point
divide by zero error an EZeroDivide should be raised (or no exception
at all depending on the FPU exception mask) and only in case of
integer divide by zero (div operation) an EDivByZero should be raised?
This was my original question and no one seems to have an opinion on that.

I don't have an opinion on that either. I'm not sure how it would be
technically possible with our current setup, as SysUtils intercepts
run time errors and then converts these into exceptions based on their
number. Since all divide-by-zero run time errors are 200, there's no
way for sysutils to differentiate between the two at this time (which
is probably the reason why it always raises EDivByZero).


Jonas

[Bart]

Jonas,

I might be just babbling, since I really have no idea how it's done, but:

In the code that "translates" the "you did something nasty" to
runerror(some_exit_code) it may be possible to differentiate between
float and integer fault?
This info could then be stored in some variable/object that can be
used by the code that intercepts the run-time error and raises the
actual fpc exception?

Bart

[Jonas Maebe]

Post by Bart
In the code that "translates" the "you did something nasty" to
runerror(some_exit_code) it may be possible to differentiate between
float and integer fault?

It is (on some platforms with some hacking).

Post by Bart
This info could then be stored in some variable/object that can be
used by the code that intercepts the run-time error and raises the
actual fpc exception?

In theory this is of course possible, but in practice this is
extremely error prone due to the nature of signal handlers.
Furthermore, adding extra global variables to the interface of the
system unit is very much frowned upon, and this would be required
since the signal handling happens in the system unit and the
translation to exceptions in the sysutils unit. Another, and much
easier, option would be to change the run time error number for either
integer or floating point div-by-zero. But that would obviously cause
its own share of compatibility problems.


Jonas

[Bart]

Jonas,

In theory this is of course possible, but in practice this is extremely
error prone due to the nature of signal handlers. Furthermore, adding extra
global variables to the interface of the system unit is very much frowned
upon, and this would be required since the signal handling happens in the
system unit and the translation to exceptions in the sysutils unit.

I think this is how it's done in Delphi though.

Another,
and much easier, option would be to change the run time error number for
either integer or floating point div-by-zero. But that would obviously cause
its own share of compatibility problems.

Yes, a different runtime error is also not Delphi compatible ;-(
(But probably not likely te be relied on in existing code?)

One other thought. In the code that translates the signal to a
runerror, could we pass the signal info as an extra, optional
parameter to the runerror procedure, or would that also be "invisible"
to the sysutils unit, when intercepting runerrors?

If we cannot reliably differentiate in the signal handler between
integer divide by zero and floating point divide by zero, then we'd
better not try to fix it. Better be reliably inconsistent with delphi
than not being able to rely on what exception is raised.
In that case we should probably document this behaviour somewhere?

Bart

[Gustavo Enrique Jimenez]

Post by Tom Verhoeff
A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk.  The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0.  When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities.  Also see
       <http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf>
for other examples and further motivation.

Mmm... the formula R = 1/(1/R1 + 1/R2 + ... + 1/Rk) is only valid if
none of Rn=0.
Programmers can and must take care of that situation.

Gustavo

[Frank Peelo]

Post by Gustavo Enrique Jimenez

Post by Tom Verhoeff
A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk. The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0. When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities. Also see
<http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf>
for other examples and further motivation.

Mmm... the formula R = 1/(1/R1 + 1/R2 + ... + 1/Rk) is only valid if
none of Rn=0.
Programmers can and must take care of that situation.

Agree: treating 1/0 as you would treat a finite number makes me feel
uncomfortable. See for instance
http://www.cocoa.uk.com/?p=63

Frank

[Tom Verhoeff]

Post by Frank Peelo

Post by Gustavo Enrique Jimenez

Post by Tom Verhoeff
A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk. The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0. When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities. Also see
<http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf>
for other examples and further motivation.

Mmm... the formula R = 1/(1/R1 + 1/R2 + ... + 1/Rk) is only valid if
none of Rn=0.
Programmers can and must take care of that situation.

Agree: treating 1/0 as you would treat a finite number makes me feel
uncomfortable. See for instance
http://www.cocoa.uk.com/?p=63

One final comment here. The field of scientific computing, including
numerical methods, is a specialism by itself, and there are many
misunderstandings about it. Respected experts from this field (including
William Kahan) have developed a floating-point standard and have been
able to convince hardware builders to implement this standard in their
(co)processors. Compiler writers would do well to trust those experts,
especially if they do not understand them.

So, yes 1/0 is not a finite number, but a (signed) infinity: 1/+0 = +Inf,
1/-0 = -Inf. Yes, there are two 0s among the floating-point numbers,
and they are all well-behaved, and some operators work (1 + +Inf = +Inf),
and others do not (+Inf + -Inf is undefined). With good reasons.

Not everybody needs to be a floating-point expert, but compiler users
that do know about this, would want to use the IEEE 754 standard, because
it makes life so much simpler for them. That may be the reason why
FORTRAN is still used in scientific computing (robust numerical libraries
and good support in the compiler; but otherwise not a language I can like).

Tom

[Tom Verhoeff]

Post by Gustavo Enrique Jimenez

Post by Tom Verhoeff
A simple example is the situation where one needs to calculate
the replacement resistor value R for parallel resistors having
values R1, R2, ..., Rk.  The formula is R = 1/(1/R1 + 1/R2 + ... + 1/Rk).
The formula gives a divide-by-zero if one of the resistors has value 0.
But in that case, the replacement value R also equals 0.  When allowing
infinities, it just works out fine (infinity + x = infinity, 1/infty = 0).
That is precisely why IEEE 754 has infinities.  Also see
       <http://www.cs.berkeley.edu/~wkahan/ieee754status/why-ieee.pdf>
for other examples and further motivation.

Mmm... the formula R = 1/(1/R1 + 1/R2 + ... + 1/Rk) is only valid if
none of Rn=0.
Programmers can and must take care of that situation.

No, that is just the point. The "easy" formula is actually not not
correct for zero resistors when adopting classical evaluation rules.
However, it works out fine with infinities. Take the case of two resistors

1/(1/R1 + 1/R2) = R1 R2 / (R1 + R2)

The left-hand side is continuous at R1=0 and at R2=0, but strictly
speaking not defined in those points. The formula on the right is the
same for R1, R2 not zero and it does not have a problem at R1=0 or R2=0
(unless R1=R2=0). For more than two resistors, rewriting as in the
right-hand side is (practically) out of the question, but the same holds.

Point is that for such rational functions, the use of infinities is
well-defined and works well. Numericists want to evaluate the
left-hand side, rather than doing a case analysis or rewriting it.

Tom