[Cfp-interest] continued fraction example with fe_try fe_catch
David Hough CFP
pcfp at oakapple.net
Thu May 29 14:52:51 PDT 2014
Let's try fe_try fe_catch on the whole loop, and do the funny business
only when needed
fe_try {
for (j = N - 1; j >= 0; j--) {
d = x + f;
q = b[j] / d;
f = a[j] + q;
d1 = 1.0 + f1;
/* CRITICAL STEP */
/* division trapped if invalid */
f1 = -(d1 / d) * q;
}
}
fe_catch_excep( int e, FE_INVALID) {
/* exceptional case when invalid division arose though it also catches
signaling NaNs etc. catching FE_INVALID_DIV looks cleaner but is
apt to be much slower */
for (j = N - 1; j >= 0; j--) {
f1j2 = f1j1; /* Save this guy for recomputation -
won't be defined until j=N-2 but that's the
earliest it will be needed.. */
f1j1 = f1; /* Intermediate save. */
d = x + f;
q = b[j] / d;
f = a[j] + q;
d1 = 1.0 + f1;
/* CRITICAL STEP */
r = d1 / d
if (isnan(r) == 0) {
/* r is a normal number */
f1 = -r * q;
} else {
/* r is NaN from 0/0 or inf/inf */
if (d1 == 0) {
/* f1 = (0/0)*inf so return infinity */
f1 = q;
} else {
/* f1 = (inf/inf)*0 so return limit value */
f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];
}
}
}
}
OK, so what about presubstitution? The first try loop remains the
same, but the catch loop becomes itself a try loop and contains a third
level:
fe_try {
for (j = N - 1; j >= 0; j--) {
d = x + f;
q = b[j] / d;
f = a[j] + q;
d1 = 1.0 + f1;
/* CRITICAL STEP */
/* division trapped if invalid */
f1 = -(d1 / d) * q;
}
}
fe_catch_excep( int e, FE_INVALID) {
/* exceptional case when invalid division arose though it also catches
signaling NaNs etc. catching FE_INVALID_DIV looks cleaner but is
apt to be much slower */
fe_try {
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
bj = b[j];
q = bj / d;
r = d1 / d ;
/* CRITICAL STEP - third level fe_try - p not defined in this try block */
fe_try { f1 = -r * q ; }
fe_catch_sub(int e, FE_INVALID_MUL) { return(p); }
f = a[j] + q;
p = bj1 *d1 / bj ;
bj1 = bj ;
}
}
}
fe_catch_sub (int e, FE_INVALID_DIV) {
/* catch 0/0 and inf/inf cases computing q or r -
both can be infinity here */
/* have to catch subexception because FE_INVALID_MUL is handled
differently */
return( INFINITY ); /* "return" is just decorative */
}
}
I'm not convinced that nesting three layers of try/catch is a win for
debuggability and provability.
More information about the Cfp-interest
mailing list