Generally, a user program examines or clears the accrued exception bits. The following example is a C program that examines the accrued exception flags.
Example 15 Examining the Accrued Exception Flags#include <sunmath.h>
#include <sys/ieeefp.h>
int main()
{
int code, inexact, division, underflow, overflow, invalid;
double x;
char *out;
/* cause an underflow exception */
x = max_subnormal() / 2.0;
/* this statement insures that the previous */
/* statement is not optimized away */
printf("x = %g\n",x);
/* find out which exceptions are raised */
code = ieee_flags("get", "exception", "", &out);
/* decode the return value */
inexact = (code >> fp_inexact) & 0x1;
underflow = (code >> fp_underflow) & 0x1;
division = (code >> fp_division) & 0x1;
overflow = (code >> fp_overflow) & 0x1;
invalid = (code >> fp_invalid) & 0x1;
/* "out" is the raised exception with the highest priority */
printf(" Highest priority exception is: %s\n", out);
/* The value 1 means the exception is raised, */
/* 0 means it isn't. */
printf("%d %d %d %d %d\n", invalid, overflow, division,
underflow, inexact);
ieee_retrospective_();
return 0;
}
The output from running the previous program is as follows:
demo% a.out
x = 1.11254e-308
Highest priority exception is: underflow
0 0 0 1 1
Note:IEEE floating-point exception flags raised:
Inexact; Underflow;
See the Numerical Computation Guide, ieee_flags(3M)
The same can be done from Fortran:
Example 16 Examining the Accrued Exception Flags – Fortran/*
A Fortran example that:
* causes an underflow exception
* uses ieee_flags to determine which exceptions are raised
* decodes the integer value returned by ieee_flags
* clears all outstanding exceptions
Remember to save this program in a file with the suffix .F, so that
the c preprocessor is invoked to bring in the header file
floatingpoint.h.
*/
#include <floatingpoint.h>
program decode_accrued_exceptions
double precision x
integer accrued, inx, div, under, over, inv
character*16 out
double precision d_max_subnormal
c Cause an underflow exception
x = d_max_subnormal() / 2.0
c Find out which exceptions are raised
accrued = ieee_flags('get', 'exception', '', out)
c Decode value returned by ieee_flags using bit-shift intrinsics
inx = and(rshift(accrued, fp_inexact) , 1)
under = and(rshift(accrued, fp_underflow), 1)
div = and(rshift(accrued, fp_division) , 1)
over = and(rshift(accrued, fp_overflow) , 1)
inv = and(rshift(accrued, fp_invalid) , 1)
c The exception with the highest priority is returned in "out"
print *, "Highest priority exception is ", out
c The value 1 means the exception is raised; 0 means it is not
print *, inv, over, div, under, inx
c Clear all outstanding exceptions
i = ieee_flags('clear', 'exception', 'all', out)
end
The output is as follows:
Highest priority exception is underflow 0 0 0 1 1
While it is unusual for a user program to set exception flags, it can be done. This is demonstrated in the following C example.
#include <sunmath.h>
int main()
{
int code;
char *out;
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf("could not clear exceptions\n");
if (ieee_flags("set", "exception", "division", &out) != 0)
printf("could not set exception\n");
code = ieee_flags("get", "exception", "", &out);
printf("out is: %s , fp exception code is: %X \n",
out, code);
return 0;
}
On SPARC, the output from the preceding program is:
out is: division , fp exception code is: 2
On x86, the output is:
out is: division , fp exception code is: 4
The following is a Fortran program that installs a signal handler to locate an exception, for SPARC-based systems only:
Example 17 Trap on Underflow (SPARC) program demo
c declare signal handler function
external fp_exc_hdl
double precision d_min_normal
double precision x
c set up signal handler
i = ieee_handler('set', 'common', fp_exc_hdl)
if (i.ne.0) print *, 'ieee trapping not supported here'
c cause an underflow exception (it will not be trapped)
x = d_min_normal() / 13.0
print *, 'd_min_normal() / 13.0 = ', x
c cause an overflow exception
c the value printed out is unrelated to the result
x = 1.0d300
x = x * x
print *, '1.0d300*1.0d300 = ', x
end
c
c the floating-point exception handling function
c
integer function fp_exc_hdl(sig, sip, uap)
integer sig, code, addr
character label*16
c
c The structure /siginfo/ is a translation of siginfo_t
c from <sys/siginfo.h>
c
structure /fault/
integer address
end structure
structure /siginfo/
integer si_signo
integer si_code
integer si_errno
record /fault/ fault
end structure
record /siginfo/ sip
c See <sys/machsig.h> for list of FPE codes
c Figure out the name of the SIGFPE
code = sip.si_code
if (code.eq.3) label = 'division'
if (code.eq.4) label = 'overflow'
if (code.eq.5) label = 'underflow'
if (code.eq.6) label = 'inexact'
if (code.eq.7) label = 'invalid'
addr = sip.fault.address
c Print information about the signal that happened
write (*,77) code, label, addr
77 format ('floating-point exception code ', i2, ',',
* a17, ',', ' at address ', z8 )
end
When the previous code is compiled with –f77, the output is as follows:
d_min_normal() / 13.0 = 1.7115952757748-309
floating-point exception code 4, overflow , at address 1131C
1.0d300*1.0d300 = 1.0000000000000+300
Note: IEEE floating-point exception flags raised:
Inexact; Underflow;
IEEE floating-point exception traps enabled:
overflow; division by zero; invalid operation;
See the Numerical Computation Guide, ieee_flags(3M),
ieee_handler(3M)
The following is a more complex C example on a SPARC-based system:
Example 18 Trap on Invalid, Division by 0, Overflow, Underflow, and Inexact (SPARC)/*
* Generate the 5 IEEE exceptions: invalid, division,
* overflow, underflow and inexact.
*
* Trap on any floating point exception, print a message,
* and continue.*
* Note that you could also inquire about raised exceptions by
* i = ieee("get","exception","",&out);* where out contains the name of the highest exception
* raised, and i can be decoded to find out about all the
* exceptions raised.
*/
#include <sunmath.h>
#include <signal.h>
#include <siginfo.h>
#include <ucontext.h>
extern void trap_all_fp_exc(int sig, siginfo_t *sip,
ucontext_t *uap);
int main()
{
double x, y, z;
char *out;
/*
* Use ieee_handler to establish "trap_all_fp_exc"
* as the signal handler to use whenever any floating
* point exception occurs.
*/
if (ieee_handler("set", "all", trap_all_fp_exc) != 0)
printf(" IEEE trapping not supported here.\n");
/* disable trapping (uninteresting) inexact exceptions */
if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
printf("Trap handler for inexact not cleared.\n");
/* raise invalid */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("1. Invalid: signaling_nan(0) * 2.5\n");
x = signaling_nan(0);
y = 2.5;
z = x * y;
/* raise division */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("2. Div0: 1.0 / 0.0\n");
x = 1.0;
y = 0.0;
z = x / y;
/* raise overflow */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("3. Overflow: -max_normal() - 1.0e294\n");
x = -max_normal();
y = -1.0e294;
z = x + y;
/* raise underflow */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("4. Underflow: min_normal() * min_normal()\n");
x = min_normal();
y = x;
z = x * y;
/* enable trapping on inexact exception */
if (ieee_handler("set", "inexact", trap_all_fp_exc) != 0)
printf("Could not set trap handler for inexact.\n");
/* raise inexact */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("5. Inexact: 2.0 / 3.0\n");
x = 2.0;
y = 3.0;
z = x / y;
/* don't trap on inexact */
if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
printf(" could not reset inexact trap\n");
/* check that we're not trapping on inexact anymore */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("6. Inexact trapping disabled; 2.0 / 3.0\n");
x = 2.0;
y = 3.0;
z = x / y;
/* find out if there are any outstanding exceptions */
ieee_retrospective_();
/* exit gracefully */
return 0;
}
void trap_all_fp_exc(int sig, siginfo_t *sip, ucontext_t *uap) {
char *label = "undefined";
/* see /usr/include/sys/machsig.h for SIGFPE codes */
switch (sip->si_code) {
case FPE_FLTRES:
label = "inexact";
break;
case FPE_FLTDIV:
label = "division";
break;
case FPE_FLTUND:
label = "underflow";
break;
case FPE_FLTINV:
label = "invalid";
break;
case FPE_FLTOVF:
label = "overflow";
break;
}
printf(
" signal %d, sigfpe code %d: %s exception at address %x\n",
sig, sip->si_code, label, sip->__data.__fault.__addr);
}
The output is similar to the following:
1. Invalid: signaling_nan(0) * 2.5 signal 8, sigfpe code 7: invalid exception at address 10da8 2. Div0: 1.0 / 0.0 signal 8, sigfpe code 3: division exception at address 10e44 3. Overflow: -max_normal() - 1.0e294 signal 8, sigfpe code 4: overflow exception at address 10ee8 4. Underflow: min_normal() * min_normal() signal 8, sigfpe code 5: underflow exception at address 10f80 5. Inexact: 2.0 / 3.0 signal 8, sigfpe code 6: inexact exception at address 1106c 6. Inexact trapping disabled; 2.0 / 3.0 Note: IEEE floating-point exception traps enabled: underflow; overflow; division by zero; invalid operation; See the Numerical Computation Guide, ieee_handler(3M)
The following code shows how you can use ieee_handler and the include files to modify the default result of certain exceptional situations on SPARC:
Example 19 Modifying the Default Result of Exceptional Situations/*
* Cause a division by zero exception and use the
* signal handler to substitute MAXDOUBLE (or MAXFLOAT)
* as the result.
*
* compile with the flag -Xa
*/
#include <values.h>
#include <siginfo.h>
#include <ucontext.h>
void division_handler(int sig, siginfo_t *sip, ucontext_t *uap);
int main() {
double x, y, z;
float r, s, t;
char *out;
/*
* Use ieee_handler to establish division_handler as the
* signal handler to use for the IEEE exception division.
*/
if (ieee_handler("set","division",division_handler)!=0) {
printf(" IEEE trapping not supported here.\n");
}
/* Cause a division-by-zero exception */
x = 1.0;
y = 0.0;
z = x / y;
/*
* Check to see that the user-supplied value, MAXDOUBLE,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
printf("double precision division: %g/%g = %g \n",x,y,z);
/* Cause a division-by-zero exception */
r = 1.0;
s = 0.0;
t = r / s;
/*
* Check to see that the user-supplied value, MAXFLOAT,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
printf("single precision division: %g/%g = %g \n",r,s,t);
ieee_retrospective_();
return 0;
}
void division_handler(int sig, siginfo_t *sip, ucontext_t *uap) {
int inst;
unsigned rd, mask, single_prec=0;
float f_val = MAXFLOAT;
double d_val = MAXDOUBLE;
long *f_val_p = (long *) &f_val;
/* Get instruction that caused exception. */
inst = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;
/*
* Decode the destination register. Bits 29:25 encode the
* destination register for any SPARC floating point
* instruction.
*/
mask = 0x1f;
rd = (mask & (inst >> 25));
/*
* Is this a single precision or double precision
* instruction? Bits 5:6 encode the precision of the
* opcode; if bit 5 is 1, it's sp, else, dp.
*/
mask = 0x1;
single_prec = (mask & (inst >> 5));
/* put user-defined value into destination register */
if (single_prec) {
uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[rd] =
f_val_p[0];
} else {
uap->uc_mcontext.fpregs.fpu_fr.fpu_dregs[rd/2] = d_val;
}
}
The following output is as expected:
double precision division: 1/0 = 1.79769e+308 single precision division: 1/0 = 3.40282e+38 Note: IEEE floating-point exception traps enabled: division by zero; See the Numerical Computation Guide, ieee_handler(3M)
You can use ieee_handler to force a program to abort in case of certain floating-point exceptions:
#include <floatingpoint.h>
program abort
c
ieeer = ieee_handler('set', 'division', SIGFPE_ABORT)
if (ieeer .ne. 0) print *, ' ieee trapping not supported'
r = 14.2
s = 0.0
r = r/s
c
print *, 'you should not see this; system should abort'
c
end
The following examples show how to use some of the exception handling features provided by libm. The first example is based on the following task: given a number x and coefficients a0, a1,..., aN, and b0, b1,..., bN-1, evaluate the function f(x) and its first derivative f'(x), where f() is the continued fraction:
f(x) =a0 + b0 / (x + a1 + b1 / (x + ... / (x + aN-1 + bN-1 / (x + aN))...)).
Computing f() is straightforward in IEEE arithmetic: even if one of the intermediate divisions overflows or divides by zero, the default value specified by the standard (a correctly signed infinity) turns out to yield the correct result. Computing f'(), on the other hand, can be more difficult because the simplest form for evaluating it can have removable singularities. If the computation encounters one of these singularities, it will attempt to evaluate one of the indeterminate forms 0/0, 0*infinity, or infinity/infinity, all of which raise invalid operation exceptions. W. Kahan has proposed a method for handling these exceptions via a feature called presubstitution.
Presubstitution is an extension of the IEEE default response to exceptions that lets the user specify in advance the value to be substituted for the result of an exceptional operation. Using the exception handling facilities in libm, a program can implement presubstitution easily by installing a handler in the FEX_CUSTOM exception handling mode. This mode allows the handler to supply any value for the result of an exceptional operation simply by storing that value in the data structure pointed to by the info parameter passed to the handler. The following example is a sample program to compute the continued fraction and its derivative using presubstitution implemented with a FEX_CUSTOM handler.
Example 20 Computing the Continued Fraction and Its Derivative Using the FEX_CUSTOM Handler#include <stdio.h>
#include <sunmath.h>
#include <fenv.h>
volatile double p;
void handler(int ex, fex_info_t *info)
{
info->res.type = fex_double;
if (ex == FEX_INV_ZMI)
info->res.val.d = p;
else
info->res.val.d = infinity();
}
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *pf1
*/
void continued_fraction(int N, double *a, double *b,
double x, double *pf, double *pf1)
{
fex_handler_t oldhdl; /* for saving/restoring handlers */
volatile double t;
double f, f1, d, d1, q;
int j;
fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
fex_set_handling(FEX_DIVBYZERO, FEX_NONSTOP, NULL);
fex_set_handling(FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI,
FEX_CUSTOM, handler);
f1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
/* the following assignment to the volatile variable t
is needed to maintain the correct sequencing between
assignments to p and evaluation of f1 */
t = f1 = (-d1 / d) * q;
p = b[j-1] * d1 / b[j];
f = a[j] + q;
}
fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
*pf = f;
*pf1 = f1;
}
/* For the following coefficients, x = -3, 1, 4, and 5 will all
encounter intermediate exceptions */
double a[] = { -1.0, 2.0, -3.0, 4.0, -5.0 };
double b[] = { 2.0, 4.0, 6.0, 8.0 };
int main()
{
double x, f, f1;
int i;
feraiseexcept(FE_INEXACT); /* prevent logging of inexact */
fex_set_log(stdout);
fex_set_handling(FEX_COMMON, FEX_ABORT, NULL);
for (i = -5; i <= 5; i++) {
x = i;
continued_fraction(4, a, b, x, &f, &f1);
printf("f(% g) = %12g, f'(% g) = %12g\n", x, f, x, f1);
}
return 0;
}
Several comments about the program are in order. On entry, the function continued_fraction saves the current exception handling modes for division by zero and all invalid operation exceptions. It then establishes nonstop exception handling for division by zero and a FEX_CUSTOM handler for the three indeterminate forms. This handler will substitute infinity for both 0/0 and infinity/infinity, but it will substitute the value of the global variable p for 0*infinity. Note that p must be recomputed each time through the loop that evaluates the function in order to supply the correct value to substitute for a subsequent 0*infinity invalid operation. Note also that p must be declared volatile to prevent the compiler from eliminating it, since it is not explicitly mentioned elsewhere in the loop. Finally, to prevent the compiler from moving the assignment to p above or below the computation that can incur the exception for which p provides the presubstitution value, the result of that computation is also assigned to a volatile variable (called t in the program). The final call to fex_setexcepthandler restores the original handling modes for division by zero and the invalid operations.
The main program enables logging of retrospective diagnostics by calling the fex_set_log function. Before it does so, it raises the inexact flag; this has the effect of preventing the logging of inexact exceptions. Recall that in FEX_NONSTOP mode, an exception is not logged if its flag is raised, as explained in the section Retrospective Diagnostics. The main program also establishes FEX_ABORT mode for the common exceptions to ensure that any unusual exceptions not explicitly handled by continued_fraction will cause program termination. Finally, the program evaluates a particular continued fraction at several different points. As the following sample output shows, the computation does indeed encounter intermediate exceptions:
f(-5) = -1.59649, f'(-5) = -0.1818 f(-4) = -1.87302, f'(-4) = -0.428193 Floating point division by zero at 0x08048dbe continued_fraction, nonstop mode 0x08048dc1 continued_fraction 0x08048eda main Floating point invalid operation (inf/inf) at 0x08048dcf continued_fraction, handler: handler 0x08048dd2 continued_fraction 0x08048eda main Floating point invalid operation (0*inf) at 0x08048dd2 continued_fraction, handler: handler 0x08048dd8 continued_fraction 0x08048eda main f(-3) = -3, f'(-3) = -3.16667 f(-2) = -4.44089e-16, f'(-2) = -3.41667 f(-1) = -1.22222, f'(-1) = -0.444444 f( 0) = -1.33333, f'( 0) = 0.203704 f( 1) = -1, f'( 1) = 0.333333 f( 2) = -0.777778, f'( 2) = 0.12037 f( 3) = -0.714286, f'( 3) = 0.0272109 f( 4) = -0.666667, f'( 4) = 0.203704 f( 5) = -0.777778, f'( 5) = 0.0185185
The exceptions that occur in the computation of f'(x) at x = 1, 4, and 5 do not result in retrospective diagnostic messages because they occur at the same site in the program as the exceptions that occur when x = –3.
The preceding program might not represent the most efficient way to handle the exceptions that can occur in the evaluation of a continued fraction and its derivative. One reason is that the presubstitution value must be recomputed in each iteration of the loop regardless of whether or not it is needed. In this case, the computation of the presubstitution value involves a floating-point division, and on modern SPARC and x86 processors, floating-point division is a relatively slow operation. Moreover, the loop itself already involves two divisions, and because most SPARC and x86 processors cannot overlap the execution of two different division operations, divisions are likely to be a bottleneck in the loop; adding another division would exacerbate the bottleneck.
It is possible to rewrite the loop so that only one division is needed, and in particular, the computation of the presubstitution value need not involve a division. To rewrite the loop in this way, one must precompute the ratios of adjacent elements of the coefficients in the b array. This would remove the bottleneck of multiple division operations, but it would not eliminate all of the arithmetic operations involved in the computation of the presubstitution value. Furthermore, the need to assign both the presubstitution value and the result of the operation to be presubstituted to volatile variables introduces additional memory operations that slow the program. While those assignments are necessary to prevent the compiler from reordering certain key operations, they effectively prevent the compiler from reordering other unrelated operations, too. Thus, handling the exceptions in this example via presubstitution requires additional memory operations and precludes some optimizations that might otherwise be possible. Can these exceptions be handled more efficiently?
In the absence of special hardware support for fast presubstitution, the most efficient way to handle exceptions in this example may be to use flags, as the following version does:
Example 21 Using Flags to Handle Exceptions (Continued)#include <stdio.h>
#include <math.h>
#include <fenv.h>
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *pf1
*/
void continued_fraction(int N, double *a, double *b,
double x, double *pf, double *pf1)
{
fex_handler_t oldhdl;
fexcept_t oldinvflag;
double f, f1, d, d1, pd1, q;
int j;
fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
fegetexceptflag(&oldinvflag, FE_INVALID);
fex_set_handling(FEX_DIVBYZERO | FEX_INV_ZDZ | FEX_INV_IDI |
FEX_INV_ZMI, FEX_NONSTOP, NULL);
feclearexcept(FE_INVALID);
f1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
f1 = (-d1 / d) * q;
f = a[j] + q;
}
if (fetestexcept(FE_INVALID)) {
/* recompute and test for NaN */
f1 = pd1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
f1 = (-d1 / d) * q;
if (isnan(f1))
f1 = b[j] * pd1 / b[j+1];
pd1 = d1;
f = a[j] + q;
}
}
fesetexceptflag(&oldinvflag, FE_INVALID);
fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
*pf = f;
*pf1 = f1;
}
In this version, the first loop attempts the computation of f(x) and f'(x) in the default nonstop mode. If the invalid flag is raised, the second loop recomputes f(x) and f'(x) explicitly testing for the appearance of a NaN. Usually, no invalid operation exception occurs, so the program only executes the first loop. This loop has no references to volatile variables and no extra arithmetic operations, so it will run as fast as the compiler can make it go. The cost of this efficiency is the need to write a second loop nearly identical to the first to handle the case when an exception occurs. This trade-off is typical of the dilemmas that floating-point exception handling can pose.
The exception handling facilities in libm are primarily intended to be used from C/C++ programs, but by using the Sun Fortran language interoperability features, you can call some libm functions from Fortran programs as well.
The following example shows a Fortran version of the program to evaluate a continued fraction and its derivative using presubstitution (SPARC only):
Example 22 Evaluating a Continued Fraction and Its Derivative Using Presubstitution (SPARC)c
c Presubstitution handler
c
subroutine handler(ex, info)
structure /fex_numeric_t/
integer type
union
map
integer i
end map
map
integer*8 l
end map
map
real f
end map
map
real*8 d
end map
map
real*16 q
end map
end union
end structure
structure /fex_info_t/
integer op, flags
record /fex_numeric_t/ op1, op2, res
end structure
integer ex
record /fex_info_t/ info
common /presub/ p
double precision p, d_infinity
volatile p
c 4 = fex_double; see <fenv.h> for this and other constants
info.res.type = 4
c x'80' = FEX_INV_ZMI
if (loc(ex) .eq. x'80') then
info.res.d = p
else
info.res.d = d_infinity()
endif
return
end
c
c Evaluate the continued fraction given by coefficients a(j) and
c b(j) at the point x; return the function value in f and the
c derivative in f1
c
subroutine continued_fraction(n, a, b, x, f, f1)
integer n
double precision a(*), b(*), x, f, f1
common /presub/ p
integer j, oldhdl
dimension oldhdl(24)
double precision d, d1, q, p, t
volatile p, t
data ixff2/x'ff2'/
data ix2/x'2'/
data ixb0/x'b0'/
external fex_getexcepthandler, fex_setexcepthandler
external fex_set_handling, handler
c$pragma c(fex_getexcepthandler, fex_setexcepthandler)
c$pragma c(fex_set_handling)
c x'ff2' = FEX_DIVBYZERO | FEX_INVALID
call fex_getexcepthandler(oldhdl, %val(ixff2))
c x'2' = FEX_DIVBYZERO, 0 = FEX_NONSTOP
call fex_set_handling(%val(ix2), %val(0), %val(0))
c x'b0' = FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI, 3 = FEX_CUSTOM
call fex_set_handling(%val(ixb0), %val(3), handler)
f1 = 0.0d0
f = a(n+1)
do j = n, 1, -1
d = x + f
d1 = 1.0d0 + f1
q = b(j) / d
f1 = (-d1 / d) * q
c
c the following assignment to the volatile variable t
c is needed to maintain the correct sequencing between
c assignments to p and evaluation of f1
t = f1
p = b(j-1) * d1 / b(j)
f = a(j) + q
end do
call fex_setexcepthandler(oldhdl, %val(ixff2))
return
end
c Main program
c
program cf
integer i
double precision a, b, x, f, f1
dimension a(5), b(4)
data a /-1.0d0, 2.0d0, -3.0d0, 4.0d0, -5.0d0/
data b /2.0d0, 4.0d0, 6.0d0, 8.0d0/
data ixffa/x'ffa'/
external fex_set_handling
c$pragma c(fex_set_handling)
c x'ffa' = FEX_COMMON, 1 = FEX_ABORT
call fex_set_handling(%val(ixffa), %val(1), %val(0))
do i = -5, 5
x = dble(i)
call continued_fraction(4, a, b, x, f, f1)
write (*, 1) i, f, i, f1
end do
1 format('f(', I2, ') = ', G12.6, ', f''(', I2, ') = ', G12.6)
end
The output from this program compiled with the –f77 flag reads as follows:
f(-5) = -1.59649 , f'(-5) = -.181800
f(-4) = -1.87302 , f'(-4) = -.428193
f(-3) = -3.00000 , f'(-3) = -3.16667
f(-2) = -.444089E-15, f'(-2) = -3.41667
f(-1) = -1.22222 , f'(-1) = -.444444
f( 0) = -1.33333 , f'( 0) = 0.203704
f( 1) = -1.00000 , f'( 1) = 0.333333
f( 2) = -.777778 , f'( 2) = 0.120370
f( 3) = -.714286 , f'( 3) = 0.272109E-01
f( 4) = -.666667 , f'( 4) = 0.203704
f( 5) = -.777778 , f'( 5) = 0.185185E-01
Note: IEEE floating-point exception flags raised:
Inexact; Division by Zero; Underflow; Invalid Operation;
IEEE floating-point exception traps enabled:
overflow; division by zero; invalid operation;
See the Numerical Computation Guide, ieee_flags(3M),
ieee_handler(3M)