Historically, most numerical software has been written without regard to exceptions, and many programmers have become accustomed to environments in which exceptions cause a program to abort immediately. High-quality software packages such as LAPACK are carefully designed to avoid exceptions such as division by zero and invalid operations and to scale their inputs aggressively to preclude overflow and potentially harmful underflow. Neither of these approaches to dealing with exceptions is appropriate in every situation. However, ignoring exceptions can pose problems when you write a program or subroutine that is intended to be used by someone else, for example, who might not have access to the source code. Attempting to avoid all exceptions can require many defensive tests and branches and carry a significant cost. For more information, see Demmel and Li, “Faster Numerical Algorithms via Exception Handling,” IEEE Trans. Comput. 43 (1994), pp. 983–992.
The default exception response, status flags, and optional trapping facility of IEEE arithmetic are intended to provide a third alternative: continuing a computation in the presence of exceptions and either detecting them after the fact or intercepting and handling them as they occur. As described above, ieee_flags or the C99 floating-point environment functions can be used to detect exceptions after the fact, and ieee_handler or fex_set_handling can be used to enable trapping and install a handler to intercept exceptions as they occur. In order to continue the computation, however, the IEEE standard recommends that a trap handler be able to provide a result for the operation that incurred an exception. A SIGFPE handler installed via ieee_handler or fex_set_handling in FEX_SIGNAL mode can accomplish this using the uap parameter supplied to a signal handler by the Solaris operating environment. An FEX_CUSTOM mode handler installed via fex_set_handling can provide a result using the info parameter supplied to such a handler.
Recall that a SIGFPE signal handler can be declared in C as follows:
#include <siginfo.h>
#include <ucontext.h>
void handler(int sig, siginfo_t *sip, ucontext_t *uap)
{
...
}
When a SIGFPE signal handler is invoked as a result of a trapped floating-point exception, the uap parameter points to a data structure that contains a copy of the machine's integer and floating-point registers as well as other system-dependent information describing the exception. If the signal handler returns normally, the saved data are restored and the program resumes execution at the point at which the trap was taken. Thus, by accessing and decoding the information in the data structure that describes the exception and possibly modifying the saved data, a SIGFPE handler can substitute a user-supplied value for the result of an exceptional operation and continue computation.
An FEX_CUSTOM mode handler can be declared as follows:
#include <fenv.h>
void handler(int ex, fex_info_t *info)
{
...
}
When a FEX_CUSTOM handler is invoked, the ex parameter indicates which type of exception occurred (it is one of the values listed in Table 34) and the info parameter points to a data structure that contains more information about the exception. Specifically, this structure contains a code representing the arithmetic operation that caused the exception and structures recording the operands, if they are available. It also contains a structure recording the default result that would have been substituted if the exception were not trapped and an integer value holding the bitwise “or” of the exception flags that would have accrued. The handler can modify the latter members of the structure to substitute a different result or change the set of flags that are accrued. (Note that if the handler returns without modifying these data, the program will continue with the default untrapped result and flags just as if the exception were not trapped.)
As an illustration, the following section shows how to substitute a scaled result for an operation that underflows or overflows. See Examples for further examples.
The IEEE standard recommends that when underflow and overflow are trapped, the system should provide a way for a trap handler to substitute an exponent-wrapped result, i.e., a value that agrees with what would have been the rounded result of the operation that underflowed or overflowed except that the exponent is wrapped around the end of its usual range, thereby effectively scaling the result by a power of two. The scale factor is chosen to map underflowed and overflowed results as nearly as possible to the middle of the exponent range so that subsequent computations will be less likely to underflow or overflow further. By keeping track of the number of underflows and overflows that occur, a program can scale the final result to compensate for the exponent wrapping. This under/overflow “counting mode” can be used to produce accurate results in computations that would otherwise exceed the range of the available floating-point formats. See P. Sterbenz, Floating-Point Computation for more information.
On SPARC-based systems, when a floating-point instruction incurs a trapped exception, the system leaves the destination register unchanged. Thus, in order to substitute the exponent-wrapped result, an under/overflow handler must decode the instruction, examine the operand registers, and generate the scaled result itself. The following example shows a handler that performs these steps.
Example 4 Substituting IEEE Trapped Under/Overflow Handler Results for SPARC-Based Systems#include <stdio.h>
#include <ieeefp.h>
#include <math.h>
#include <sunmath.h>
#include <siginfo.h>
#include <ucontext.h>
#ifdef V8PLUS
/* The upper 32 floating point registers are stored in an area
pointed to by uap->uc_mcontext.xrs.xrs_prt. Note that this
pointer is valid ONLY when uap->uc_mcontext.xrs.xrs_id ==
XRS_ID (defined in sys/procfs.h). */
#include <assert.h>
#include <sys/procfs.h>
#define FPxreg(x) ((prxregset_t*)uap->uc_mcontext.xrs.xrs_ptr)
->pr_un.pr_v8p.pr_xfr.pr_regs[(x)]
#endif
#define FPreg(x) uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[(x)]
/*
* Supply the IEEE 754 default result for trapped under/overflow
*/
void
ieee_trapped_default(int sig, siginfo_t *sip, ucontext_t *uap)
{
unsigned instr, opf, rs1, rs2, rd;
long double qs1, qs2, qd, qscl;
double ds1, ds2, dd, dscl;
float fs1, fs2, fd, fscl;
/* get the instruction that caused the exception */
instr = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;
/* extract the opcode and source and destination register
numbers */
opf = (instr >> 5) & 0x1ff;
rs1 = (instr >> 14) & 0x1f;
rs2 = instr & 0x1f;
rd = (instr >> 25) & 0x1f;
/* get the operands */
switch (opf & 3) {
case 1: /* single precision */
fs1 = *(float*)&FPreg(rs1);
fs2 = *(float*)&FPreg(rs2);
break;
case 2: /* double precision */
#ifdef V8PLUS
if (rs1 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
ds1 = *(double*)&FPxreg(rs1 & 0x1e);
}
else
ds1 = *(double*)&FPreg(rs1);
if (rs2 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
ds2 = *(double*)&FPxreg(rs2 & 0x1e);
}
else
ds2 = *(double*)&FPreg(rs2);
#else
ds1 = *(double*)&FPreg(rs1);
ds2 = *(double*)&FPreg(rs2);
#endif
break;
case 3: /* quad precision */
#ifdef V8PLUS
if (rs1 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
qs1 = *(long double*)&FPxreg(rs1 & 0x1e);
}
else
qs1 = *(long double*)&FPreg(rs1);
if (rs2 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
qs2 = *(long double*)&FPxreg(rs2 & 0x1e);
}
else
qs2 = *(long double*)&FPreg(rs2);
#else
qs1 = *(long double*)&FPreg(rs1);
qs2 = *(long double*)&FPreg(rs2);
#endif
break;
}
/* set up scale factors */
if (sip->si_code == FPE_FLTOVF) {
fscl = scalbnf(1.0f, -96);
dscl = scalbn(1.0, -768);
qscl = scalbnl(1.0, -12288);
} else {
fscl = scalbnf(1.0f, 96);
dscl = scalbn(1.0, 768);
qscl = scalbnl(1.0, 12288);
}
/* disable traps and generate the scaled result */
fpsetmask(0);
switch (opf) {
case 0x41: /* add single */
fd = fscl * (fscl * fs1 + fscl * fs2);
break;
case 0x42: /* add double */
dd = dscl * (dscl * ds1 + dscl * ds2);
break;
case 0x43: /* add quad */
qd = qscl * (qscl * qs1 + qscl * qs2);
break;
case 0x45: /* subtract single */
fd = fscl * (fscl * fs1 - fscl * fs2);
break;
case 0x46: /* subtract double */
dd = dscl * (dscl * ds1 - dscl * ds2);
break;
case 0x47: /* subtract quad */
qd = qscl * (qscl * qs1 - qscl * qs2);
break;
case 0x49: /* multiply single */
fd = (fscl * fs1) * (fscl * fs2);
break;
case 0x4a: /* multiply double */
dd = (dscl * ds1) * (dscl * ds2);
break;
case 0x4b: /* multiply quad */
qd = (qscl * qs1) * (qscl * qs2);
break;
case 0x4d: /* divide single */
fd = (fscl * fs1) / (fs2 / fscl);
break;
case 0x4e: /* divide double */
dd = (dscl * ds1) / (ds2 / dscl);
break;
case 0x4f: /* divide quad */
qd = (qscl * qs1) / (qs2 / dscl);
break;
case 0xc6: /* convert double to single */
fd = (float) (fscl * (fscl * ds1));
break;
case 0xc7: /* convert quad to single */
fd = (float) (fscl * (fscl * qs1));
break;
case 0xcb: /* convert quad to double */
dd = (double) (dscl * (dscl * qs1));
break;
}
/* store the result in the destination */
if (opf & 0x80) {
/* conversion operation */
if (opf == 0xcb) {
/* convert quad to double */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(double*)&FPxreg(rd & 0x1e) = dd;
}
else
*(double*)&FPreg(rd) = dd;
#else
*(double*)&FPreg(rd) = dd;
#endif
} else
/* convert quad/double to single */
*(float*)&FPreg(rd) = fd;
} else {
/* arithmetic operation */
switch (opf & 3) {
case 1: /* single precision */
*(float*)&FPreg(rd) = fd;
break;
case 2: /* double precision */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(double*)&FPxreg(rd & 0x1e) = dd;
}
else
*(double*)&FPreg(rd) = dd;
#else
*(double*)&FPreg(rd) = dd;
#endif
break;
case 3: /* quad precision */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(long double*)&FPxreg(rd & 0x1e) = qd;
}
else
*(long double*)&FPreg(rd & 0x1e) = qd;
#else
*(long double*)&FPreg(rd & 0x1e) = qd;
#endif
break;
}
}
}
int
main()
{
volatile float a, b;
volatile double x, y;
ieee_handler("set", "underflow", ieee_trapped_default);
ieee_handler("set", "overflow", ieee_trapped_default);
a = b = 1.0e30f;
a *= b; /* overflow; will be wrapped to a moderate number */
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
a /= b; /* underflow; will wrap back */
printf( "%g\n", a );
x = y = 1.0e300;
x *= y; /* overflow; will be wrapped to a moderate number */
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
x /= y; /* underflow; will wrap back */
printf( "%g\n", x );
ieee_retrospective(stdout);
return 0;
}
In this example, the variables a, b, x, and y have been declared volatile only to prevent the compiler from evaluating a * b, etc., at compile time. In typical usage, the volatile declarations would not be needed.
The output from the preceding program is as follows:
159.309
1.59309e-28
1
4.14884e+137
4.14884e-163
1
Note: IEEE floating-point exception traps enabled:
underflow; overflow;
See the Numerical Computation Guide, ieee_handler(3M)
On x86-based systems, the floating-point hardware provides the exponent-wrapped result when a floating-point instruction incurs a trapped underflow or overflow and its destination is a register. When trapped underflow or overflow occurs on a floating-point store instruction, however, the hardware traps without completing the store and without popping the stack, if the store instruction is a store-and-pop. Thus, in order to implement counting mode, an under/overflow handler must generate the scaled result and fix up the stack when a trap occurs on a store instruction. Example 5 illustrates such a handler.
Example 5 Substituting IEEE Trapped Under/Overflow Handler Results for x86-Based Systems#include <stdio.h>
#include <ieeefp.h>
#include <math.h>
#include <sunmath.h>
#include <siginfo.h>
#include <ucontext.h>
/* offsets into the saved fp environment */
#define CW 0 /* control word */
#define SW 1 /* status word */
#define TW 2 /* tag word */
#define OP 4 /* opcode */
#define EA 5 /* operand address */
#define FPenv(x) uap->uc_mcontext.fpregs.fp_reg_set.
fpchip_state.state[(x)]
#define FPreg(x) *(long double *)(10*(x)+(char*)&uap->
uc_mcontext.fpregs.fp_reg_set.fpchip_state.state[7])/*
* Supply the IEEE 754 default result for trapped under/overflow
*/
void
ieee_trapped_default(int sig, siginfo_t *sip, ucontext_t *uap)
{
double dscl;
float fscl;
unsigned sw, op, top;
int mask, e;
/* preserve flags for untrapped exceptions */
sw = uap->uc_mcontext.fpregs.fp_reg_set.fpchip_state.status;
FPenv(SW) |= (sw & (FPenv(CW) & 0x3f));
/* if the excepting instruction is a store, scale the stack
top, store it, and pop the stack if need be */
fpsetmask(0);
op = FPenv(OP) >> 16;
switch (op & 0x7f8) {
case 0x110:
case 0x118:
case 0x150:
case 0x158:
case 0x190:
case 0x198:
fscl = scalbnf(1.0f, (sip->si_code == FPE_FLTOVF)?
-96 : 96);
*(float *)FPenv(EA) = (FPreg(0) * fscl) * fscl;
if (op & 8) {
/* pop the stack */
FPreg(0) = FPreg(1);
FPreg(1) = FPreg(2);
FPreg(2) = FPreg(3);
FPreg(3) = FPreg(4);
FPreg(4) = FPreg(5);
FPreg(5) = FPreg(6);
FPreg(6) = FPreg(7);
top = (FPenv(SW) >> 10) & 0xe;
FPenv(TW) |= (3 << top);
top = (top + 2) & 0xe;
FPenv(SW) = (FPenv(SW) & ~0x3800) | (top << 10);
}
break;
case 0x510:
case 0x518:
case 0x550:
case 0x558:
case 0x590:
case 0x598:
dscl = scalbn(1.0, (sip->si_code == FPE_FLTOVF)?
-768 : 768);
*(double *)FPenv(EA) = (FPreg(0) * dscl) * dscl;
if (op & 8) {
/* pop the stack */
FPreg(0) = FPreg(1);
FPreg(1) = FPreg(2);
FPreg(2) = FPreg(3);
FPreg(3) = FPreg(4);
FPreg(4) = FPreg(5);
FPreg(5) = FPreg(6);
FPreg(6) = FPreg(7);
top = (FPenv(SW) >> 10) & 0xe;
FPenv(TW) |= (3 << top);
top = (top + 2) & 0xe;
FPenv(SW) = (FPenv(SW) & ~0x3800) | (top << 10);
}
break;
}
}
int main()
{
volatile float a, b;
volatile double x, y;
ieee_handler("set", "underflow", ieee_trapped_default);
ieee_handler("set", "overflow", ieee_trapped_default);
a = b = 1.0e30f;
a *= b;
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
x = y = 1.0e300;
x *= y;
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
ieee_retrospective(stdout);
return 0;
}
As on SPARC-based systems and compiled with –xarch=386, the output from the preceding program on x86 is:
159.309
1.59309e-28
1
4.14884e+137
4.14884e-163
1
Note: IEEE floating-point exception traps enabled:
underflow; overflow;
See the Numerical Computation Guide, ieee_handler(3M)
C/C++ programs can use the fex_set_handling function in libm to install a FEX_CUSTOM handler for underflow and overflow. On SPARC-based systems, the information supplied to such a handler always includes the operation that caused the exception and the operands, and this information is sufficient to allow the handler to compute the IEEE exponent-wrapped result, as shown above. On x86-based systems, the available information might not always indicate which particular operation caused the exception; when the exception is raised by one of the transcendental instructions, for example, the info->op parameter is set to fex_other. (See the fenv.h file for definitions.) Moreover, the x86 hardware delivers an exponent-wrapped result automatically, and this can overwrite one of the operands if the destination of the excepting instruction is a floating-point register.
Fortunately, the fex_set_handling feature provides a simple way for a handler installed in FEX_CUSTOM mode to substitute the IEEE exponent-wrapped result for an operation that underflows or overflows. When either of these exceptions is trapped, the handler can set
info->res.type = fex_nodata;
to indicate that the exponent-wrapped result should be delivered. The following is an example showing such a handler:
#include <stdio.h>
#include <fenv.h>
void handler(int ex, fex_info_t *info) {
info->res.type = fex_nodata;
}
int main()
{
volatile float a, b;
volatile double x, y;
fex_set_log(stderr);
fex_set_handling(FEX_UNDERFLOW | FEX_OVERFLOW, FEX_CUSTOM,
handler);
a = b = 1.0e30f;
a *= b; /* overflow; will be wrapped to a moderate number */
printf("%g\n", a);
a /= b;
printf("%g\n", a);
a /= b; /* underflow; will wrap back */
printf("%g\n", a);
x = y = 1.0e300;
x *= y; /* overflow; will be wrapped to a moderate number */
printf("%g\n", x);
x /= y;
printf("%g\n", x);
x /= y; /* underflow; will wrap back */
printf("%g\n", x);
return 0;
}
The output from the preceding program resembles the following:
Floating point overflow at 0x00010924 main, handler: handler 0x00010928 main 159.309 1.59309e-28 Floating point underflow at 0x00010994 main, handler: handler 0x00010998 main 1 Floating point overflow at 0x000109e4 main, handler: handler 0x000109e8 main 4.14884e+137 4.14884e-163 Floating point underflow at 0x00010a4c main, handler: handler 0x00010a50 main 1
The previous program example works on SPARC, on x86 with –xarch=sse2, and on x86 with –xarch=386.