[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