<html><body>
<p><font size="2" face="sans-serif">An interesting statement is</font><br>
<tt><font size="2"> presubstitute() needs to be faster than a conditional branch in order to be useful. Thus it shouldn't be a synchronizing operation that stops the pipe</font></tt><br>
<br>
<font size="2" face="sans-serif">On many architectures that really means faster than a floating-point compare and conditional branch.</font><br>
<br>
<font size="2" face="sans-serif">One way to beat that would be by the suggestion of having "substitution registers" containing values to be substituted for the normal results whenever an exception occurs. Another would be by having an exception trap to a handler whenever an exception occurs, with zero cost otherwise. On many systems initializing and checking control and status registers is a synchronizing operation that stops the floating-point pipeline so wouldn't meet the performance goal.</font><br>
<br>
<font size="2" face="sans-serif">Exception handling may be faster for other examples, and could also be useful because it is so much more usable than redesigning algorithms.</font><br>
<br>
<font size="2" face="sans-serif">- Ian McIntosh IBM Canada Lab Compiler Back End Support and Development<br>
</font><br>
<br>
<img width="16" height="16" src="cid:1__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt="Inactive hide details for David Hough CFP ---2014-05-13 04:35:23 PM---Sorry I forgot to actually include the example last time:"><font size="2" color="#424282" face="sans-serif">David Hough CFP ---2014-05-13 04:35:23 PM---Sorry I forgot to actually include the example last time: ==========</font><br>
<br>
<table width="100%" border="0" cellspacing="0" cellpadding="0">
<tr valign="top"><td width="1%"><img width="96" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<ul style="padding-left: 4pt"><font size="1" color="#5F5F5F" face="sans-serif">From:</font></ul>
</td><td width="100%"><img width="1" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<font size="1" face="sans-serif">David Hough CFP <pcfp@oakapple.net></font></td></tr>
<tr valign="top"><td width="1%"><img width="96" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<ul style="padding-left: 4pt"><font size="1" color="#5F5F5F" face="sans-serif">To:</font></ul>
</td><td width="100%"><img width="1" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<font size="1" face="sans-serif">cfp-interest@ucbtest.org, </font></td></tr>
<tr valign="top"><td width="1%"><img width="96" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<ul style="padding-left: 4pt"><font size="1" color="#5F5F5F" face="sans-serif">Date:</font></ul>
</td><td width="100%"><img width="1" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<font size="1" face="sans-serif">2014-05-13 04:35 PM</font></td></tr>
<tr valign="top"><td width="1%"><img width="96" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<ul style="padding-left: 4pt"><font size="1" color="#5F5F5F" face="sans-serif">Subject:</font></ul>
</td><td width="100%"><img width="1" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<font size="1" face="sans-serif">[Cfp-interest] exception handling example included this time</font></td></tr>
<tr valign="top"><td width="1%"><img width="96" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<ul style="padding-left: 4pt"><font size="1" color="#5F5F5F" face="sans-serif">Sent by:</font></ul>
</td><td width="100%"><img width="1" height="1" src="cid:2__=0ABBF64BDFFF51C38f9e8a93df938@ca.ibm.com" border="0" alt=""><br>
<font size="1" face="sans-serif">cfp-interest-bounces@oakapple.net</font></td></tr>
</table>
<hr width="100%" size="2" align="left" noshade style="color:#8091A5; "><br>
<br>
<br>
<tt><font size="2"><br>
Sorry I forgot to actually include the example last time:<br>
<br>
==========<br>
<br>
From dgh Tue Oct 8 14:56:16 1991<br>
To: nceg@cray.com<br>
Subject: Exception Handling I: The Critical Example<br>
<br>
In his proposals for floating-point exception handling, W. Kahan uses the<br>
following example of computing a continued fraction and its derivative.<br>
It's a critical case for exception handling proposals because the desired<br>
exception handling for the 0*inf case changes on each iteration.<br>
Here's the underlying code with no exception handling:<br>
<br>
void<br>
continued_fraction(N, a, b, x, pf, pf1)<br>
int N;<br>
double *a, *b, x, *pf, *pf1;<br>
{<br>
double f, f1, d, d1, q;<br>
int j;<br>
<br>
/*<br>
* Evaluate function f at x by continued fraction a[j] and b[j];<br>
* function value to *pf, derivative to *pf1<br>
*/<br>
<br>
/* Assume<br>
aj finite<br>
bj finite and nonzero<br>
x finite<br>
<br>
CRITICAL STEP is evaluating f1= -(d1/d)*q<br>
<br>
If d == 0 and d1 != 0 on step j, then<br>
on the next step j-1 you want f1 = b(j-1)/b(j) * d1(j),<br>
a result determined by taking limits as d->0 in step j.<br>
Instead by IEEE default, <br>
on the next step j-1 you will get f1 = -(inf/inf)*0<br>
or NaN.<br>
<br>
if d == 0 and d1 == 0 on step j, then <br>
on the next step j-1 you want f1 = 0 .<br>
Instead by IEEE default, <br>
on this step j you get f1 = -(0/0)*inf<br>
or NaN which pollutes all subsequent f1's.<br>
*/<br>
<br>
f1 = 0;<br>
f = a[N];<br>
<br>
for (j = N - 1; j >= 0; j--) {<br>
d = x + f;<br>
q = b[j] / d;<br>
f = a[j] + q;<br>
d1 = 1.0 + f1;<br>
f1 = -(d1 / d) * q; /* CRITICAL STEP */<br>
}<br>
<br>
*pf = f;<br>
*pf1 = f1;<br>
return;<br>
}<br>
<br>
Thus IEEE default results for invalid exceptions<br>
are unhelpful for this particular example, producing NaN's for the derivative<br>
f1. The key point is that the correct result depends on what went on before<br>
in a complicated way. The straightforward way to code the needed test in<br>
the absence of language support looks like<br>
<br>
for (j = N - 1; j >= 0; j--) {<br>
f1j2 = f1j1; /* Save this guy for recomputation -<br>
won't be defined until j=N-2 but that's the<br>
earliest it will be needed.. */<br>
f1j1 = f1; /* Intermediate save. */<br>
d = x + f;<br>
q = b[j] / d;<br>
f = a[j] + q;<br>
d1 = 1.0 + f1;<br>
r = d1 / d; /* CRITICAL STEP */<br>
if (r != r) { /* r is NaN from 0/0 or inf/inf */<br>
if (d1 == 0) {<br>
/* f1 = (0/0)*inf so return infinity */<br>
f1 = q;<br>
} else {<br>
/* f1 = (inf/inf)*0 so return limit value */<br>
f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];<br>
}<br>
}<br>
else { /* r is a normal number */<br>
f1 = -r * q;<br>
} <br>
}<br>
<br>
<br>
The "normal" case is slowed down by the necessity of retaining f1j2,<br>
the value of f1 from two iterations back, which can't be avoided, <br>
and by having to test<br>
the result of d1/d and branch conditionally. All the exception-handling<br>
mechanism to be discussed below is intended to avoid that conditional<br>
branch or reduce its cost. <br>
<br>
On non-IEEE systems on which division by zero was fatal or infinity or NaN could<br>
not be distinguished from finite values, producing robust results for this<br>
computation would be much more expensive. As they also would be on IEEE systems<br>
which stubbornly optimized away if (r!=r), which could also be expressed<br>
as isnan(r) or unordered(r,r) at somewhat higher cost. <br>
Note that there is nothing<br>
truly exception or abnormal about the case d == 0; all that's required is that<br>
the value of the input parameter x happen to be the negative of one of the<br>
partial function values f, which is mostly a matter of luck.<br>
<br>
==========<br>
<br>
From dgh Tue Oct 8 15:36:13 1991<br>
To: nceg@cray.com<br>
Subject: Exception Handling II: Presubstitution<br>
<br>
<br>
The previous message demonstrated a loop to evaluate a continued fraction<br>
and its derivative, and showed how points where the formula for the<br>
derivative breaks down can be bypassed by conventional means, at least in<br>
IEEE arithmetic.<br>
<br>
The next step is to consider Kahan's proposal for presubstitution, which turns<br>
the loop in question into the form:<br>
<br>
extern double infinity(void);<br>
extern void presubstitute(char *, double);<br>
<br>
void<br>
continued_fraction(N, a, b, x, pf, pf1)<br>
int N;<br>
double *a, *b, x, *pf, *pf1;<br>
{<br>
double f, f1, d, d1, q, bj, bj1, p;<br>
int j;<br>
<br>
/*<br>
* Evaluate function f at x by continued fraction a[j] and b[j];<br>
* function value to *pf, derivative to *pf1<br>
*/<br>
<br>
/* Assume<br>
aj finite<br>
bj finite and nonzero<br>
x finite<br>
<br>
Critical step is evaluating f1= -(d1/d)*q<br>
<br>
If d == 0 and d1 != 0 on step j, then<br>
on the next step j-1 you want f1 = p = b(j-1)/b(j) * d1(j),<br>
a result determined by taking limits in step j.<br>
Since f1(j-1) is nominally computed as (inf/inf)*0,<br>
the desired result can be obtained by substituting <br>
inf for inf/inf and p for inf * 0.<br>
<br>
if d == 0 and d1 == 0 on step j, then <br>
on the next step j-1 you want f1 = p = 0 .<br>
Since f1(j) is nominally computed as (0/0)*inf, and<br>
since f1(j-1) is nominally computed as (f1(j)/inf)*0,<br>
the desired result can be obtained by substituting<br>
inf for 0/0 and p (= 0) for inf * 0.<br>
*/<br>
<br>
presubstitute("0/0", infinity());<br>
presubstitute("inf/inf", infinity());<br>
<br>
f1 = 0;<br>
f = a[N];<br>
bj1 = b[N]; <br>
<br>
for (j = N - 1; j >= 0; j--) {<br>
d = x + f;<br>
d1 = 1.0 + f1;<br>
bj=b[j];<br>
q = bj / d;<br>
f1 = -(d1 / d) * q; /* CRITICAL STEP */<br>
f = a[j] + q;<br>
p = bj1 * d1 / bj;<br>
presubstitute("0*inf", p);<br>
bj1 = bj;<br>
}<br>
<br>
*pf = f;<br>
*pf1 = f1;<br>
return;<br>
}<br>
<br>
With j counting downward, bj1 means b[j+1],<br>
the b[j] that was in effect during the<br>
previous iteration.<br>
<br>
Note that there is an error in Kahan's 17 March 1987 draft of<br>
"Presubstitution, and Continued Fractions" -<br>
the presubstitution value should be written<br>
<br>
b * d' / b<br>
j+1 j<br>
<br>
The first subscript appears uniformly in the draft as j-1,<br>
perhaps a relic of an earlier draft in which the loop ran from j=0 to j=N-1.<br>
<br>
In contrast to the "conventional method" proposed previously, which required<br>
allocating two registers to store previous-iteration data and more importantly<br>
a conditional branch to test the result of d1/d on every iteration,<br>
the presubstitution method requires setting up two presubstitutions outside<br>
the loop to cause 0/0 and inf/inf to be replaced by inf, <br>
and another within the loop that causes 0*inf to be replaced by p, which<br>
must be computed on each iteration and inserted into the arithmetic pipe<br>
in case 0*inf arises on the next iteration. <br>
Thus updating the presubstitution value must be<br>
a fast operation, even though actually doing the presubstitution in the<br>
relatively rare event of an exception might be slow.<br>
<br>
There are a couple of things to note about the function presubstitute():<br>
<br>
1) If implemented as hardware fpop codes, presubstitute()<br>
needs to be faster than a<br>
conditional branch in order to be useful. Thus it shouldn't be a synchronizing<br>
operation that stops the pipe; updates to the presubstituted values have<br>
to flow through the pipe with other fpops.<br>
<br>
Furthermore there must be some presubstitution registers, at least<br>
one each for 0/0, inf/inf, 0*inf, other invalid, finite/0, overflow, and<br>
underflow; on a V8 SPARC, for instance, presubstitution registers must be<br>
able to accommodate 128-bit long double data types, so there are already at<br>
least 7 * 4 = 28 32-bit words of presubstitution registers, compared to<br>
the 32 32-bit floating-point data registers used in normal operation;<br>
all must be saved and restored during context switches. I'd like to hear<br>
from hardware designers about the cost of this approach.<br>
<br>
Alternatively presubstitution registers could be implemented<br>
as 5-bit pointers to the 32 floating-point data registers. <br>
Then data registers used for presubstitution values<br>
(4 32-bit registers in this example) would be subtracted from the pool of <br>
conventional<br>
registers available for optimizer allocation. Some compiler cooperation would<br>
be required to understand this - presubstitution would be an operation known<br>
to compilers. Hardware paths would have to be available<br>
to bring presubstitution values from the appropriate data registers to the<br>
destination register when needed. This operation (inserting the presubstituted<br>
value) could be slow because it's rare, unlike setting up the presubstituted <br>
value, which must be fast.<br>
<br>
On vector machines, a vector version of the continued fraction program<br>
might have x a vector parameter rather than a scalar.<br>
Presubstitution registers would have to be the<br>
same size as vector registers. Or else this kind of code simply wouldn't<br>
be vectorized.<br>
<br>
2) Presubstitute() can't be implemented as a system call. That would<br>
take longer than the conditional compare.<br>
<br>
3) Presubstitute() could be implemented simply by storing presubstituted<br>
values into a normal array in user address space memory. Programs planning<br>
to use presubstitution first register their intent at program initialization<br>
by a system call that tells the kernel where to find this array of <br>
presubstitution values in user space. Thus the array can be updated relatively<br>
quickly - it's just a (probably cached) write. <br>
Performing the substitution when an<br>
exception arises requires that the floating-point hardware be trappable<br>
and continuable after traps. The system call that initializes presubstitution<br>
would also initialize the hardware in this case to trap on every invalid<br>
operation. The kernel trap handler would examine each such trapping fpop<br>
and obtain the value to be used from the user-space table. It is not<br>
necessary for this purpose that user-mode code handle the trap, although<br>
that's one possible implementation. Hopefully the user-mode code won't<br>
actually be written by the end user.<br>
<br>
This appears to be most likely way of implementing presubstitution,<br>
should it prove worthwhile.<br>
<br>
==========<br>
<br>
From dgh Tue Oct 8 16:46:50 1991<br>
To: nceg@cray.com<br>
Subject: Exception Handling III: Operator Modifiers<br>
<br>
<br>
In my first posting, I described a loop to evaluate a continued fraction<br>
and its derivative and how it might be made robust in the face of 0/0 or<br>
inf/inf by conventional means costing a conditional branch on every iteration.<br>
<br>
In my second posting, I described how presubstitution might work. In the best<br>
possible implementation, an extra fp multiply, divide, and "presubstitute"<br>
are required on every iteration. The "presubstitute" might be a load to<br>
a hardware presubstitution register or a store to a user-space memory location.<br>
<br>
I now turn to another approach which is also based on conditional branches -<br>
and thus is probably no cheaper than the conventional approach in this<br>
example. The payoff is in other situations where the exceptionalness<br>
can't be so readily determined by doing just one floating-point comparison,<br>
such as when an inexact must be detected.<br>
<br>
The idea is that instead of testing the numerical operands or result of an<br>
operation, or depending on hardware presubstitution registers or recovering<br>
from an asynchronous trap, the code tests the accrued IEEE exception bits.<br>
It's intended to permit an important optimization on IEEE systems that<br>
maintain separate "current exception bits" for the last fpop: if there is<br>
only one fpop that can generate the exception of interest, then only the<br>
current exception bits for that fpop need be tested. Abstractly this is not<br>
inherently more or less synchronizing than a conditional branch on <br>
just-computed floating-point data, although current implementations might<br>
do a better job on the latter.<br>
<br>
Conventionally, if you want to test an fp expression for a particular IEEE<br>
exception such as FE_INVALID, you would do something like<br>
<br>
tinv = get_accrued(FE_INVALID);<br>
clear_accrued(FE_INVALID);<br>
dest = expression;<br>
tinv2 = get_accrued(FE_INVALID);<br>
set_accrued(FE_INVALID, tinv | tinv2);<br>
if (tinv2 == 0) {<br>
/* code for normal case that no exception arises */<br>
} else {<br>
/* code for unusual case that exception arises */<br>
}<br>
<br>
So that's at least four function calls, worst case, resulting best case in<br>
two stores and two loads of the accrued exception register; plus a conditional<br>
branch. On most current<br>
high-performance IEEE systems those loads, stores, and branch will be much more<br>
expensive than the "dest=expression" that one would most typically want to<br>
test.<br>
<br>
The proposal below is for language support for the foregoing construct, with<br>
the intent that it would often be optimized to <br>
<br>
dest = expression;<br>
if (get_current(FE_INVALID) == 0) ... else ...<br>
<br>
and that consequently there would in time be hardware support for either the<br>
five instructions<br>
<br>
"branch conditionally if the previous fpop generated FE_X"<br>
<br>
for each IEEE exception INVALID, OVERFLOW, etc., or else one instruction<br>
<br>
"branch conditionally if the previous fpop generated FE_X|FE_Y|..."<br>
<br>
where any subset of the five IEEE exceptions can be OR'd together<br>
to specify the desired condition; and thus after further time fpop<br>
implementations might send early information to branch units that certain<br>
IEEE exceptions are guaranteed not to arise on this fpop code - a MIPSy idea<br>
much appreciated by Kahan - releasing control flow past the conditional <br>
branch.<br>
<br>
In my experience the most common case, by far, is that only one operation<br>
matters, so my original idea was to attach the branch condition as a modifier<br>
directly to the operation:<br>
<br>
z = x / y ;<br>
<br>
becoming<br>
<br>
z = x / [on FE_INVALID goto label] y ;<br>
<br>
indicating that the fpop(s) implementing / be tested for invalid and when<br>
encountered, control would flow to a local labeled statement. <br>
A goto was chosen with malice aforethought to insure a lightweight<br>
implementation burden: if the exception arises the compiled code is neither <br>
required to have assigned to z nor to have left z untouched.<br>
There is no possible way to return back into the middle of the expression.<br>
<br>
A semantically similar but syntactically quite different approach, inspired by<br>
Stallman, is to attach a cast-like modifier to an expression:<br>
<br>
(__oninvalid__ label) (x/y)<br>
<br>
or maybe even<br>
<br>
(__oninvalid__ label) (z=x/y) <br>
<br>
The distinction between these two<br>
is important on extended-precision-based systems that<br>
typically won't detect overflow or underflow on the divide but only on the<br>
assignment forcing conversion to storage precision.<br>
<br>
Anyway if the tested code is an expression rather than a single operation,<br>
compilers should eventually recognize possibilities for optimization.<br>
<br>
An even more general approach is inspired by C++:<br>
<br>
try { z=x/y ; }<br>
catch (FE_INVALID) { /* code for exceptional case */ } ;<br>
/* followed by code for non-exceptional case <br>
that will also be executed if exceptional case falls through */<br>
<br>
There can be multiple catch clauses for different exceptions.<br>
<br>
I don't think the C++ try/catch mechanism is appropriate for C - too<br>
heavyweight: it allows for user-defined exceptions and requires dynamic<br>
scoping of handlers: if FE_OVERFLOW is not caught here, is should be caught<br>
by any dynamically enclosing catch even in some other procedure. This is<br>
too much like a long jump - costly and not helpful for what I want to do,<br>
which is very local.<br>
<br>
Worse, the main purported advantage - <br>
a cleaner syntax devoid of gratuitous gotos and labels -<br>
will frequently prove illusory, since it will frequently be the case<br>
that exceptional-case code should not fall through to normal-case code as<br>
it would appear to do in the example above.<br>
<br>
Thus a truly satisfactory syntax awaits discovery and promulgation.<br>
I'm eager for good suggestions.<br>
<br>
I'll close by illustrating the foregoing with the continued fraction example<br>
from the previous postings. First with the operator or expression<br>
modifier syntax:<br>
<br>
void<br>
continued_fraction(N, a, b, x, pf, pf1)<br>
int N;<br>
double *a, *b, x, *pf, *pf1;<br>
{<br>
double f, f1, d, d1, q, r, f1j1, f1j2;<br>
int j;<br>
<br>
/*<br>
* Evaluate function f at x by continued fraction a[j] and b[j];<br>
* function value to *pf, derivative to *pf1<br>
*/<br>
<br>
/*<br>
* Assume aj finite bj finite and nonzero x finite<br>
* <br>
* Critical step is evaluating f1= -(d1/d)*q<br>
* <br>
* If d == 0 and d1 != 0 on step j, then on the next step j-1 you want<br>
* f1 = p = b(j-1)/b(j) * d1(j), a result determined by taking limits<br>
* in step j. Since f1(j-1) is nominally computed as (inf/inf)*0, the<br>
* desired result can be obtained by substituting inf for inf/inf and<br>
* p for inf * 0.<br>
* <br>
* if d == 0 and d1 == 0 on step j, then on the next step j-1 you want<br>
* f1 = p = 0 . Since f1(j) is nominally computed as (0/0)*inf, and<br>
* since f1(j-1) is nominally computed as (f1(j)/inf)*0, the desired<br>
* result can be obtained by substituting inf for 0/0 and p (= 0) for<br>
* inf * 0.<br>
*/<br>
<br>
f1 = 0;<br>
f = a[N];<br>
<br>
for (j = N - 1; j >= 0; j--) {<br>
f1j2 = f1j1; /* Save this guy for recomputation -<br>
won't be defined until j=N-2 but that's the<br>
earliest it will be needed.. */<br>
f1j1 = f1; /* Intermediate save. */<br>
d = x + f;<br>
q = b[j] / d;<br>
f = a[j] + q;<br>
d1 = 1.0 + f1;<br>
#ifdef operatorsyntax<br>
r = d1 / [on FE_INVALID goto invdiv] d;<br>
#else /* cast expression syntax */<br>
r = (__oninvalid__ invdiv) (d1 / d);<br>
#endif<br>
f1 = -r * q; /* non-exceptional case */<br>
<br>
endofloop: { /* some kind of null statement required <br>
syntactically here */<br>
}<br>
}<br>
<br>
*pf = f;<br>
*pf1 = f1;<br>
return;<br>
<br>
invdiv: /* invalid division 0/0 or inf/inf detected */<br>
<br>
if (d1 == 0) {<br>
/* f1 = (0/0)*inf so return infinity */<br>
f1 = q;<br>
} else {<br>
/* f1 = (inf/inf)*0 so return "p" */<br>
f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];<br>
}<br>
goto endofloop;<br>
}<br>
<br>
<br>
The corresponding try/catch looks a little nicer:<br>
<br>
void<br>
continued_fraction(N, a, b, x, pf, pf1)<br>
int N;<br>
double *a, *b, x, *pf, *pf1;<br>
{<br>
double f, f1, d, d1, q, r, f1j1, f1j2;<br>
int j;<br>
<br>
/*<br>
* Evaluate function f at x by continued fraction a[j] and b[j];<br>
* function value to *pf, derivative to *pf1<br>
*/<br>
<br>
/*<br>
* Assume aj finite bj finite and nonzero x finite<br>
* <br>
* Critical step is evaluating f1= -(d1/d)*q<br>
* <br>
* If d == 0 and d1 != 0 on step j, then on the next step j-1 you want<br>
* f1 = p = b(j-1)/b(j) * d1(j), a result determined by taking limits<br>
* in step j. Since f1(j-1) is nominally computed as (inf/inf)*0, the<br>
* desired result can be obtained by substituting inf for inf/inf and<br>
* p for inf * 0.<br>
* <br>
* if d == 0 and d1 == 0 on step j, then on the next step j-1 you want<br>
* f1 = p = 0 . Since f1(j) is nominally computed as (0/0)*inf, and<br>
* since f1(j-1) is nominally computed as (f1(j)/inf)*0, the desired<br>
* result can be obtained by substituting inf for 0/0 and p (= 0) for<br>
* inf * 0.<br>
*/<br>
<br>
f1 = 0;<br>
f = a[N];<br>
<br>
for (j = N - 1; j >= 0; j--) {<br>
f1j2 = f1j1; /* Save this guy for recomputation -<br>
won't be defined until j=N-2 but that's the<br>
earliest it will be needed.. */<br>
f1j1 = f1; /* Intermediate save. */<br>
d = x + f;<br>
q = b[j] / d;<br>
f = a[j] + q;<br>
d1 = 1.0 + f1;<br>
try { r = d1/d ; }<br>
catch (FE_INVALID) { /* invalid division 0/0 or inf/inf <br>
detected */<br>
if (d1 == 0) { /* f1 = (0/0)*inf so return infinity */<br>
f1 = q;<br>
} else { /* f1 = (inf/inf)*0 so return "p" */<br>
f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];<br>
}<br>
goto endofloop;<br>
};<br>
<br>
f1 = -r * q;<br>
<br>
endofloop: { /* some kind of statement required here */<br>
}<br>
}<br>
<br>
*pf = f;<br>
*pf1 = f1;<br>
return;<br>
}<br>
<br>
By the way, figuring out that the important exceptional cases are<br>
(0/0)*inf and (inf/inf)*0, given the input constraints on x, a, and b,<br>
is an interesting exercise.<br>
<br>
==========<br>
<br>
From dgh Tue Oct 8 17:39:45 1991<br>
To: nceg@cray.com<br>
Subject: Exception Handling IV: try/catch a switch?<br>
<br>
<br>
try/catch is a sort of like a switch. Maybe if it were more so, the<br>
problems with labels I mentioned previously would go away, e.g. a<br>
somewhat different syntax for the example I have been using:<br>
<br>
__eswitch__ (r = d1/d) {<br>
case FE_INVALID: /* invalid division 0/0 or inf/inf detected */<br>
if (d1 == 0) { /* f1 = (0/0)*inf so return infinity */<br>
f1 = q;<br>
} else { /* f1 = (inf/inf)*0 so return "p" */<br>
f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];<br>
}<br>
break; /* instead of goto... */<br>
default: /* unexceptional case */<br>
f1 = -r * q;<br>
break;<br>
}<br>
<br>
Maybe that's an improvement: an __eswitch__ is like a switch, except the<br>
cases are cases of exceptions rather than cases of expression values.<br>
<br>
_______________________________________________<br>
Cfp-interest mailing list<br>
Cfp-interest@oakapple.net<br>
</font></tt><tt><font size="2"><a href="http://mailman.oakapple.net/mailman/listinfo/cfp-interest">http://mailman.oakapple.net/mailman/listinfo/cfp-interest</a></font></tt><tt><font size="2"><br>
<br>
</font></tt><br>
<br>
</body></html>