From: Kbwms AT aol DOT com
Message-ID: <1d9.f96c51c.2c767308@aol.com>
Date: Thu, 21 Aug 2003 15:10:00 EDT
Subject: Arithmetic Exceptions in C99
To: djgpp-workers AT delorie DOT com
MIME-Version: 1.0
Content-Type: multipart/mixed; boundary="part1_1d9.f96c51c.2c767308_boundary"
X-Mailer: 8.0 for Windows sub 6011
Reply-To: djgpp-workers AT delorie DOT com
--part1_1d9.f96c51c.2c767308_boundary
Content-Type: multipart/alternative;
boundary="part1_1d9.f96c51c.2c767308_alt_boundary"
--part1_1d9.f96c51c.2c767308_alt_boundary
Content-Type: text/plain; charset="US-ASCII"
Content-Transfer-Encoding: 7bit
In a recent communication with Eli Zaretskii, the problem of raising
arithmetic exceptions, as described in C99, arose. His latest thoughts are appended
to this email as a postscript. Paragraph F.9 of C99 is attached to this email
for those who do not have C99 on their system.
In C99, functions nearbyint are specifically prohibited from endangering the
setting of the inexact exception. I am not aware of any other functions where
a similar restriction exists. So, what do we do about exceptions? Except
for nearbyintl(), my functions are letting it all hang out -- whatever happens
rules. In some cases, such as acoshl() when the argument is less than 1, errno
is set to EDOM, the return value is set to NaN, and the invalid exception
flag is raised. It seems to me like that is enough noise.
What do you think?
For those who have been pining for a copy of C99, here is a link:
http://anubis.dkuug.dk/JTC1/SC22/WG14/www/docs/n869/
KB Williams
PS
In a message dated 8/21/2003 9:57:50 AM Eastern Standard Time,
eliz AT elta DOT co DOT il writes:
> >From: Kbwms AT aol DOT com
> >Date: Wed, 20 Aug 2003 14:57:01 EDT
> >
> >The relevant fragments are spread all over the place. Fortunately,
> >paragraph F9 of C99 has everything regarding exceptions in one place
> >-- along with a lot of other stuff. It is for installations that
> >are, or intend to be, IEC 60559 compliant.
> >
> >Is this the kind of stuff that djgpp-workers will look at?
>
> Yes.
>
> >Perhaps, it would be better if it were posted somewhere. How to I do that?
>
> I think you could simply mention F9, as almost everybody on the list
> have their copy of C99 and could read the text if they don't remember.
>
> The main issue, as far as I'm concerned, is this: when C99 says
> ``raises such-and-such exception'', does that mean that we need to
> actually trigger a SIGFPE, or merely that the appropriate bit in
> fexcept_t should be set?
>
--part1_1d9.f96c51c.2c767308_alt_boundary
Content-Type: text/html; charset="US-ASCII"
Content-Transfer-Encoding: quoted-printable
In a recent communication with Eli Zaretskii, the proble=
m of raising arithmetic exceptions, as described in C99, arose. His la=
test thoughts are appended to this email as a postscript. Paragraph F.=
9 of C99 is attached to this email for those who do not have C99 on their sy=
stem.
In C99, functions nearbyint are specifically prohibited from endangering the=
setting of the inexact exception. I am not aware of any other functio=
ns where a similar restriction exists. So, what do we do about excepti=
ons? Except for nearbyintl(), my functions are letting it all hang out=
-- whatever happens rules. In some cases, such as acoshl() when the a=
rgument is less than 1, errno is set to EDOM, the return value is set to NaN=
, and the invalid exception flag is raised. It seems to me like that i=
s enough noise.
What do you think?
For those who have been pining for a copy of C99, here is a link:
http://anu=
bis.dkuug.dk/JTC1/SC22/WG14/www/docs/n869/
KB Williams
PS
In a message dated 8/21/2003 9:57:50 AM Eastern Standard Time, eliz AT elta DOT co.=
il writes:
>From: Kbwms AT aol DOT com
>Date: Wed, 20 Aug 2003 14:57:01 EDT
>
>The relevant fragments are spread all over the place. Fortunately,=
>paragraph F9 of C99 has everything regarding exceptions in one place
>-- along with a lot of other stuff. It is for installations that
>are, or intend to be, IEC 60559 compliant.
>
>Is this the kind of stuff that djgpp-workers will look at?
Yes.
>Perhaps, it would be better if it were posted somewhere. How to I do tha=
t?
I think you could simply mention F9, as almost everybody on the list
have their copy of C99 and could read the text if they don't remember.
The main issue, as far as I'm concerned, is this: when C99 says
``raises such-and-such exception'', does that mean that we need to
actually trigger a SIGFPE, or merely that the appropriate bit in
fexcept_t should be set?
--part1_1d9.f96c51c.2c767308_alt_boundary--
--part1_1d9.f96c51c.2c767308_boundary
Content-Type: text/plain; name="C99-F9.txt"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: inline; filename="C99-F9.txt"
505 Committee Draft -- January 18, 1999 WG14/N869
F.9 Mathematics
[#1] This subclause contains specifications of
facilities that are particularly suited for IEC 60559
implementations.
[#2] The Standard C macro HUGE_VAL and its float and long
double analogs, HUGE_VALF and HUGE_VALL, expand to
expressions whose values are positive infinities.
[#3] Special cases for functions in are covered
directly or indirectly by IEC 60559. The functions that IEC
60559 specifies directly are identified in F.3. The other
F.9 IEC 60559 floating-point arithmetic F.9
=0C
506 Committee Draft -- January 18, 1999 WG14/N869
functions in treat infinities, NaNs, signed zeros,
subnormals, and (provided the state of the FENV_ACCESS
pragma is on) the exception flags in a manner consistent
with the basic arithmetic operations covered by IEC 60559.
[#4] The invalid and divide-by-zero exceptions are raised as
specified in subsequent subclauses of this annex.
[#5] The overflow exception is raised whenever an infinity
-- or, because of rounding direction, a maximal-magnitude
finite number -- is returned in lieu of a value whose
magnitude is too large.
[#6] The underflow exception is raised whenever a result is
tiny (essentially subnormal or zero) and suffers loss of
accuracy.298)
[#7] Whether or when the trigonometric, hyperbolic, base-e
exponential, base-e logarithmic, error, and log gamma
functions raise the inexact exception is implementation-
defined. For other functions, the inexact exception is
raised whenever the rounded result is not identical to the
mathematical result.
[#8] Whether the inexact exception may be raised when the
rounded result actually does equal the mathematical result
is implementation-defined. Whether the underflow (and
inexact) exception may be raised when a result is tiny but
not inexact is implementation-defined.299) Otherwise, as
implied by F.7.6, the functions do not raise
spurious exceptions (detectable by the user).
[#9] Whether the functions honor the rounding direction mode
is implementation-defined.
[#10] Functions with a NaN argument return a NaN result and
raise no exception, except where stated otherwise.
[#11] The specifications in the following subclauses append
to the definitions in . For families of functions,
the specifications apply to all of the functions even though
only the principal function is shown.
Recommended practice
____________________
298IEC 60559 allows different definitions of underflow.
They all result in the same values, but differ on when
the exception is raised.
299It is intended that undeserved underflow and inexact
exceptions are raised only if determining inexactness
would be too costly.
F.9 IEC 60559 floating-point arithmetic F.9
=0C
WG14/N869 Committee Draft -- January 18, 1999 507
[#12] If a function with one or more NaN arguments returns a
NaN result, the result should be the same as one of the NaN
arguments (after possible type conversion), except perhaps |
for the sign.
F.9.1 Trigonometric functions
F.9.1.1 The acos functions
[#1]
-- acos(1) returns +0.
-- acos(x) returns a NaN and raises the invalid exception
for |x|>1.
F.9.1.2 The asin functions
[#1]
-- asin(=B10) returns =B10.
-- asin(x) returns a NaN and raises the invalid exception
for |x|>1.
F.9.1.3 The atan functions
[#1]
-- atan(=B10) returns =B10.
-- atan(=B1) returns =B1pi/2.
F.9.1.4 The atan2 functions
[#1]
-- atan2(=B10, x) returns =B10, for x>0.
-- atan2(=B10, +0) returns =B10.300)
-- atan2(=B10, x) returns =B1pi, for x<0.
-- atan2(=B10, -0) returns =B1pi.
-- atan2(y, =B10) returns pi/2 for y>0.
-- atan2(y, =B10) returns -pi/2 for y<0.
____________________
300atan2(0, 0) does not raise the invalid exception, nor
does atan2(y, 0) raise the divide-by-zero exception.
F.9 IEC 60559 floating-point arithmetic F.9.1.4
=0C
508 Committee Draft -- January 18, 1999 WG14/N869
-- atan2(=B1y, ) returns =B10, for finite y>0.
-- atan2(=B1, x) returns =B1pi/2, for finite x.
-- atan2(=B1y, -) returns =B1pi, for finite y>0.
-- atan2(=B1, ) returns =B1pi/4.
-- atan2(=B1, -) returns =B13pi/4.
F.9.1.5 The cos functions
[#1]
-- cos(=B10) returns 1.
-- cos(=B1) returns a NaN and raises the invalid exception.
F.9.1.6 The sin functions
[#1]
-- sin(=B10) returns =B10.
-- sin(=B1) returns a NaN and raises the invalid exception.
F.9.1.7 The tan functions
[#1]
-- tan(=B10) returns =B10.
-- tan(=B1) returns a NaN and raises the invalid exception.
F.9.2 Hyperbolic functions
F.9.2.1 The acosh functions
[#1]
-- acosh(1) returns +0.
-- acosh(+) returns +.
-- acosh(x) returns a NaN and raises the invalid exception
if x<1.
F.9.1.4 IEC 60559 floating-point arithmetic F.9.2.1
=0C
WG14/N869 Committee Draft -- January 18, 1999 509
F.9.2.2 The asinh functions
[#1]
-- asinh(=B10) returns =B10.
-- asinh(=B1) returns =B1.
F.9.2.3 The atanh functions
[#1]
-- atanh(=B10) returns =B10.
-- atanh(=B11) returns =B1 and raises the divide-by-zero
exception.
-- atanh(x) returns a NaN and raises the invalid exception
if |x|>1.
F.9.2.4 The cosh functions
[#1]
-- cosh(=B10) returns 1.
-- cosh(=B1) returns +.
F.9.2.5 The sinh functions
[#1]
-- sinh(=B10) returns =B10.
-- sinh(=B1) returns =B1.
F.9.2.6 The tanh functions
[#1]
-- tanh(=B10) returns =B10.
-- tanh(=B1) returns =B11.
F.9.3 Exponential and logarithmic functions
F.9.2.1 IEC 60559 floating-point arithmetic F.9.3
=0C
510 Committee Draft -- January 18, 1999 WG14/N869
F.9.3.1 The exp functions
[#1]
-- exp(=B10) returns 1.
-- exp(+) returns +.
-- exp(-) returns +0.
F.9.3.2 The exp2 functions
[#1]
-- exp2(=B10) returns 1.
-- exp2(+) returns +.
-- exp2(-) returns +0.
F.9.3.3 The expm1 functions
[#1]
-- expm1(=B10) returns =B10.
-- expm1(+) returns +.
-- expm1(-) returns -1.
F.9.3.4 The frexp functions
[#1]
-- frexp(=B10, exp) returns =B10, and stores 0 in the object
pointed to by exp.
-- frexp(=B1, exp) returns =B1, and stores an unspecified
value in the object pointed to by exp.
-- frexp(x, exp) stores an unspecified value in the object
pointed to by exp (and returns a NaN) when x is a NaN.
-- frexp raises no exception.
[#2] On a binary system, the body of the frexp function |
might be
{ |
*exp =3D (value =3D=3D 0) ? 0 : (int)(1 + logb(value)=
);|
return scalbn(value, -(*exp)); |
} |
F.9.3 IEC 60559 floating-point arithmetic F.9.3.4
=0C
WG14/N869 Committee Draft -- January 18, 1999 511
F.9.3.5 The ilogb functions
[#1] No additional requirements.
F.9.3.6 The ldexp functions
[#1] On a binary system, ldexp(x, exp) is equivalent to
scalbn(x, exp)
F.9.3.7 The log functions
[#1]
-- log(=B10) returns - and raises the divide-by-zero
exception.
-- log(1) returns +0.
-- log(x) returns a NaN and raises the invalid exception
if x<0.
-- log(+) returns +.
F.9.3.8 The log10 functions
[#1]
-- log10(=B10) returns - and raises the divide-by-zero
exception.
-- log10(1) returns +0.
-- log10(x) returns a NaN and raises the invalid exception
if x < 0.
-- log10(+) returns +.
F.9.3.9 The log1p functions
[#1]
-- log1p(=B10) returns =B10.
-- log1p(-1) returns - and raises the divide-by-zero
exception.
-- log1p(x) returns a NaN and raises the invalid exception
if x<-1.
-- log1p(+) returns +.
F.9.3.4 IEC 60559 floating-point arithmetic F.9.3.9
=0C
512 Committee Draft -- January 18, 1999 WG14/N869
F.9.3.10 The log2 functions
[#1]
-- log2(=B10) returns - and raises the divide-by-zero
exception.
-- log2(x) returns a NaN and raises the invalid exception
if x<0.
-- log2(+) returns +.
F.9.3.11 The logb functions
[#1]
-- logb(=B1) returns +.
-- logb(=B10) returns - and raises the divide-by-zero
exception.
F.9.3.12 The modf functions
[#1]
-- modf(value, iptr) returns a result with the same sign
as the argument value.
-- modf(=B1, iptr) returns =B10 and stores =B1 in the object
pointed to by iptr.
-- modf of a NaN argument stores a NaN in the object
pointed to by iptr (and returns a NaN).
[#2] modf behaves as though implemented by
#include
#include
#pragma STDC FENV_ACCESS ON
double modf(double value, double *iptr)
{
int save_round =3D fegetround();
fesetround(FE_TOWARDZERO);
*iptr =3D nearbyint(value);
fesetround(save_round);
return copysign(
isinf(value) ? 0.0 :
value - (*iptr), value);
}
F.9.3.9 IEC 60559 floating-point arithmetic F.9.3.12
=0C
WG14/N869 Committee Draft -- January 18, 1999 513
F.9.3.13 The scalbn and scalbln functions
[#1]
-- scalbn(x, n) returns x if x is infinite or zero.
-- scalbn(x, 0) returns x.
F.9.4 Power and absolute value functions
F.9.4.1 The cbrt functions
[#1]
-- cbrt(=B1) returns =B1.
-- cbrt(=B10) returns =B10.
F.9.4.2 The fabs functions
[#1]
-- fabs(=B10) returns +0.
-- fabs(=B1) returns +.
F.9.4.3 The hypot functions
[#1]
-- hypot(x, y), hypot(y, x), and hypot(x, -y) are
equivalent.
-- hypot(x, y) returns + if x is infinite, even if y is a
NaN.
-- hypot(x, =B10) is equivalent to fabs(x).
F.9.4.4 The pow functions
[#1]
-- pow(x, =B10) returns 1 for any x, even a NaN.
-- pow(x, +) returns + for |x|>1.
-- pow(x, +) returns +0 for |x|<1.
-- pow(x, -) returns +0 for |x|>1.
-- pow(x, -) returns + for |x|<1.
-- pow(+, y) returns + for y>0.
F.9.3.13 IEC 60559 floating-point arithmetic F.9.4.4
=0C
514 Committee Draft -- January 18, 1999 WG14/N869
-- pow(+, y) returns +0 for y<0.
-- pow(-, y) returns - for y an odd integer > 0.
-- pow(-, y) returns + for y>0 and not an odd integer.
-- pow(-, y) returns -0 for y an odd integer < 0.
-- pow(-, y) returns +0 for y<0 and not an odd integer.
-- pow(=B11, =B1) returns a NaN and raises the invalid
exception.
-- pow(x, y) returns a NaN and raises the invalid
exception for finite x<0 and finite non-integer y.
-- pow(=B10, y) returns =B1 and raises the divide-by-zero
exception for y an odd integer < 0.
-- pow(=B10, y) returns + and raises the divide-by-zero
exception for y<0 and not an odd integer.
-- pow(=B10, y) returns =B10 for y an odd integer > 0.
-- pow(=B10, y) returns +0 for y>0 and not an odd integer.
F.9.4.5 The sqrt functions
[#1] sqrt is fully specified as a basic arithmetic operation
in IEC 60559.
F.9.5 Error and gamma functions
F.9.5.1 The erf functions
[#1]
-- erf(=B10) returns =B10.
-- erf(=B1) returns =B11.
F.9.4.4 IEC 60559 floating-point arithmetic F.9.5.1
=0C
WG14/N869 Committee Draft -- January 18, 1999 515
F.9.5.2 The erfc functions
[#1]
-- erfc(+) returns +0.
-- erfc(-) returns 2.
F.9.5.3 The lgamma functions
[#1]
-- lgamma(+) returns +.
-- lgamma(x) returns + and raises the divide-by-zero
exception if x is a negative integer or zero.
-- lgamma(-) returns +.
F.9.5.4 The tgamma functions
[#1]
-- tgamma(+) returns +.
-- tgamma(x) returns a NaN and raises the invalid
exception if x is a negative integer or zero.
-- tgamma(-) returns a NaN and raises the invalid
exception.
F.9.6 Nearest integer functions
F.9.6.1 The ceil functions
[#1]
-- ceil(x) returns x if x is =B1 or =B10.
The double version of ceil behaves as though implemented by
#include
#include
#pragma STDC FENV_ACCESS ON
double ceil(double x)
{
double result;
int save_round =3D fegetround();
fesetround(FE_UPWARD);
result =3D rint(x); // or nearbyint instead of rint
fesetround(save_round);
return result;
}
F.9.5.1 IEC 60559 floating-point arithmetic F.9.6.1
=0C
516 Committee Draft -- January 18, 1999 WG14/N869
F.9.6.2 The floor functions
[#1]
-- floor(x) returns x if x is =B1 or =B10.
See the sample implementation for ceil in F.9.6.1.
F.9.6.3 The nearbyint functions
[#1] The nearbyint functions use IEC 60559 rounding
according to the current rounding direction. They do not
raise the inexact exception if the result differs in value
from the argument.
F.9.6.4 The rint functions
[#1] The rint functions differ from the nearbyint functions
only in that they do raise the inexact exception if the
result differs in value from the argument.
-- rint(=B10) returns =B10 (for all rounding directions).
-- rint(=B1) returns =B1 (for all rounding directions).
F.9.6.5 The lrint and llrint functions
[#1] The lrint and llrint functions provide floating-to-
integer conversion as prescribed by IEC 60559. They round
according to the current rounding direction. If the rounded
value is outside the range of the return type, the numeric
result is unspecified and the invalid exception is raised.
When they raise no other exception and the result differs
from the argument, they raise the inexact exception.
F.9.6.6 The round functions
[#1] The double version of round behaves as though
implemented by
F.9.6.1 IEC 60559 floating-point arithmetic F.9.6.6
=0C
WG14/N869 Committee Draft -- January 18, 1999 517
#include
#include
#pragma STDC FENV_ACCESS ON
double round(double x)
{
double result;
fenv_t save_env;
feholdexcept(&save_env);
result =3D rint(x);
if (fetestexcept(FE_INEXACT)) {
fesetround(FE_TOWARDZERO);
result =3D rint(copysign(0.5 + fabs(x), x));
}
feupdateenv(&save_env);
return result;
}
The round functions may, but are not required to, raise the
inexact exception for non-integer numeric arguments, as this
implementation does.
F.9.6.7 The lround and llround functions
[#1] The lround and llround functions differ from the lrint
and llrint functions with the default rounding direction
just in that the lround and llround functions round halfway
cases away from zero, and may (but need not) raise the
inexact exception for non-integer arguments that round to
within the range of the return type.
F.9.6.8 The trunc functions
[#1] The trunc functions use IEC 60559 rounding toward zero
(regardless of the current rounding direction).
F.9.7 Remainder functions
F.9.7.1 The fmod functions
[#1]
-- fmod(=B10, y) returns =B10 if y is not zero.
-- fmod(x, y) returns a NaN and raises the invalid
exception if x is infinite or y is zero.
-- fmod(x, =B1) returns x if x is not infinite.
The double version of fmod behaves as though implemented by
F.9.6.6 IEC 60559 floating-point arithmetic F.9.7.1
=0C
518 Committee Draft -- January 18, 1999 WG14/N869
#include
#include
#pragma STDC FENV_ACCESS ON
double fmod(double x, double y)
{
double result;
result =3D remainder(fabs(x), (y =3D fabs(y)));
if (signbit(result)) result +=3D y;
return copysign(result, x);
}
F.9.7.2 The remainder functions
[#1] The remainder functions are fully specified as a basic
arithmetic operation in IEC 60559.
F.9.7.3 The remquo functions
[#1] The remquo functions follow the specifications for the
remainder functions. They have no further specifications
special to IEC 60559 implementations.
F.9.8 Manipulation functions
F.9.8.1 The copysign functions
[#1] copysign is specified in the Appendix to IEC 60559.
F.9.8.2 The nan functions
[#1] All IEC 60559 implementations support quiet NaNs, in
all floating formats.
F.9.8.3 The nextafter functions
[#1]
-- nextafter(x, y) raises the overflow and inexact
exceptions if x is finite and the function value is
infinite.
-- nextafter(x, y) raises the underflow and inexact
exceptions if the function value is subnormal or zero
and x!=3Dy.
F.9.7.1 IEC 60559 floating-point arithmetic F.9.8.3
=0C
WG14/N869 Committee Draft -- January 18, 1999 519
F.9.8.4 The nexttoward functions |
No additional requirements.
F.9.9 Maximum, minimum, and positive difference functions
F.9.9.1 The fdim functions
[#1] No additional requirements.
F.9.9.2 The fmax functions
[#1]
-- If just one argument is a NaN, the fmax functions
return the other argument (if both arguments are NaNs,
the functions return a NaN).
The body of the fmax function might be301)
{ return (isgreaterequal(x, y) ||
isnan(y)) ? x : y; }
F.9.9.3 The fmin functions
[#1] The fmin functions are analogous to the fmax functions.
See F.9.9.2.
F.9.10 Floating point multiply-add
F.9.10.1 The fma functions
[#1]
-- fma(x, y, z) computes the sum z plus the product x
times y, correctly rounded once.
-- fma(x, y, z) returns a NaN and optionally raises the
invalid exception if one of x and y is infinite, the
other is zero, and z is a NaN.
-- fma(x, y, z) returns a NaN and raises the invalid
exception if one of x and y is infinite, the other is
zero, and z is not a NaN.
-- fma(x, y, z) returns a NaN and raises the invalid
exception if x times y is an exact infinity and z is
also an infinity but with the opposite sign.
____________________
301Ideally, fmax would be sensitive to the sign of zero, for
example fmax(-0.0, +0.0) would return +0; however,
implementation in software might be impractical.
F.9.8.3 IEC 60559 floating-point arithmetic F.9.10.1
=0C
520 Committee Draft -- January 18, 1999 WG14/N869
--part1_1d9.f96c51c.2c767308_boundary--