LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 no. 782-787 CO p . 2 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 9 RK1I L161 — O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/statisticaltheor787park Vn? Report No. UIUCDCS-R-76-787 THE STATISTICAL THEORY OF RELATIVE ERRORS IN FLOATING-POINT COMPUTATION by NSF-0CA-DCR73-07980 A02-000018 i? Douglass Stott Parker, J» March 1976 Report No. UIUCDCS-R-76-787 THE STATISTICAL THEORY OF RELATIVE ERRORS IN FLOATING-POINT COMPUTATION by Douglass Stott Parker, Jr, March 1976 Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 61801 * This work was supported in part by the National Science Foundation under Grant No. US NSF DCR73-07980 A02 and was submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science, March 1976. n TABLE OF CONTENTS Page 1. INTRODUCTION 1 2. RELATIVE ERRORS AS PROBABILITY DENSITIES 3 3. WHY ACCUMULATED ERROR DENSITIES ARE ALMOST NORMAL .... 9 4. PROOF THAT STATISTICAL BOUNDS HOLD FOR ROUNDOFF ERRORS. . 20 5. CORRELATED ERRORS 26 6. THE SNORT SYSTEM AND PROBLEMS OF IMPLEMENTATION 28 7. CONCLUSIONS 61 REFERENCES 62 1. INTRODUCTION For the man on the street, roundoff analysis of floating- point computation in synonymous with keeping bounds on accumulated errors: at any stage in some computation, the computed approximation fl(A) to a partial result A satisfies fl(A) = A(l + 6) where 6 = (f 1 (A) - A)/A is the relative error of the computation [16]. The idea is to find limits on the magnitude of 6. Recently, however, a number of papers have appeared adding probabilistic considerations to the distribution of roundoff errors in single arithmetic operations and of--what almost amounts to the same- errors in representing the reals as floating-point numbers ([3], [6], [14]). They show that relative errors are not uniformly distributed within their bounds, but instead cluster around zero. The error bounds themselves are \/ery pessimistic . In this paper relative errors 6 are viewed as random variables with suitable probability distributions. The results of applying classical probability theory to roundoff analysis in this fashion are: (1) A proof of the generally accepted statement that accumulated errors approach normal distributions after several operations ([13, p. 113], [11, p. 104], [4, p. 306]). (2) An appreciation of the rate of conver- gence to normal distributions and how floating-point addition and corre- lated errors can slow this convergence down. (3) A new method of automated error analysis comparable to (but in most ways inferior to) Interval Analysis, but in ewery way superior to Wilkinson bound analysis We should comment that it is not "unrealistic" to treat rela- tive errors as random variables. To do so is just to follow the basic philosophy of Bayesian statistics : if you don't know the exact value of a variable, but know what its possible values are and what the probability of taking each of those values on is, then you should treat the variable as random with the corresponding probability distribution. Treating roundoff errors in this way makes our work fairly clean, plus lets it be general and rigorous as well. 2. RELATIVE ERRORS AS PROBABILITY DENSITIES In the analysis of Wilkinson [16] the relative error 6. = (fl(A) - A)/A in a computed result A is rarely known exactly- usual ly one says that it lies within known bounds. These bounds are proportional to the machine precision 3 , where 3 is the arithmetic base and t is the number of mantissa digits used in representing floating-point numbers. Computation of bounds is simple: if the floating-point approximations to A and B involve the relative errors 6. and 6 B , respectively, then the floating-point approximation to A-B is fl(A-B) = [A(l + 6 A )][B(1 + 6 B )](1 + 6 R ) - AB(1 + 6 A + 6 B + 6 R + 6 A 6 B + 8^ R where 6 D is the rounding error incurred in the multiplication and is bounded by the machine precision. Since hopefully |6^| « 1, |6g| « 1, we have fl(AB) - AB(1 + 6 A + 6 B + 6 R ) (2.2) to a yery good approximation. Then the bound for the error on fl(AB) is the sum of the bounds for 6«, 6g, and 6 R . Up to this point we have treated each relative error 6 as an unknown element of a known interval. It is not a big step to view 6 as being smeared over this interval, having at each point a probability of occurrence: Lower Bound Upper Bound Thus we equate 6 with an underlying probability density p(6) so that Probability (5 £ x ) ° p(x) dx rX Lower Bound p(x) dx (2.3) Recall that a probability density f is a function on 1R satisfying (i) f(x) £ V x dR (ii) f(x) dx = 1 ■ -00 (2.4) and that the probability distribution F corresponding to the density f is defined by rX F(x) J-c f(t) dt (2.5) The question arises as to what expressions like "A(l + 6«)" mean if we view 6« as a probability density. This is resolved by the following definition and explanation. Def : A smeared number (or fuzzy number ) is a pair N = [A,f] where feF and f is a probability density on 1R. The probability that N takes on a value ^ A«C for any CdR is given by f C Pr(N i A-C) = ' f(t) dt J -oo (2.6) The idea is that every floating-point approximation to a number is a smeared number . Thus fl(A) = A(l + 6.) = [A,(l + OL We clarify this concept with three observations: Observation 1. The constant 1 in expressions of the form (1 + 6) should be thought of as an impulse function A (Dirac delta) at 1. This is natural since 1 should represent a density with all its support at 1 . In the following we will use 1 and A interchangeably when no confusion should result. Observation 2. If 6, and 6~ are independent random variables with corre- sponding densities p, and Pp, then 6 = 6, + 6 corre- sponds to the convolution density p = p, * p ? . Since <5 is the density which gives the probability at each point x that ^-i + <$? = x ' Pr(6 * x) = p^y) p 2 (z) dydz y + z = x •°° fX - z p-|(y) P 2 (z) dydz J -00 J -CD Differentiating with respect to x, P(x) = i-c p, (x - z) pJz) dz = p, * p 9 (x) Note from observations 1 and 2 that if p. is the density corresponding to 6., then fl(A) = A(l + 6 A ) = [A, (1 + 6 A )] = [A, A* p A ] Observation 3. If 6, and 6~ are independent random variables with corresponding densities p, and p ? , then 6 = 0) is tfx, = _J_ e -(x-u) 2 /2o 2 /2tF a so its probability distribution $(x) is (3.5) *(x) = 1 /2? a J rx /. x2 /0 2 e -(t - y) /2a dt (3.6) It is easy to show that the mean and variance of \p are y and a , respectively. 11 The Central Limit Theorem may be stated as follows: given a set {6,, ..., 6*.} of random variables having respective probability 2 densities p., means u., and variances a. (i = 1 , . . . , N) then under "certain loose conditions" the sum variable 6 = 6-, + 6 ? + . . . + 6 N (having density p = p, * p 2 * ... * p N ) converges as N -*• °° to a normal 2 2 2 density with mean u = u-. + vu + • • • + y N and variance a = a-i + a 2 ♦ ... + 4 Proof of the Central Limit Theorem is dependent of course on the "loose conditions" surrounding the summand variables and the type of convergence used. Good surveys may be found in [1] and [12]. We will omit any proof here for roundoff errors because the convergence to a normal distribution, although useful intuitively, is difficult to prove and unnecessary for obtaining statistical bounds (see Section 4). The idea is that the smeared number A(l +6, + 6 ? + ... 6 N ) is equal to A(l + 6) where 6 has a density that is almost normal. Since: 2 1. We can compute fi's mean u and variance a easily. 2. Assuming 6 normal, the probability that |6 - u| = 3a is 0.002798 (i.e., a 99.7% confidence interval on 6's value is [y - 3a, y + 3a]). 3. For N larger than 2 or 3 the error bound [u - 3a, u + 3a] is inside the crude Wilkinsonian error bound [-NB 1_t /2, N6 1-t /2] for Rounding or [-Ne 1_t , 0] for truncation 12 --apparently without much effort we can get much tighter, more realistic error bounds and improve our error estimates. This gives us our motivation: statistical error bounds seem tight. But we need more details. We compute the mean and variance for the rounding error density p R defined by (2.11). For rounding , we have ^Rounding = "Rounding = ^"W/S - ^12 (3.7) For chopping , similarly Chopping " B /2 a 2 - (e^) 2 - (-g 1 ^) 2 = e 2-2t , . Chopping 3 2 p /u ^- i5) Therefore o D d l f a D ,. = a r . . = 3 1-t //T2 (3.9) R Rounding Chopping ' v ' Notice that from (2.9) and (2.10) expressions involving only repeated multipliation, repeated division, or mixed multiplication and division 5L___— — will result in smeared numbers of the form [A, 1 * p R * p R * p R * ... p R ] where the number m of p R 's in the convolution product depends of course on the number of operands in the expression and the accumulated error for each operand. For example if A,B,C are irrational then fl(A(B*C)) ■ [ABC, 1 * p R * (p R * P R * P R ) * P R ] since each number has a representa- tion error bounded by p , and the two multiplication operations have K 13 rounding (truncation) errors bounded by p R . The Wilkinson bounds for the 5 1-t 5 1-t relative error in ABC would be 6 ABC e [--^ 3 , j 3 ] for rounding and ^ABC e ^"^ " ' ^ ^ or truncatl0n ' If we assume (Central Limit Theorem) that (p R * p R * p R * p R * p R ) 2 2 1-t is normal with a = 5a R , y = for rounding or y = 5(-l/23 ) for chopping, then the error in ABC should satisfy f ["Vyf ^ 1 " t »Vrl ^"^ Rounding ^ £[M - 3o ' y + 3o] A(^f-3^-t, (^s/^ Chopping 99.7% of the time. Now 3^-rp- * 1.94 < 5/2, and we have improved our error bounds by 20% after only two operations. As the number of operations N increases we will get better improvement since Wilkinson's bound are proportional to N while the 3a bounds are proportional to •ft. It is not unreasonable to assume that (p R * p R * p R * p R * p R ) is normal, either. Figures 1-4 show the result of repeated convolution with the error density for Rounding. Already (p R * p R * p R * p R ) is essentially normal (Figure 4). On a superficial basis at the least, convergence to a normal density seems rapid. It should be pointed out that the use of division introduces "checked" densities p R in the convolution product (c.f. (2.10)). This is not a problem because p R (x) for Rounding Pp(x) = Pp(-x) ={ -, t p R (x - 6 ) for Chopping (3.10) 14 a = 6 1_t //T2" »t-l -26 1-t 23 1_t -23 1_t Figure 1 . Figure 2. p R for Rounding and its normal approximation Pd*Pr an d normal approximation o=B 1 - t /2 ■2B" " 26 Figure 3. p R *p R *p R and normal approximation 1-t a = ^' X //3 Figure 4. p*p*p*p R and normal approximation 15 Obviously what was said in the previous paragraph holds without modi- fication for Rounding. For chopping we must simply be more careful in v 1 1-t computation of the mean u, since the mean of p R is + «-3 instead of - iS as it was before. For Chopping, division errors tend to cancel out multiplication errors since they are of the opposite sign--this is reflected here by the fact that division makes the mean of the final error more positive, while multiplication makes it more negative. We commented in Section 2 that addition and subtraction were not as nice as multiplication and division insofar as their results were not easily expressible as smeared numbers. For the rest of this section we address the problem of reducing a sum of smeared numbers to a single smeared number. It is expedient to work in absolute errors instead of the relative errors we have been working with until now. The following definition and lemma provide the necessary groundwork. Def : For a bounded density / and nonzero real number A, the stretched density / is given by / A (x) ■ —/(f) (3.1D |A| A A We say that / is / stretched by the value A. Here we have restricted our definition to bounded densities to avoid the problem introduced by impulse functions. Lemma 1 If / is a bounded density, then [A,/] = [1,/ ]. 16 Proof Since f A-C rA-C 1 r(t)dt = ] ^/(|)dt| y=t/A = (_ /(U)dy. [1,/ ] is equivalent to [A,/] under the definition of smeared numbers. Notice that [1,/ ] is an "absolute" smeared number, no longer relative to A. Effectively [l,/ A ] = / A . Theorem 2 If fl (A + B) = [A, 1 * p ] + [B, 1 * p £ ] and A f 0, B/0, A + B ,* then f 1(A + B) = [A + B, 1 * p] where p(x) = p{ A/A+B ) * p< B/A+B) (x) (3.12) Proof Brute force. Note that (1 */)(x) =/(x - 1) and (1 */) (x) *TK[flj - *>■ Now: fl(A + B) = [A, 1 * p,] + [B, 1 * p 9 ] = [i, (i * Pl n + [i, (i * p 2 n = (i * p/ *.(i * p 2 ) b (Lemma 1 (This follows from the remark that [l,/] = /, and observation 2, Part 2.) t=(A+B)y J _oo II I I A + B AB p 2 ( (A * B ^ " B ) % y=v+(B/A+B) A + B AB r°° , x - (A + B)v - (A + B) H* A ) /(A + B)vx . P 9 (- i — 5 — H dv 17 A + B L D (B/ A AB B A + B 'A+B), A + B Pp v v p (A/A+B) ( A + B 1 - v) 1 JA/A+B) * n (B/A+B), x n hTBTI p 2 l A + B " u A + B H 1 i * D (A/A+B) * (B/A+B), x , |A + B| p l p 2 l A + B J (1 * p (A/A+B) * p (B/A+B) } A+B [1. (1 *p| A / A+B ) * p (B/ A + B) )A+B] = [A + B. 1 * p{ A/A+B) * 4 B/A+B) ] Notice that A - B = A + (-B) so Theorem 2 applies to sub- traction also. Theorem 2 leaves us a few steps short of finishing off all the problems of addition. Corollary 1 . n If fl(_Z A.) = [A r 1 * p 1 l + [A 2 , 1 * p 2 ] + n n then fl( J A^ = 11 A r 1 * p] . + [A , 1 * p ] n r n where 1=1 P = P i=l (A^ZA.) (A /^A.) T i 2' i (A /ZA.) n i (3.13) ( Proof . Trivial. Naturally we now require that A. f (i = 1 , . . . , n) n and I A. f 0.) i=l n 18 2 Theorem 3 . Let the mean and variance for the density / be y and a , a 2 2 respectively. Then the mean and variance for / are (Ay) and (A a ) , respectively. ( Proof . Trivial.) From coral lary 1 and Theorem 3 it is immediately obvious how to apply the central limit theorem to addition and subtraction. 2 If the relative error density p, of A has mean y, and variance a 1 2 and the density p~ of B has mean y~ and variance a ? , then the density p of A + B (as in Theorem 1) has A B mean y = ( A + B )y ] + ( A + B ) y 2 (3.14) 2 , A \2 2 , B ,2 2 and variance a ( /TTb) a l + ( FT~B ) a 2 (3J5) After many additions the error density will converge to a normal probability density—but more operations are needed because the stretched (A./ZA.) densities p. J "" will not be identical up to translation as they were for multiplication: accumulated errors in multiplication become normally distributed much more quickly than those in addition . The reader who has researched the Central Limit Theorem has probably dis- covered that most of the existing proofs assume identically distributed summand variables. Unidentical distributions can slow down convergence-- especially when the variables differ greatly. Thus when stretching factors like A/(A - B) became large (implying lots of cancellation and loss of significance in the computation of A - B), convergence toward normal error densities is slowed. In the next section we show that we can guarantee confidence intervals almost as good as [y - 3a, u + 3a] all of the time, even with 19 slowed convergence. Thus convergence and the Central Limit Theorem in general should be used only for an intuitive feel of what is going on, while the growth of a is worth serious attention. Note that when stretching factors like A/(A - B) get large, convergence not only slows down but variances get large also (Theorem 3). This is a serious weakness of statistical error analysis. 20 4. PROOF THAT STATISTICAL BOUNDS HOLD FOR ROUNDOFF ERRORS Section 3 introduced the Central Limit Theorem, showed that it could be used to show how accumulated roundoff errors have asymptotically normal distributions, and discussed the rate of convergence to normal distributions. The literature surrounding the Central Limit Theorem is labyrinthine since there are any number of assumptions one can make about the summand variables which simplify or aggravate the proof. How- ever, as we mentioned earlier we are not so much concerned with proving convergence as being able to say, for any roundoff variable 6, some- thing like Pr( |6 - y| * 3a) ^ .0028 (4.1) as we can with normally distributed random variables. In general given a fixed small probability p we want to find a corresponding X such that Pr(|5 - u| = Xa) - P. (4.2) Since 6 will correspond to a lengthy convolution of error densities, determining X exactly will be unfeasible. Instead we can get bounds on A. It would be simpler to just assume that 6 is normally distributed and look up A in normal distribution tables. However, this is not always legitimate: see e.g. [5]. Although three or four convolutions can approximate a normal density well (see Figures 1-4) large stretching factors can upset this approximation. The purpose of this section is to bound how much these factors can upset it, by obtaining bounds for A in (4.2). We start with Chebyshev's Inequality . 21 Theorem 4 . Let 6 be a random variable with underlying density p(t), 2 having mean and variance a . Then Pr( 1 6 | = Aa) = 1/A 2 Proof , a 2 = E[6 2 ] J_c r p(t) dt => L> t2p(t) dt P (t) dt '|t|=Ac (Aa) 2 |t|=Aa = (Aa) 2 Pr(|6| = Aa) Pr ( | 6 1 = Aa) = -j , as desired, A This says that if we want Pr(|6| = Aa) = .0028 then we will be guaranteed this relationship provided 1 0028 =s> A = 18.89 Clearly this is not very good since A = 3 works for normally distributed variables, and after many convolutions we expect 6 to be very close to normal. It turns out we can obtain a much tighter limit on A if we go through the analysis used to obtain the Chernoff bound . N Theorem 5 . Let 6 = I 6. be a roundoff error, i.e. s <5 is represented by i=l n 22 the convolution of N (possibly stretched) copies of p R , the rounding density. Then Pr(|6| = Xo) = 2e 2 where a is the variance of 6. Proof . Set d(x) = {° J f *° Then d(x) = e(x) — exp(L(x-Xa)) for any positive L, This implies Pr(6 = \o) = E[d(6)] = E[e(6)] = E[exp(L6) exp(-LAa)] = E[exp(L6)] e" LAa Now if p(t) is the density associated with 6 then E[exp(L6)] = exp(Lt) p(t) dt J -00 00 f] fCO I — j- t p(t) dt [p is of compact support] n=0 n - J-" oo m L n = y J] n=0 where M is the nth moment of p(t). Since p(t) is symmetric about (i.e., even) the odd terms vanish and we are left with 23 oo M L 2n E[exp(LS)] = I 2" n=0 K '' We now need the following lemma. Lemma 2 . For roundoff densities p(t) M 2n < (2n)l (M 2 ) n n!2 n Proof . It can be shown that if p(t) = p, * p~ * ... * P N (t) and M . is the nth moment of p. , that the nth moment M of p is k n n , N M n = I n 'n ' n ' , n M n k n Q n r n 2 . ... n N . k=] n k ,K n 1 +...+n N =n This is done as in Theorem 1 by comparing the characteristic functions (essentially Fourier transforms) of both sides of the equation N p = II * p . We then use the fact that the p, are stretched versions k=l k k of the roundoff density p R (2.11), with moments [A k (e 1_t /2)] n n,k n even n odd (A ) where p. = p D , i.e., A. is a nonzero stretching factor. Pulling all k K k of this material together leads to the bound of the lemma. It is interesting that if p(t) is normal with mean 0, 24 (M 2 ) n a 2 " n! 2 n (M 2 = ° ' --so we have equality with the stated bound. Therefore: 1. The bound is fairly tight, especially as N gets large. 2. The lemma holds for other densities besides stretched versions of p R . 3. There is probably an elegant proof of this lemma. It now follows that oo L M 2 2 E[exp(L6)] = I ] -r (-- T ^) n = exp(^-f-) n=0 n ' c c and therefore > < L 2 2 P(6 = Xo) = exp( — y — LAa) for all positive L. The Chernoff Bound is now obtained by choosing L such that the right hand side of this expression is optimal. By differentiating it is easy to see that L = A/a is the correct choice. It now follows immediately that 2 Pr(5 = Xo) = e" A /2 and 2 Pr(|6| = Xo) = 2e" X /2 as stated. The Chernoff Bound on the probability is much tighter than the one obtained by Chebyshev's Inequality, and is almost as good as having a normal distribution. For comparison we tabulate the A's obtained from each method: 25 Chernoff Bound for Chebyshev Bound for > X X X Pr(|6| = Xo) (6 is normal) (6 is roundoff error) (6 arbitrary) 2.448 3.162 2.716 4.472 3.035 7.072 3.255 10.00 3.625 18.90 3.717 22.36 3.899 31.62 4.292 70.71 4.799 223.6 5.257 707.2 10 1.645 05 1.960 02 2.326 01 2.576 002798 3 002 3.090 001 3.290 0002 3.719 00002 4.265 000002 4.753 Table 1. Bounds for X given information about 6 Table 1 shows us that 99.7% confidence intervals on the magnitude of 6 are (-3.625a, 3.625a); this result compares quite closely with the (-3a, + 3a) bounds we could use if it were known that 6 were normal, and has been used in the program that actually does statistical error analysis (Section 6) . 26 5. CORRELATED ERRORS Until now we have been assuming that all our relative errors were independent random variables. Now we consider the case in which they are not all independent, i.e., some of the errors are related to others in a common expression. This is the case, for example, in the evaluation of fl(A-A). Assuming that fl(A) = A(l + 6 A ) = [A,l * p A ] (5.1) then our prior analysis would claim fl(A-A) = A-A(l + 5 A + 6 A + 6 R ) = [A-A.1 * p A * p A * p R ]. (5.2) This is not correct, however; it is too optimistic. The correct ex- pression is fl(A-A) = A-A(l + 26 A + 6 R ) = [A-A.A * p^ 2) * p R ] (5.3) (2) where pi ' is p A stretched by the factor 2, as in (3.11); i.e., p A 2)(t) = lP A (t/2) - (5 - 4) We remark that the rounding error 6 R (with density p R ) is always initially an independent random variable . Correlation comes from repeated use of an existing error as in the case with 6. above; new rounding errors are always uncorrelated with anything. 27 The generalization of the above analysis to the general case is straightforward. Consider a smeared number A (c-,) (c ? ) (c N ) A = [A,l * P] * p 2 * ... * P N N ] (5.5) where p. is a probability density (k = 1, ..., N) and the C. are stretching factors. Suppose that two distinct error densities p., p. correspond to the same error variable 6. Without loss of generality (convolution is commutative) we can assume these densities are p, and p ? ; A can then be rewritten as (Ci + < A = [A,l * p 1 ' * p 3 J * ... * P N n ] (5.6) (c, + c ) (c ) (c N ) This corresponds to merging the two references to 6 into one. Other references to 6 should also be merged. Note that, because convolution of densities (or addition of random variables) is commutative, correlated errors must be merged throughout computation. For example, in computing X = fl(---(((A + B) + C) + D) + E) + ... +Z) + A) (5.7) the error 6. involved in the last step of the addition will be slightly correlated with the error of the rest of the sum; the resulting smeared number is X - [X, 1 * P< 2A/X) * P< B/X) * ... *P< Z/X) ] (5.8) where X=2A+B+C+...+Z. 28 6. THE SNORT SYSTEM AND PROBLEMS OF IMPLEMENTATION As this project of statistical error analysis progressed, it was hoped that the final results could be automated to perform a_ posteriori error bound calculation in the spirit of Interval Analysis. Since for small relative errors 6 = (fl(x) - x)/x we have x z fl(x) (1 - 6) it seems natural to evaluate fl(x) in the natural way and use statisti- cal bounds on 6 to bound the error in fl(x). While Interval Analysis will give results like x e [x low' x high ] we would give results like x e [fl(x) (1 + u - B), fl(x) (1 + u + B)] 2 where B = (3.625a) is the 99.7%-confidence bound on 6, u and a being the computed mean and variance of 6 viewed as a random variable. Automation of the statistical analysis was in fact completed by development of SNORTRAN, a PL/I-like language. It is similar in philosophy to PL/I-FORMAC in that SNORTRAN source is preprocessed and expanced into PL/I, which is then handled by the PL/I compiler. The preprocessor, a SNOBOL routine included at the end of this section, is unsophisticated: current need did not warrant sophistication. The results of several test programs (also included at the end of this section) reflect negatively on continued use of statistical 29 relative error analysis for almost all problems . There are at least four reasons for the discouragement: 1. There is no way to represent as a smeared number. 2. In subtraction of two numbers with similar values, can- cellation causes enormous increase in error variance (stretching factors get large). One or two such can- cellations will result in a value of a so large that "useless" error bounds result. 3. In only the most well-conditioned, stable computations (e.g., repeated multiplication) are statistical bounds competitive with Interval Analysis. 4. Correlation between errors cannot be easily taken into account (as shown in Section 5) so statistical "bounds" ignoring correlation may use smaller stretching factors than they should, and give incorrectly small results. The actual implementation uses Interval Analysis methods to generate worst-case values of a and y, hence statistical bounds gen- erated by the program should be bounds (ignoring correlation effects) There are four basic SNORT commands, all of them keyword-driven for simplicity: TMTTTni T7r j ROUNDED) INITIALIZE {cHOPPEDJ' SMEAR [,]; EVAL = ; PRINT [,]; 30 INITIALIZE causes the preprocessor to copy SNORT system declarations and routines into the translated program output, depending on the type of arithmetic being simulated. If rounded arithmetic is being used, the preprocessor copies in the file SNORT. RND, otherwise copying SNORT. CHP. The former is reproduced at the end of this section; the latter, requiring full Interval Analysis to compute error means, was never written. SMEAR results in declarations for which permit it to be used as a smeared number. "" may be sub- scripted if desired (e.g., SMEAR A (10,0:5)) and a list of s generates a list of appropriate declarations. EVAL = causes an SLR(1 ) parser to break down into SNORT system subroutine calls. Currently should only involve the operators +,-,/,* although this could easily be expanded or modified by changing the parser tables. All variables in and should be SMEARed--the pre- processor does not currently check. PRINT provides formatted output of smeared variables. It lists computed values, 3.625a bounds and the value bounds they imply, and Interval Analysis value bounds. Four sample runs are included in the following pages, in- volving algorithms which highlight progressively more serious problems with statistical relative error analysis. The first shows that even in the best case of serial multiplication, statistical analysis does not give better bounds than interval analysis. The other three show 31 respectively what happens when addition, subtraction, and computed zeroes are added. Preliminary example: Computation of sin (1.0) The first test of the SNORTRAN system attempted to compute sin(l) = .84147098... via the little-used relationship 9 n=l 2 n This concentrates on the roundoff properties of multiplication (ignoring computation of the cosines) which we know statistical analysis handles well. Unfortunately the results are not encouraging. Statisti- cal analysis obtains the bounds sin(l) e (.84146956, .84147235) after 32 factors, while Interval Analysis gets sin(l) e (.8414694, .8414710). 32 » < LU CD ** « H * * * # ♦ ». t 1 Z t 1 1 •• • • • 1 — CC — • Of • 1 # •o c • 1 » »- -4 QC # If. | V. or CO vO | Ul < X LLI a: «n 1 O LL 4~ □ 1 O 1 Cl z t- in 1 O CO a. X a O 1 O O 1-4 u — < -0 I O O I i- LL 0> | O CO — < 00 | O # z CD t- LU i- 1 O c V 2: z O 1 O 1- mi x UJ 1-4 1 O *- t- ■~- 1/1 l/"> CO 1 O 03 < i/, LU -* 1 ^ O | _j c or CO I at O < LU o a. <* 1 0T UJ (M 1 H O a. * » — a i-O 1 uoo or LU sT r- 1 < O I Lli X •^ a. ^ ! 9W a. z a K X z LU - 1 tt |o « at — •_ . H >* 1 CO O CO o O CO > LU CO 1 1— # LU 2 o < < • 1 ON II or on LJ X O O I 1 < LL O _ — *- mm r> — IO Uj 1 ILU ►- CO *- Z IN < K- co 1 CO O Z O cc *— »» "^ < *- 1 J t < — d < 1 l cc LU •— 1 1 K X) II •— 1 < 1 ~H CO CO < a a CO t- 1 l< II 1 c OC Ul X Ul 00 1 LL IZ -J -1 .« LL 1— *V z LU LU ( O a. -J < O O •» u LU < — *— — a 1 1 O O < > Z IX a ISI ~ X CO a: OT | z IO LU LU Z — r> — -J mm a LU c 1 »— " z — -1 o -J < Z CO 1 1/1 1— • CO Ui LU < > •— 1 ( a. y — UJ CO 1 1 1 or >- a. O K 1 1 t 1 LU UL Z a lu ♦ I UJ LU »• O IA h» 0> CO -4 *\ ITl -4 O -^ l/l r-4 00 -1 1 ~o ift LU m m O m (m — 1 o» h. r- a> >t >»■ r» * * LU a a 0000000000 00 000 000 000000000000 000000 -^rgro^ino^aocO-H icgro^irtot^oo 000000000" -1 "*•""• m+ ~+ r+ r* I -M —< 000000000000 000 000 o o * I Ui LU <\J CO o o f- •»• 00 ift rO 0> -J -1 0* -41 o in >*■ o I ooim 1 LU -O >0 LU ao ia a> o -o c> -o o ft o> «* «»• (Ni ^ »-l ^ o r» >r <* ao • 1 t O f\l co as ut I O o r» >*• — — — • or CO o or or LU •• "UJQ UJ > Z O — 3 -1 »- O < < CO > _l UJ UJ IO OKUQ uu z. z. t- Z IU 3 0000 a. » ca O O Z CO O Z O ►- Z3 O CO O I >■ co »* -J < • z x o» < O 0" -• -J --. » LU > aj u> —1 ac I C\j _J UJ < v0 CL (— m m m> 3 O <1J u S- Z3 O CO or: I- or: o sz co CD Q. E (O X a. z' 33 Example 0: Computation of exp(8.0) This program involves the use of addition, multiplication, and division by computing exp(8) with the usual power series, truncated o at 40 terms. To ten digits e = 2980.957987, and at the end statistical analysis claims e 8 e (2980.9549..., 2980.961...) while Interval Analysis gets e 8 e (2980.955 , 2980.958) - a slightly tighter bound 34 * CO ♦ < Q or •a- I Ul X r» O in >o o o» r» in s — i r« ■► 0" O uj COS O >C in h- 0> 0* coo ^ a co -COO • • • — i rvj <\J o o «1 CC r>- -* C "\ O IA o> -« OM r^ op O O in o ♦ r» o ui C l/N in O O CO o o f- CO CO • • • -H «M (M o o ♦ ♦ LU Uj m —j -1 O O f»- m -* •3 ^ O O M .* + s O uj in ~< a cc o flD ♦ 0- Ouj -O — M CO r- -r m moo >r o o n CO CO r» O O • • • _. (M IM +-> O- 4-> 3 o fl- its CD a. t- X o Ul ,, z < u. • » UJ — 1/1 •^ or «» t- z z UJ X < ^ ♦ X o a 3 z of X -J I/O z + o UJ h- o X ex z LU z < X < U- EL X —• oo • * t~ # LU z * O ir, z II — • » » ^* • * >* z X II •• < o z o • • c z X X uj 1- LU o c o M 1— X o- — c 1_> .. o UJ H K o a X «/> z < o t o LU z < X LU Z 3 u. UJ ~> • •— i uo X u. UJ o o • • o -J n *~ m or » o • II z _l _J _J —I Z »« >— z LU *-< II _J <: < > > CL Z O i^J II 1- X D U LU LU LU Q- LU •— • » II U a. ■ • UJ _J X z > > > z O — V) LU LU LU UJ UJ of a. • • o X OOOOOOOOOOOOOOO QOOOOOOOOOOOOOO OOOOOOOOO-^^*-^"-'!-''^ OOOOOOOOOOOOOOO -J fO o o •♦■ J> Ul U' o ^- cm rn o m r>- vo s >f m MOO * f- + •> O LL' u m in 0> in in r- o* J> o- o o n nJ- co co Ofl CO 1 + • • • LU -H IM (NJ C I c lA CM CO r<- « o or o o oo or m or a> uj •• O 30 ui c o> > z • — 3 nj k O o c c _l — CC < 1/ u > r z «/i z c •— C D l_J i>1 Uj O I >- ►- CC »" — I 13 r~ < o. < • z a: i c> «a c l; o» o — -I wr. in <-i or (M — I ui •o a i- • i z -« rri O O 4- + LU Ul (*1 -^ f^ CO O O O CT> o> o o ~r m co m o in cvj + SO UJ <*• >n in m in in ino 1 0~ eo O o <*\ r~ oo "o o vo a- o^ ■♦■ • • • UJ -^ IN CSJ O I o in nj CO r>- .. O or o o CO Q. m or o> lu •• O w 00 UJ O o» > z • — 3 c\l ►- C «* CO -J •• UJ UJ uo a <_> O z z IJJ 7 L' D root? _J — CC < \r. u > Q -T t/» z o — OZ.i_.ijn lu C I >- »- cc »• -J _> r«- »i a < • z i 5. o- < C O U O — _J uo C < * Ul > m •— o: rvj —t u. vO a t- • X Z ■— J «"| o o 4- + UJ U) rO 0> ro 0^ O o> mo INJ N- ♦ C 0> Ul in >sr m co m m -coo c^ O CO CO O r~ o» j> + . t • lu ^ i\: LU o OE 5> IjO lu C > Z • — 3 M-G < CO _J •• 111 Ui uo or o o •• z ^ Ul 2 Uj Z> ■D C CC -> —X «3 l/l u. > r o B u> CT l- X n. <3 > X o c l_) ►- -> Z u- o — l_) I > M _l r- < t z a- < a _i i/1 Q < # Ul > m — or nj _i ui a- ♦ • t • Ul -* CN4 ui •• O I/) no ui O a > z • — z> Nl-O < CO _J » uj uj tr ■X U CZ _5 O C O _J — cf < r?w O 3 O i/) UD I > I- C K _l D r- O — -I «/. 3 < • ui > in — a- IN -J Ul O O. K- • X z r^i — — C lO (U o 3 o oo < o 00 CO D. X O a> Q. E IT3 X K a x Q. X 35 Example 1: Computation of cos(4tt) This problem involves all four arithmetic operations (using again the usual power series) and is not well-disposed to statistical analysis since there is eventually tremendous cancellation between terms. Note that to the first ten terms the computed sum is roughly -309, and the eleventh term is around +387.6, resulting in a factor of 16 increase in variance. Statistical analysis can no longer give decent results: it provides the bounds cos(4tt) e (-.4183606..., 2.4383515...) whereas Interval Analysis gets cos(4tt) e (.7754054 , 1.227057). Neither of these results are good but the statistical bound is terrible. The relative error bound is so big that the assumption that second order effects are negligible is in question. 36 « < < < <5 < < < ►- >- I- >- k- t- >- 1— UJ UJ UJ LU Ul UJ UJ LU >»» is CO CO CO CO X X X >v » *"* • '» » # • * * * # at. O -^ INI m -J- J- r^ -J- ~J- IM m < O O "5 O 3 3 O O 3 3 3 3 O a X ♦ * * + ♦ + ♦ ♦ ♦ ♦ + ♦ ♦ ♦ J> 4 »— UJ LU UJ UJ LU LU LU LU UJ UJ LU UJ UJ UJ UJ LU in ao m 0- -J- ru 1^1 AJ in X O LU >• r- X O -0 CM O in m m — < 3 r- •e _ AJ ^ o I < 3 i*> r~ -»■ ro -J- X r^ O —* AJ 3 cr O 3 r- X 3 —4 X JJ AJ O o — z -4 r^ 4 n + f-j 5< + to 3 * rj» CC ■*■ in J- X J- A X ^J *1 x -0 LU rr> -c UJ 3 UJ O UJ -0 O UJ ■0 Al UJ in -A UJ X *> LU ' * —t %l ^vJ "A f^ r^ 3 — < X O -0 3 *^. -A Al * X UA A- n n m JJ * O — * 9- ro t> J- m O n 3 *r -J- in 3 lT 3 cn 9 CT -O in Al — 2: < •r at O •T ■O O -i- r^ ^ -3 *H -* ,0 -r •*■ 3 _. UJ LU T\ n in O ^ X X O •O wa\ O O INI 0- CT- -4 M r~ 3 3 n O S> * -4 ~4 3 AJ AJ n X O — AJ i/> > t » r» f~ O CI ^J■ in n 3 3 3 0" O O n O O n —4 p4 ♦ -0 LU 3 Z — < .£ a I r- f~ 1- -* a- LT- ~ ' -T vT AJ Nri -H AJ m ~* r — ' ■* r^ — ■ — * -* -0 z 1 1 1 1 1 1 1 1 ■ » UJ LO • ► — H- >- «■" Z 3. -J -1 1— • «■ ♦ a. -< * O X v}" z < O -^ —* tM p4 ^ -* -T -* -r -^ * —« J" AJ «"l uO O -J m z u_ + ♦ * O + * * 3 ♦ 3 + 3 + 3 * 3 + 3 3 + 3 + 3 X • ■ :n II 1/) •^ L/A ■>> UJ UJ LU LU UJ LU Ul LU UJ IU UJ LU LU UJ IU XI a. o Z z in —* ^ CO 0> - 3 AJ X AJ O B 1^ O >}■ •r c in O- Al * 0" CJ * ^ ^ "v Z ^ r^ r«- ~D X 3 AJ r^- M -* -r I-* 0> X m ■x. z. • » l/l -2 + in O O r^ O X m C7- ? ■T LA 3 n m O O a uo N ^, ,» ►— r* -- — * - r\J O X r*. r- TD O CT -0 3 — * -O CiJ > -J O p ^ .-A 3 r-« ^~ JJ -j 1 * • —i T + n "3 + >NI -0 + O AJ + -> n 4. -A r«1 X st * ♦ AJ x * x M 3 > «-_" X < uO -T UJ -*> •JJ O O XI *0 T- LU ^ri LU JA LU in 3 UJ X CT X) <_: •~~l Cij < # UL u 3 x> n f\; —1 -0 P^ X -r — < m in CO J JJ "A rA X X uj in JA CJ> M K ►- z 0- r> O) T ■n IM 3 3 T T m O 3 3 X 1 — "* i* t> AJ co m <»• ■O ■0 •r r^ f^ X — 1 -H O J- ~r Nl O O JA _< ~-4 r* r- O z U — O i/i X II in ^ ■n 3 O *- -0 •n J- —4 —* O ~r ■T •n AJ AJ 0> O ■■4 0> c» m . • < tm vf z 1 II • • »4 N O ^ r\j •^ -- m p^- 3 vr in 0> a- -T -r -j -* * 3 AJ AJ >» n X X m X ■^ —* < □ LL O ,„ z X O r^ r~ '3 H •O O -J- n 3 3 3 CT O Q 3 LA ■O O n ~4 — 1 O J- •0 X aJ p» LU O a II >— X l/l + UJ r» • r» ♦ LU ^ a- O LU _I •T -r ♦ • UJ ^J J J + • UJ AJ J J + LU * J, J, + UJ A. z _! ♦ UJ — < x M a ■i l*\ tm O LU r- 1— on □ tn 1 1 1 O 1 O 1 1 1 1 O 1 1 1 3 1 3 | 1 1 3 1 1/1 z z JO O » O 'JJ z < a O f\l O O 3 3 3 3 9 z d X 4 ■* LLi ^-< t > 7> OL z ■b n "j J^ ■3 O) CI O -0 l3 j N II II t~ X u LU UJ LU CL JU •r Qt r- X -r -x ^ U. 1*1 JC 1— j: 3 LL r~ JL a. II L_J un •» m ^ r^ X cr -r or 3 -r: aj or ■— a t X O LU r- UJ LU -r LU jj UJ in UJ u jj -J V ~>J -? •r oo AJ 1/1 O l^> J> LO ac 4 X X X u_ X O* JJ O _< UJ •u 3 -r UJ J -* IU c> AJ IU 3 X ■XI 3 ■^ LU 3 3 p^ a LU N. >• 7. -0 -> £ ■n > ;» 3 > *r z <0 > 2 —* > J» > z 3 3 3 3 3 3 3 3 h- < _l —J _J _j _i r- K- O 3> t- ~J -r 1— a ^ >_ J p^ »— -1 prf h" O _, H J 1- a UJ Prt (.j -1 < < < < O 1 < -n < 1 •-[ < 1 < "A 4 O 1 -> > ^» > z _/ • • _ f • • —1 •• _! •• — 1 •• —1 • • «j _i •• o — • y> LJ LU LU LU u LU LU or aj I/) c 'JJ LU '-> Li IU of UJ Ll u0 UJ or IU L-1 3 UJ UJ u 3 XI or UJ UJ UJ O ^1 c UJ LU 1/1 a St #t Z ^ •• <_' ^ •• ir z. • • 3 ^ •• Z z • • .-: z • • it: •• jr z Q. UJ Z Ul ■D u ? U.. ~ LU ^ i'-' 3 ' UJ ■» ■ u 3 IU L7 lU ~> 'U » U ' ^ 'XI r* x. _i LU "^ 'U 3 3 O O ■3 3 rj r3 ■J 3 ■D 3 ■i a O 3 a a 3 3 3 3 ,3 i_i _A O lJ 3 -J — i ^ — s _ 1 — X _J •«* CO ^ — • an ^ _ c_ J — n _l — ^ < •> u. < 1^ u_ •a 1/1 LL 4 \s m. < 1^1 u_ 4 i/> u. < ^) IL 4 1" l^ — * > n T ^ ^» 7T L/l > ^ \n > O — m > 3 7" I/ 1 ■> 3 j^ y > *A 7" L0 > O T in «" ^ *-> •^ ?• LJ ^-1 2 c MM s. 5 -* S. O !-• •• J •— \ LU ■« z 3 •— IU O 3 O ^ O rj J ^) ^? a u 1/1 .3 3 _) ./> UJ 3 u,' lO »J 3 L-, 3 z. J l/< CJ JJJI JJ C 1 > IU a J_ >- U>> ■3 1 1» UJ O 1 v 1U T 1 > UJ C) l_ >- LU J 1 > UJ 3 1 > t- Tj •* -J i— ,ii _l 1- ^ •-r -J ►- ■c r» _j H- ..j J» -1 H- J -1 »- £ r _J •— ^,, p* J 3 r» < 3 f>« 1 D ^~ 4 3 *~ « 3 f- t :> h- ^ i ^. J> ■4 > ^. J- O 9 u* CJ ■J o> J L? CT 9 CT a O -j> OOO O O O O O a l_> -J U — _l O _l l_> ^- _i O _l u -1 _) _l -1 a o O OOO O O O O O O O i* ^2 j < l^ O 4 ■ n 'J 4 1/1 CI < r> O -t LO -J 4 -« 'J l»l r lA O r- T) P ~t r\l 1*1 r m O • IU > • n it! > • n u > X « LA Xi > lY « A UJ > n UJ > • A UJ > IX LA ui > — X o o O OOO '-J ^J -^ — < — < W-* ■4 — 1 '■j -J u- IN _i UJ aj ^ IU IM -1 IU 'U -J •j _J Lu IM -J UJ IM -1 UJ o o o OOO O O O O O O T. t- 3 1 H- jj v 1— u a h- O 3. t- a. t— O Ci ►- a c. >- • * ^ • jj. Z 1 jr ^ • T. £ • Z z • 1 z • I 2 • X z r^ K X X < X X X X ^1 1/1 i/> .1 vt ul •/I in O 'J ' 1 !■) O O O >-> .) O O O 37 Pg 01 ISI -1 HI nj * -« * — in o ■0 -0 o -c 3 J3 3 P"0 P- 3 o o 3 3 3 3 o o o 3 3 o o o o o 3 O 3 3 O 3 O O ♦ * * ♦ ♦ + ♦ ♦ * * * ♦ * 1 * ♦ + ♦ * ♦ + + ♦ ♦ UJ III UJ UJ UJ UJ LU LU UJ UJ UJ UJ UJ UJ LL UJ LU LU UJ Ul LU LU UJ UJ pg UN .0 o r*^ — < f*i o c t> O —4 3 rr, 3 ♦ o -r o 0- O -< 3 ci- oo -« JN — 1 — l«| — o m ■e 3 X O 9 O »* 3 in 3 o m O in o —■ H in -rt 9 33 -O ISI ^j 3 p- o -b o PsJ O in O r*> oo O X * -si n r- n 3 r- f\) — * o in in O P- 3 m 3 HI o 4- OP- O •0 1-\ J- »1 o ^* o o m *^» n- X o NJ O -*• 3 X 3 3 o o O - 3 *> n — i m 9 > i\i n >vj -^ 3 —j —J "l 3 3 o m -^ >n m 3 3 pg 3 3 9 3 ot» 3 n ■n o J- (•1 3 o — < O -c ^ o n aj 3 ^1 T 3 0< i^- ,A * 3 rr, ■• X ■o * Nl — -U Nl cr UJ O ^J UJ s0 in UJ in X UJ rr\ CT- UJ -a n ai ■0 J" UJ O Nl LU -n pg LU p- -!■ UJ CP n UJ o T x •NJ s -3 o •— ■O -< ^i B" —4 -!• CD m M s» "- J3 n ^p X X pg in HI p- p» Ml ~r O * o -^ o * y> ^f~ ■o -vj 3 * p- X X l»3 9 ~r * — t NJ -r o J- r^ UN O X P- J- o n m X u .N M ^ 3- j^ ~J •o -■ s n -0 -J" P* X vT t 3 P- o p- n o —j 3 3 f^ ■t 3 -» '3* JN o HI o l«3 *» "— » m J^ "V ? ^ X —i — * pg N 4-> Z3 Q. 4-> o M Ml Nl -n «l rsi •t -* -»■ W4\ in o O ■0 J3 ^^ J3 — 1 p--i p- -4 o O o o 3 3 o o O 3 3 3 O 3 o 3 O o 3 3 oo O o « * ♦ * * ♦ ♦ ■♦ ♦ ♦ + ♦ * 1 * 1 ♦ 1 ♦ 1 + 1 * 1 ■o UJ uj u. UJ aJ UJ UJ UJ UJ UJ Uj UJ UJ UJ UJ ul Ul UJ JJ LU LULU UJ UJ -* 3 CO 1*1 IN ~4 3 n 3 * 3 -J- 3 3 u 3 om o in * N -A 9 r* -r r~ ^- -^ X in rPi 3 ni 3 ■", 3 n 3 ■0 O-i o ^ CD rn UN ■A r^ * •0 O 9 ro m i^ ^J U M O LP 3 UN 3 UN O o J3 U Sm 9 c <3 n O m T j- ^1 m im ~* O in —i HI n 3 3 X M UN --1 ■4 3 3 ~j 3 -* M| TM M| m UN _i ■f ■ 3 C -^ o X 3 O m 9 o ■-M -a 3 CP o 3 n 3 o P- X O 3 L"l 3 ,N 3 O O X o UN a gr B- *• -»• O ♦ ~g + -f -o ♦ X 3 * jn NJ + p- PH 1 X ■M I 3 m 1 r^. m 1 O 3 1 X X I 3 NJ X UJ -v o UJ in UJ O 6 UJ in o UJ ■o ON UJ -n n UJ -0 —4 UJ O X LU rO J3 JJ P- i* UJ CP -J- Ul O s n Nl M o nj in -n p- ^3 * -* NJ •T X -H 1 1 O L» -si ~J CO J- p* o p- ^ pj LP ^>• NJ Nl 4- •^ -1 r* -3 ^* X X —j o 3 m X *~ 9 M o LI n .* — • — 9 .J-v •O ■T X X * X 3 m -4 f^ —1 P- — . ■J3 t> HI UN Mi J- X CP J3 JO 3 •Pi r- "- 9 J^ rt |\| 3 *M co O rr\ m4 -T -^ ^/ p- r- NJ -r i«| P^ J- ^ t 9» HI 3 UN J- -3- pg NJ NJ H| ■» z 1*1 X — — » --i •^ -J- * IM — < 9 0« -* ^3 n n •^ 3 ■& X o ■^ 14* -r —» i»l X 3 J3 —4 -0 O r*3 3 a 3 O O •u 3 CP P- LN o r— a m Qu o ~i J -1 J- 3> o >»■ o 3 O % O O o rg J> X On iTi in 3 o 3 -0 o -* X ^ O JN in -r O O LT J3 3 ■f vT * 3 J- O P- o m ~4 r- N n. Uj > _ _ LU i»l «1 0\ UJ ^ r*» r* UJ O -< -H UJ NJ -T "1 UJ ■J I", PJ LU 1*1 X X UJ ■O O P- JJ 3~ J- p- UJ -4-BI P- UJ Nl ■* p» t— O 1 1 1 o 1 3 I 1 I 3 1 X 1 1 | X 1 + 1 0» 1 -* 1 — 4 1 PSI 1 -T 1 1 (_>-. o in O 3 X X o 3 HI -D 3 JN o 2: 3 ~ n 3 —4 X -o m J» U3 NJ o m •^ O f~ T vj- —4 j- 0~ U3 O n —• 3 1 a LA ■0 3 p- o a> in n -- C/5 N 3 •^ r^ D rr\ —4 3 X o- pg n -4 M J3 H r^ • • -^ •■ HI .. r- • • 9 ■ i p4 *• h4 • • X .. m M 3 • • ■*> * *r ^ « J- i — Jf t*\ -c (SI JL X a J- J. 3 JL p- 3. J3 JC t= O 3 N 3 X 3 n O •>JO ISI O o o -4 o bn o -PI 3 > o -r 3 f«. a IP -t X ac rA a -fUL m ex. oa JL J LO X L0 ISJ LA P- in o ia O in —4 V^ gP in 0P J"> -* — 3 ■t u 2 •50 OJ o n uj O ■n JJ o 1} LU O P«J IJJ a -13 JJ 3 3 u c -4 U.' 3 o jj ~ 3 u 3 o o» > Z 3 > z O > z X 2S n > ^ 3 ^> z — > 3 r> > Z O > :r 3 -> ^ o >z 3 > 3 u • -• 3 • — • 3 • M 3 • • M 3 • M4 3 • -4 3 • — . 3 t Ml 3 • M| 3 • — 3 • Ml rsi » 3 -« — o ro »— -) ■^ i-O -* ^. a •* u- LU JS h- 3 —4 1— a -- I— o —4 t— 3 -J 1- J *-g h- 3 1 4 3. < -0 1 'Oj UJ Li UJ UJ LA UJ UJ l/» U' LU L0 UJ u_ *s UJ jj u*> UJ UJ 1/3 u- Lu u0 JJ UJ in u> UJ 1/3 M-— . u: a a U •3 ( o o or o D or O o or o g cr U o a o 3 Ji o c cr- ■_> l3 or o a Or 3 3 ■• Z *" •• ^; -i •• z - * 3 a> » u," 3 LI Z UJ 3 UJ ^ UJ 3 -u 2 UJ 3 UJ z Ul 3 u 2 u, 3 u, jj i_j 3 UL' z jj 3 .u Z _^. 3 u. Z u. 3 LU .' U' 3 ■3 Z o J 3 3 3 c 3 3 '3 3 3 33 a 3 3 3 3 3 3 o J 3 3 3 3 3 3 3 '_j 3 lTj 3 3 3 z Q a 3 3(3 n 3 C 3 O r— -1 — UC -J M4 -E _f — . I j ■* 33 -J — . i „j ■4 33 -J b-. £ -J -b 33 _j — • J3 _j Ml rrj u — cc _l — ro Q. < ul Ik < ^ l^. < */) u. -3 L0 u> 4 lO u. < ^-> jj_ < LA UL 1 i/) IL -3 2 IT -> 1 ^ */* > ^ ~» »•» > o ^ L/> > 3 ^ ■n > 3 ~~ m p» 2 "■ i/l > 3 • 11 ^> c 7" -1 > 3 7" v~ > 13 7" 1/3 > 3 .* i/l Z w — » ^ O •■ ji d "- ^ 3 B- 2 3 — Z. _) — 2 J *-• z 3 •-• 2T 3 M" -JT a •-• Z w •-• ^ \Zj « IXJ ~ 3 U. w -J |/| — 3 O V 3 5 o i/l o 3u) O 1 > u; 3 J LO 13 3 c 1/3 >- 3 3 3 _> U3 > o 3 J i/3 3» O 3 u> -r\ C 3 o 1/3 v -J 3 O >■ a 3 u) X I JM i— J 1 ^ -JJ UJ UJ 1 V O 1 1 1 U i_l 3 1 LU r- _> »v -J r— £_ J L* ^ ■» -J ►— a: •» ^ •— -T 1" J 1- a ■» U i— -2 • ■ -J ►— i t •J h- .3 pv -J h- J3 JJ -j M» JJ K -J i- X •» j 3 •- < 3 -^ < -J >- 1 3 ^» < 3 N < 3 p* < 3 N 4 3 r*~ 1 3 ►«• •T 3 r* 4 3 *- 4 3 p- 4 3. 4 • Z CL I • ^f i < • ^ a 4 • ^ 0. j • £ a ** • ^ 3. < • 7 3. < • «. a < • -:• 3. < • *T 3- 4 • .T CL <: • •T J. JL J 4 ^ i. 9 < ^ i 9 ■1 t S 9 1 i- 5- ■y •a JT ^_ 9 *l K JL J" c E r 3> 4 z Z. »s -a X i CP 4 X S- 3 1 •a •. I J- 4 C J » 3 J > •3 >5 9 3 O 3> ^ 9 3 3 9 CJ O 0> JO 9 a 3 3- 3 3 LP Q OOP l_> 3 CP -J — —J o — _J u> — a -J <_> «■ -1 O _i UJ — _l -) MB -J J Ml -J u a— _l O — -J u> — J UJ — -J i/> L_ 4 I/ - . _j < •/> 3 4 >J1 _) < /i j < w0 a < ■s 3 < |J- 3 < in J *X .j-i a 4 U3Q ♦ u- > » LU > • LJ > • UJ > 1! ^j > * UJ > » UJ > • UJ > * X > » '.U > « LU > n — 1 . A — X *n ■M jr J> -• I n IM or in _ 'X a-i — ->r LA — X /I — 3C UN Ml 3. n — IT in — or J -i LU cj ^ -u AJ J U ■M -J LU -si — i LU NJ -I UJ NJ -J UJ rg •J JJ NJ -J UJ pg -J yj PsI -J UJ O a i~ ^ <^ »v c "i ^ 3 1 )»• 3 0. _ 3 a ►- 3 0l H» O Z. L- JO C ^~ o a. ►-• J 3. I— C 0. L- • Z -T • i. « • r 2: • jr i X ^ a z. • I z • 3t z • I Z • 1 Z • X z • I z -*l X X X X « X >< K X X X X l/> ^1 >/) ^-i uo L-> JS LP m in 1/1 J-l 3 £3 3 O a c o 3 3 3 3 3 U J U o j o ul u) i_> u) -i J 38 Final example: Lanczos' Algorithm Lanczos' Algorithm is a method for reducing a symmetric matrix to tridiagonal form. Without re-orthogonal ization, the vectors of the orthogonal transformation which it generates deviate gradually from orthonormality. The SNORTRAN coding of the algorithm was written before the problem with computed zeroes in statistical relative error analysis was fully appreciated, otherwise it would never have been done. It is included here as an example of how the SNORTRAN system was intended to be used (the PL/ I preprocessor-output is included, since the algorithm uses so many SNORTRAN features), and also as an indicator of how severe a problem the irrepresentabil ity of zero is. The program is initially given v.. =(1,0,0,0,...) and is asked to generate 63 other mutually orthonormal vectors. Statistical bounds for are huge because the result is so close to zero, and during computation of the program blows up because a zero is actually encountered. 39 C0100 00200 00300 00*00 00500 00600 007OO 00800 00900 01000 01100 01200 01JOO 01*00 01500 016OO 017OO 0180C 01900 C2000 02100 02300 02400 02500 02600 02700 02800 02900 03000 03100 03200 03300 03400 03500 03600 03700 03800 039OO 04000 04100 04200 04300 0*400 04500 C4600 04700 04800 C4900 05000 05100 05200 05300 05400 05500 05o00 05700 05dOO 05 900 06000 LANCZOS: PRGC OPT ICNS( MA IN) REORCEP; OECLAKE ( 1 tJtK.I I, IK, IX, IL,NLIMIT,N,P,NP) FIXED BIN, T FLCAT BINARYI53); SMEAR 0(64), V(8,64), rf(64), ALPHA, BETA, GAMMA; ShEAR VSAVE(3,64>, IPVALlb), TEMPU64), TEMP2164); DECLARE STRING (8) ChAR(13) I N IT ( • • , • • , •% «' , , «, •»,«', , «); INITIALIZE ROUNDED; N = 8; P=8; NP = N*P; /* INI T V(l) »/ CALL SETC0NST4NT(V(1,1),SJRT( I . 000000000 COOOOOOOOCOOE 08 /NP J ) ; DC <=1 TC NP; VSAVE(1,K) , VU,K) = Vll.l); END; /* /* Computation of the square root of the smeared number GAMMA INIT Ull) = A v( 1) DO 11 = 1 TO NP; TEMP1UI) = V(l,II); CALL AMbLTtO, TEMPI); MAIN LOOP DO IX=2 TO NP; IK = MODI lX-2,8) + l; DO I 1=1 TO NP; TEMPKI I) CALL INNERIALPHA, TEMPI, U) ; DO K=l TO NP; EVAL W(K) = U(K) - END; CALL INNfJP LIKE SMEARED_NUMBER, 1 ALPHA LIKE SMEARED_NUMBER , I BETA LIKE SME ARED_NUMBER , 1 GAMMA LIKE SMEARED_NUMBER; DECLARE i VSAVE(8,64) LIKE SMEARED_NUMBER, I IPVALI8) LIKE SMEARED_NUMBER, I TEMPK64) LIKE SME ARED.NUMBER, I TEMP2I64) LIKE SME ARED_NUMBER; DECLARE STRING (8) CHARI13) IN I T( • • , • • , •« f '* ,»' , , , ,«« , , « ); DECLARE 1 2 2 2 DECLARE SMEARED_NUMBER, COMPUTED_VALUE FLOAT VARIANCE FLOAT INTDATA, 3 IA_LB FLOAT 3 IA_UB FLOAT BINARYI21 ), BINARY<21), BINARYI21 ), BINARY121); /♦UNITS BETA**(-2T)*/ /* /* INTERVAL ANALYSIS DATA */ */ I REG(IO) LIKE SMEARED.NUMPER, 8ETA_MT FLOAT BINARYI21) INIT ( 1.0E-24B1 , HUGE FLOAT BINARYI21) INIT < 1 .0E+88B ) ; /*l6**(-6)*/ SMEARED_ADD: PRQC (A,B,CI; /* PERFORMS ADDITION A = B + C FOR SMEARED NUMBERS */ DECLARE I (A,B,CJ LIKE SMEARED.NUMBER , IT,T1,T2,T3) FLOAT BINARY153); T = B.COMPUTED_VALUE; A.COMPUTED_VALUE = ROUNDIT ♦ C .COMPUTED_VALUE 1 ; T2 = MAXSO(B. INTDATA); T3 = MAXSQIC. INTDATA); CALL INTERVALU. INTDATA, 1 /*ADD*/, B. INTOATA, C. I NTD ATA) ; Tl * MINSOU. INTDATA); IF KT2+T3I/T1) > HUGE THEN DO; A. VARIANCE * HUGE; RETURN; ENO; A. VARIANCE = ROUNDUP! T2/T1*8.VAR I ANCE ♦ T3/T l*C .VARI ANCE CALL ADD_ROUNDING_ERROR(A); END SMEARED_ADD; ) ; 42 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 1 2 1 2 1 2 SMEARED_NEGATE: PROC (A, B) ; /* PERFORMS NEGATION A » -B FOR SMEARED NUMBERS */ DECLARE 1 (A t B) LIKE SMEAR ED.NUMBER, T FLOAT BINARY(21); A.COMPUTED_VALUE = -B.COMPUTED.VALUE; A. VARIANCE = B. VARIANCE; T = B.IA_LB; A.IA.LB = -B.IA_UB; A.IA_UB = -T; END SMEARED.NEGATE; SMEARED.SUB: PROC (A,B,C); /* PERFORMS SUBTRACTION A = B - C FOR SMEARED NUMBERS */ DECLARE 1 (A,B,C1 LIKE SMEARED.NUMBER , HUGE THEN DO; A. VARIANCE = HUGE; RETURN; END; A. VARIANCE = ROUNDUP! T2/T1*B. VAR I ANCE ♦ T3/T 1*C . VARI ANCE ); CALL ADD_RnuNDING_ERROR(A); END SMEARED.SUB; SMEARED_MULT: PROC(A,B,C); /* PERFORMS MULTIPLICATION A = B * C FOR SMEARED NUMBERS */ DECLARE 1 (A,B,C> LIKE SMEARED_NUMBER t T FLOAT BINARYI53) ; T = B.COMPUTED_VALUE; A.COMPUTED_VALUE = ROUNDIT * C .COMPUTED_VALUE I ; CALL INTERVAL(A.INTDATA t 3 /*MULT*/ , B. I NTDATA, C. I NTD ATA ) ; T = B. VARIANCE; A. VARIANCE = ROUNDUPIT + C. VARIANCE); CALL ADD_ROUNDING_ERR0R(AI; ENO SMEARED.MULT; SMEARED.DIV: PROCIA.B.C); /* PERFORMS DIVISION A * B / C FOR SMEARED NUMBERS */ DECLARF I (A,B,C) LIKE SMEARED.NUMBER , T FLOAT BINARYI53I; T « B.COMPUTED_VALUE; A.COMPUTED_VALUE = ROUNDIT / C . COMPUTED_VALUE ) ; CALL INTERVAL(A.INTDATA,4 /*DIV*/, B . INTDATA ,C . INTDATA) ; T = B. VARIANCE; A. VARIANCE = ROUNDUPIT ♦ C. VARIANCE); CALL ADD_ROUN0ING_ERROR(A) ; END SMEARED_DIV. SETCONSTAN T : PROC(A,C); /* SETS SMEARED NUMBER A TO THE DECLARE I A LIKE SMEAR ED_NUMBER» C FLOAT BINARYC53); A.COMPUTED_VALUE = ROUND(C); A.IA_LB = ROUNDDOWN(C) ; A.IA_UB = POUNDUP(C); A. VARIANCE = OEOB; IF C -= PRECISIONS, 21) THEN CALL ADD_ROUNDING_ERROP (A) ; END SETCONSTANT; INITIAL VALUE C */ 43 i o 2 86 2 87 2 88 2 2 2 1 2 INTERVAL: PROC ( A, OP, B, C ) ; /* PERFORMS INTERVAL ARITHMETIC A = B OP C */ DECLARE 1 INTERVAL, 2 IA.LB 2 IA_UB 1 (A,B,C) LIKE INTERVAL, (SL,S2,S3,S4) FLOAT BINARY(53), I0P(4) LABEL, OP FIXED BINARYU5); SI = B.IA_LB; S2 = B.IA_UB; S3 = C.IA_LB; S4 = C.IA.UB; GO TO IOP(OP); MAXSQ FLOAT BINARY(21), FLOAT BINARY(21), IOP(l): A.IA_LB A.IA_UB RETURN; IOPC2): A.IA_LB A.IA_UB RETURN; IOP(3): A.IA.LB A.IA_UB = MAX( /* ADDITION */ ROUNDOOWNfSl + S3); ROUNDUP (S2 ♦ S4); /* SUBTRACTION */ ROUNDDOWNIS1 - S4); ROUNDUP (S2 - S3); /* MULTIPLICATION */ = MIN( ROUNDDOWN(Sl ROUNDDOWN( SI ROUNDDOWN(S2 ROUNDDOWN(S2 ROUNDUP (SI ROUNDUP (SI ROUNDUP (S2 ROUNDUP (S2 * S3), * S4) , * S3) , * S4) * S3), * S4), * S3), * S4) RETURN; I0PI4I: A.IA.LB A.IA_UB * MAX( / S3), / S4), / S3), / S4) / S3), / S4), / S3), / S4) ); /* DIVISION */ MIN( ROUNDDOWNtSl ROUNDDOWN(Sl ROUNDDOWN(S2 ROUNDDOWN(S2 ROUNDUP (SI ROUNDUP (SI ROUNDUP (S2 ROUNDUP , (Tl,T2) FLOAT BINARYI53); */ IF /* INTERVAL CONTAINS ZERO */ RETURN(Tl); END; = HUGE I A.COMPUTED_VALUE=0E0B THEN RETURN; D = MIN(MANTISSA(A.I A_LB) ,MANTISSA(A.IA_UB) ); S = IE0B / 1.0000000000000000000000E08); Y * Y * 1.0000000000000000000000E-4B; end; do while (y < 1.0000000000000000000000e-a8); y = y * l.0oooo0o0o0oooo00ooo0ooe+*8; END; RETURN(Y) ; END MANTISSA; SMEARED_PRINT: PROC(AA,A); DECLARE I A LIKE SMEARED_NUMBER , AA CHAR!*), (S,S1*S2, ABS1,ABS2,ABS3,ABS*) FLOAT BINI53); S = SORTS. VARIANCE); 51 = -3.625 * S; /* 99.71 CONFIDENCE LIMITS */ 52 = +3.625 * S; A8S3 = A.COMPUTED_VALUE * (1 ♦ S1*BETA_MT); ABS* = A.COMPUTED_VALUE * (1 ♦ S2*BETA_MT); IF A.COMPUTED_VALUE < THEN DO; ABS1=ABS3; ABS3=ABS*; ABS*=ABS1: END; PUT FILE(SYSPRINT) SKIP EDIT J AA, 'COMPUTED VALUE:* , A.COMPUTED_VALUE } ISKIP,COL(7),A,X(8),A,E(25,15,16H ( »3.625«SIGMA BOUNDS ON RELATIVE ERROR: *, '(•, SI , • ,',S2, • I * BETA**C-T)' ) ( SKIP, COL( 121, A, COL (50) ,A,2 ( E ( 25 , I 5, 16 I , A ) I 45 41 2 42 1 43 2 *4 2 <»5 2 46 2 1 47 2 1 48 2 1 49 2 I 50 2 51 2 52 2 53 1 54 2 55 2 56 2 60 2 61 1 62 2 63 2 64 2 68 2 69 71 72 73 74 1 75 1 76 79 80 o. 81 1 82 1 85 1 86 1 IMMPLIED 99.7*-C0NFIDENCE BOUNDS: 1 , •(', ABS3, • ,*, ABS4, ' )• ) (SKIP, COL (12), A, COL (50), A, 2 (E( 25, 15, 16), A)) (MNTEPVAL ANALYSIS BOUNDS: 1 , • ( • , A .1 A_LB , • f «, A.IA_UB,« l«| (SKIP, COL (12), A, COL (50), A, 2 ( E ( 16,6 ,7 ) , X.< 9) ,A) ) ; END SMEARED_PRINT; ROUND: PROC(X) RETURNS! FLOAT BINARY(21II; DECLARE (X,HALF) FLOAT BINARY(53), FLOAT BINARYI21), BIT (64), BIT(8) DEFINED E POSITION(l), BIT(l) DEFINED E POS I TIONl 33 ) , BIT(64) INIT( • 0O0O00000OO0OOOO00OO000O0OOO00O01 OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO' B) , FEXPT BIT(8> DEFINED F POSITION(l); E = UNSPEC(X); RX E EEXPT RBIT F IF /* ROUND */ RBIT THEN DO; FEXPT = EEXPT; UNSPEC(HALF) = F; RX * X ♦ HALF; END; ELSE RX = X; /* TRUNCATE */ RETURN(RX); END ROUND; ROUNDUP: PROC(X) RETURNS (FLOAT BINARY(2D); DECLARE X FLOAT BINARY(53), RX FLOAT BINARY(21) ; IF X > OEOB THEN RETURN( ROUND ( X )) ; ELSE DO; RX=X; RETURN(RX); END; END ROUNDUP; ROUNDDOWN: PROC(X) RETURNS (FLOAT BINAPY(21)); DECLARE X FLOAT BINARY(53), RX FLOAT BINARY(21); IF X < OEOB THEN RETURN( ROUND( X) ) ; ELSE DO; RX=X; RETURN(RX); END; END ROUNDDOWN; N = 8; P=8; NP = N*P; INIT V( 1) */ CALL SETCONSTANT(V(l,i) ,SQRT( 1. 0000000000000000000000E03 /NP ) ) ; DO K=l TO NP; VSAVE(1,K) , V(1,K) = V(l,l); END; INIT U( 1) = A V( 1) DO 11 = 1 to NP; TEMPKII) = V(1,II); END; CALL AMULT(U, TEMPI) ; MAIN LOOP DO IX=2 TO NP*. IK * MOD(IX-2,8) ♦ l; DO 11 = 1 TO NP; TEMPKII) * V(IK,II); END; CALL INNER(ALPHA, TEMPI, U); DO K=l TO NP; /* /* /* */ */ 46 187 188 189 190 191 192 193 195 197 198 199 200 201 2 02 205 206 207 208 209 210 212 213 214 215 218 221 222 223 224 225 226 227 228 229 230 231 2 2 2 2 2 1 z 232 233 2 34 2 2 2 z 2 2 235 2 2 2 36 2 2 237 2 2 238 239 ?4,0 2 2 2 2 1 CALL CALL CALL SMEARED_MULT{ REG< 1 ) , ALPHA, V( IK,K ) I ; SMEARED_SUB1 THEN CALL SMEARED_SUBIU(K),U(K) ,V(K-1) ); IF KN THEN CALL SMEARED_SUB(UCK) ,U(K),V(K«-1)); IF J>1 THEN CALL SMEARED_SUBCU(K) ,UCK) ,V(K-N) ) ; IF J

NI ); END; END; END AMULT; 47 UJ CD » * < UJ CD UJ CD I » < CO — o O O * * UJ UJ o m o» m U1 ec p- o> k- m o o» in >*■ tr\ O O «. O O in o o moo f*> o o in -j ♦ 4 UJ UJ i»» o — o Ml «<■ •a o o> r- >o ^o f- >o m o — i ,0 1 CO — < UJ (M 0* O mi >0 >t f>- t\j UJ co a o> mi p. in o^o — -f i C* n uj m -^ (\i a m o *. (M o- f\j inj >C i*> p- ec O O r-i • t • •a >t >* -t-> s- o •* a. — . Z z * Pm %■» a o > •• » •» t- » m* a or Q •« — * UJ UJ UJ *-» *: -i • P or CO CO X II w w p» UJ X a: *- * 3 O MP > CO DDU. ► UJ Q p 3: z z o — oc Z ■# a 3 1 1* a — 1 ► UJ OC • z a d — a. UJ a. l UJ UJ UJ •• O Ml z _^ a or a. a: O UJ » z M UJ < < < a. a »* u af UJ UJ _J II •M M4 ■ • o < s ru K — § \r> a UJ VO IA UJ Q. -1 o o a i/l UJ UJ r> a •x < I I UJ KJ u z •• UJ *— *■■* o a < ec M -J -J UJ UJ mi UJ oc cc z _i *■» <-» < < o z * * UJ UJ z m» UJ < u UJ o a IM Ml D > a: t mi -J - -1 < < UJ o o o o o m* —1 t-« o o art (VJ «NJ «M <\j m o * * * •# «* <# * * ^ Ifl N N NNN «sj rg ivj pj pg »-4 O cr in r- mi er in D_ o o (NJ —4 o o (NJ -^ + •f ♦ •f «v | + ♦ OO O IVI O UJ UJ UJ UJ UJ UJ UJ UJ 0* «t r- o O Ml f\ .*• o 0^ t— < o O (NJ rs >*■ D» IM Ml -J- O 0- m r- in O o o O P- o- cr- 2: o r~ m r- 0> r» m4 o 0> -0 o- o> C u^ o 1 o h» o I o o «i- m (M 3- i0 o CO r~ in (\J in. -* «3 00 m m o r>- w ^* s f- ^> in <«■ r» n >}■ IN h« (NJ Ml >- QJ O O P- in o o> •o O m r« Ml f^ P- o< «\ (NJ (NJ o a O ro o o> m in N o CO o t 2 in t • • CO • • • Ml • • • m • • • UJ o> ^i ~4 o *■ •». o_ o fl c* «M c -Jt o «« mi «* •* Mi Ml *m o> w *M M* rr\ mr w am UJ o o o m* vo UJ . . # •• • • • t • I • • • Ml a: or Ml Ct < 3 o 1 D D D ec Q Q- E X a a or or UJ a: a or ct z u UJ • « UJ • • UJ • « u • • □ O •■ c« at \r, •t oo »• oo M* er UJ UJ D UJ UJ Q UJ UJ Q UJ UI O ►- Q- UJ r> > Z Z> > 2 rs > Z n > r» — _j M 3 ml .— D _j *— 3 _i IM Z> a 2 ^— < ►- C < H C < t- o < V- D z Ml ro > < CO > < CO > < CO o C _J •« -I 1* -J •• _J ■t o a •i — o UJ u l/l o u UJ V D U' Ui 00 o 111 LO OO cc u_ UJ K CJ n U-l cr o c UJ or u o UJ cr L' O • (M t- ?■ 2 K 2T z h- Z Z t- *" 2 Jl O ■"> 2 U' ~> T> Z UJ ~ * Ti z UJ => 3 7 UJ -J o o a a O O 0. D O o 0. n o a c C Q O _l c 5" mt CO y ^-< CO y: m. cr s: M CO u 4> O i^ u. o i/i U- c io u. C) o" a or o r: V. ul o o ^ 00 o n Z oo o c: Z oo UJ t- i; o ^« i: O M« iT o — 2 c Ml > UJ n o l/» 3 i^ 00 3 l_) oo 13 O 00 C O0 o 1 V a 1 V C 1 > r 1 >- • u. cr M —I < co P- -J < cc K -J co _l < u. a L" O- r^ < ^^ -J A _- _J Ml mI A *— -J O i/i Q •J ma on o < OO o < -~ 1/1 e < II o A * UJ > — 1 ■» m > A • UI > ^^ « UJ > ■ (NJ «» IT *— or 1 in »—. o «■ in m. a 1 IT -~ cf U' -J tM _J UJ -5 (M _j UJ "J rsj _l UJ -> rv -J Uj o ►- w* ^ a ►- m» 49 c h- «M •o o y- *^ >0 a >- c J* a > • s; z > I r Z > i T. Z > t 3" Z (J UI a •■ tn P r: z 3 <3 *~ mi •P* ■»» o U 1 -> -> 1 -> » 1- c ^ ** •JM* ^ < UJ > > > > ►- cr f\l V V f« v V ~0 oo < Uj a 00 M It CO "> -) M. 48 SNORTRAN System Listings 1. SNORTRAN Preprocessor (SNOBOL) pp. 49-55 2. SNORTRAN System Routines (PL/I ) pp. 56-60 (Referenced as file ' SNORT. RND' in line 09000 of preprocessor) 49 00100 C0200 C0300 00400 00500 00600 00700 00800 00900 01000 01100 01200 01300 OIhOO 01500 01600 01700 01800 01900 02000 02100 02200 02300 02400 02500 02600 02700 02800 02900 03000 03100 03200 03300 03400 035C0 03600 03700 03800 03900 04000 04100 C4200 04300 04400 04500 04600 04700 04800 04900 05000 05100 05200 05300 05400 05500 C5600 05700 05800 05900 06000 »-- * * * * * — -EJECT * * A POSTERIORI ROUNDOFF ERROR LANGUAGE PREPROCESSOR 6TRACE = 1000 DEFINE('FIL..( VECTOR, STRING II ,VAL« ) DEFINED PARSE. EX PR (II t J t LOOK AHEAD. TOKEN, LOOK AHEAD. SYMBGL • I OEFINEl 'SEMANTICS! PR C0NUM1 •» DEFINE! 'GETCHAR! )« ) DEF1NE( 'GETNCNBLANKd •} DEFINE( 'SCAN! !• I DEFINE( 'ADD( !• I DEFINE! 'REMOVE! ) ■ ) OEFINE! 'EKIT(STRING)' I DEFINE! • LASTEM IT (TARGET) XtY* J DEFINE! 'ALLOCATE D.REG I II • I DEFINE! 'FREE. REGI ST ERl II' I DEFINE! 'CLEAKREGS( II' 1 KEYWORD = ('INITIALIZE' I 'SMEAR' | »EVAL' I 'PRINT'! . KEY NREGS = 10 REGISTER = ARRAYtNREGSI SEMICOLON = 1 VARIABLE" = 2 EQUALS = 3 PLUS = 4 MINUS = 5 STAR = 6 SLASH = 7 LPAREN = 8 RPAREN = 9 FACTOR = 10 TERM = 11 EXPR = 12 ASSMT = 13 ACCEPT = 14 COMMA = 20 SLR(l) EXPRESSION PARSER TABLES DEPTH = 50 STATESTACK = ARRAY(DEPTHI SYM30LSTACK = ARRAY(DEPTH) TOKENSTACK = ARRAY(DEPTH) STORED. DEPTH = 5 STORED. TOKE;M = A RR AY t STOR EC. Cb'PTH I STORED. SYMBOL = ARRAY { STOP ED. DEPTH I NPRJDS = 12; NTRANS = 76; NSTATES = 23 3, 1, 2, 4, 5, 8, 10, 11, 12, 2, d, 10, 11,2, 8, 10 ,11, 2,' TS = '2,13 1,4,5,1,4,5,6,7 •4,5,8, 10, 11,12,1,4,5,0,7,9, uttJiitti '9, 1, 4, 5, 6, 7, 9, 4, 5, 9, 2, 8, 10, 2, 8, IU, 2, 8, 10, 11,' •2, 8, 10 ,11, 1 , 4, 5, 0,7, 9, 1 , 4, 5, 6, 7, 9,' TA = ' 2,3,4,0,5,6,7,0,9, 10,11,5,8,9, 12,5,8,9,13, 5, • •6, 7, 3, 9, 10, 14, -5, -5, -5, 15, lu, -5, -2, 17, 10,-6,-6,-6,15, 16, • •-o, -7, -7, -7, 15, io,- 7, 17, 18, 19, 5, 8, 20, 5, 8, 21, 5, 8, 9, 22,' •5,8,9,23,-3,-3,-3,15,16,-3,-4,-4,-4,15,16,-4,' FT = '1,3,4,5,0,12,16,20,0,27,33,36,42,48,51,54,57,61,0,0,' 50 06100 ♦ •0,65/71,' 06200 LT = '2,3,4,11,-12,15,19,26,-10,32,35,41,47,50,53,56,60,64,-11' 06300 + ' ,-3,-9, 70,76, • 06400 LS = '14, 13, 12, 12, 12, 12, 12, 11,11, 11, 10, 10,' G6500 RL = '2,3,3,3, 1,2, 2, 3,3,1,3,1,' 06600 TRANS. SYMBCL = ARRAY ( NTRANS) ; F ILL( .TRANS. SYMBGL , TS 1 06700 TRANS. ACTION = ARRAY 1 NTRANS ) ; F I LL( .TRANS. ACT ION ,T A ) 06800 FIRST. TRANS = ARR AY( NSTATES ) ; F I LL( .F IRST. TRANS , FT) 06900 LAST. TRANS = ARRAY( NSTATES) ; F ILL( .LAST .TRANS , LT ) 07000 LHS. SYMBOL = ARRAY ( NPRODS J ; F ILL ( .LHS . SYMBOL , LS ) 07100 RHS. LENGTH = ARRAY( NPRODS ) ; F 1LL< .RHS .LENGTH, RL ) 07200 ♦ :(REAO) 07300 * 074C0 FILL 1=1+1 07500 STRING SPAN! • 0123456789- • ) . VAL ',' = :F(RETURN) 07600 ITEMUVECTOR,I) = CONVERT ( VAL ,. I NTEGER) :(FILL) 07700 * 07800 * 07900 READ CARO = INPUT :F(END) 08000 RESUME CARO dP KEYWORD :S(PREPARE) 08100 OUTPUT = DIFFER(TRIMICARD) ) CARD :(READ) C8200 * 08300 PREPARE CARD LEMPI . SNARK = 08400 OUTPUT = IFFERl TR IM( SNARK ) ) • • SNARK :($KEY) 08500 * 08600 * INITIALIZE — CUPY IN SYSTEM DECLARES AND SUBROUTINES 08700 * DETERMINE ALSO ROUNDING MODE. 08800 * 08900 INITIALIZE SCAM) 09000 INITFILE = 'SNORT. RND • 09100 INITFILE = I DENT ( SCAN (),' CHOPPED • ) 'SNORT. CHP' 09200 INPUTI.INIT, INITFILE) 09300 INITCOPY OUTPUT = INIT :S(INITCOPY) 09400 SCAN!) 09500 OUTPUT = NE ( SCAN TYPE , SEM ICCLON ) • /***** ERROR IN INITIALIZE' 09600 ♦ • STATEMENT *****/• : (RESUME) 09700 * 09800 * SMEAR — GENERATE DECLARATIONS FCR SMEARED VARIABLES 09900 * 10000 SMEAR SCAN() 10100 OUTPUT = • DECLARE' 10200 SMEARLOOP V = SCAN( ) 10300 EQ(SCANTYPE, VARIABLE) :F(SMERRCR) 10400 SCANl ) 10500 EQ(SCANTYPE, COMMA) :F(SMEARFIN) 10600 OUTPUT = ' L ' V ' LIKE SMEAREU_NUMB ER , • MSMEARLOOP) 10700 SfEARFIN EQ ( SC ANTYPE , SEM I COLGN ) :F(SMERROR) 10800 SMcNC OUTPUT = • 1 • V ' LIKE SMEAREC_NUM3ER; ' KRESUME) 109J0 SKbkRCR OUTPUT = • /***** ERKCR IN SMEAR LIST. END FORCED *****/• 1100U : (RESUME) 11100 • 11200 * PRINT — FCRMATTE> PRINT CF SMEARED VARIABLES 11300 * 11400 PRINT SCANO 11500 P> LIST V = SCAM ) 11600 E<*( SCANTYPt .VARIABLE) :F(PRERROR) 11700 OUTPUT = » CALL S ME A RED_PRI NT ( • " V "• , " V ");" llbOO SCAN( ) ■ v-0 EfcHSCANTYPE, COMMA! :S(PRLIST) 120 E0( SC^NTYPE, SEMICOLON) :S(RESUME) 51 12100 PRERROR OUTPUT = • /***** ERROR IN PRINT LIST. END FORCED *****/• 12200 MRESUME) 12300 -EJECT 12400 * 12500 * EVAL — PARSE EXPRESSICNS OF SMEARED NUMBERS AND 12600 * GENERATE APPRCPRIATE SUBROUTINE CALLS. 12700 * 12800 EVAL SCANO 12900 PARSE. EXPRO :(RESUME) 13000 * 13100 SEMANTICS :($(»PRUU* PRODNUM)) 13200 * * 13300 * l: S ==> ; * 13400 * * 13500 PRODI :(RETURN) 1370C * 2: ==> = * 13800 * * 13900 PRCD2 LASTEMIT(TGKENSTACK) 14000 CLEARREGS1) :(RETURN) 14100 * * 14200 * 3: ==> ♦ * 14300 *■ * 14400 PKLD3 FREE . REGl S TEk ( TOKENSTACK J 14500 FREE.REGISTER(TCKENSTACK) 14600 AREG = ALLOCAT ED .REG( ) 14700 EMIT! ■ CALL SMEARE D_ADD( ' AREG • ,' TOKENSTACK 14800 + •»■ TCKENSTACkXSTACKPTR + 2> •);• ) 14900 TOKENSTACK = AREG : (RETURN) 15000 * * 15100 * 4: ==> - * 15200 * * 15300 PRCD4 FREE.REGISTERJ TOKENST ACK ) 15400 FREE.REGISTERiTOKENST ACK ) 15500 AREG = ALLCCAT ED .REG ( ) 15600 EMIT! ■ CALL SMEARED.SUB ( • AREG •,' TOKENST ACK 15700 ♦ *t' TOKENSTACK •);' ) 15800 TOKENSTACK = AREG : (RETURN) 15900 * * 16000 * 5: ==> * 16100 * * 16200 PR0C5 MRETURN) 16400 * t>: = = > + * 16500 * * 16600 PFGCo TUKENSTACK = TOKENSTACK : (RETURN) 16800 * 7: ==> - * 17000 PRCD7 FkEE. REGISTER! TfiKENST ACK) 17100 AREG = ALLGCATED.R^Gt ) 17200 EMIT( • CALL SME ARED_NEGATE ( • AREG •,' 17300 + TCKENSTAC«STACKPTR ♦ 1> •);• ) 17400 TuKENSTACK = AKEG : (RETURN) 17500 * * 176C0 * b: ==> * * 17700 * * 17800 Pt-008 FKEE.KEGISTFR(TGKENSTACK) 17900 FFEE.RcGISTER( TUKENSTACK ) 18000 AREG = ALLCCATED.REGt I 52 18100 EMIT! » CALL SMEA RED.MULT ( • AREG •»• TOKENSTACK 18200 + •»• TCKENSTACK ■)!' I 18300 TOKENSTACK = AREG '.(RETURN) 18400 * * 18500 * 9: ==> / * 18700 PR0D9 FkEE.REGISTER(TUKENSTACK) 18800 FREE.P.EGISTER{TOKENSTACK) 18900 AREG = ALLOCATED . REG( ) 19000 EMITi • CALL SMEARED_DI V ( • AREG ',' TOKENST ACK 19100 + ' t' TOKENSTACK •);• ) 19200 TOKENSTACK = AREG : (RETURN) 19300 * * 19400 * 10: ==> * ' 19500 * * 19600 PRCD10 : (RETURN) 19800 * 11: ==> ( ) * 19900 * * 20000 PRGOU TCKENSTACK = TQK£NSTACK : (RETURN) 20100 * * 20200 * ==> * 20300 * * 20400 PK0D12 EQIVARINFO.O) :S(RETURN) 20500 CCNSTANT.VAK AREG = ALLCC ATEO.REG ( ) 20600 EMIT( « CALL SETCONST ANT ( • AREG •,' 20700 ♦ TOKENSTACK •);' ) 20800 TOKENSTACK = AREG KRETURN) 20900 * * 21000 ALLOCATED. REG I = LT(I.NREGS) I + 1 :F(ALL.USED) 21100 REG1STER = I DENT( REGI STER< I> ) 'X' : F( ALLOCATED .REG ) 21200 ALLOCATE ALLOCAT EC. REG = 'REGC I •)• KRETURN) 21300 ALL. USED OUTPUT = • /***** ALL REGISTERS USED *****/' :{ALLOCATE) 21400 * 21500 FkEE. REGISTER REGNAME «REGC BREAKC)') . I :F(RETURN) 21600 REGISTER = : = : (CLEARREGS) 22000 * 22100 EMIT OUTPUT = D I FFER ( SA VED .STR ING ) SAVED. STRING 22Z00 SAVED. STRING = STRING KRETURN) 22300 * 22400 LASTEMIT SAVED. STRING (BREAM'(') LEN(l)) . X BREAKC,') = 22500 ♦ X TARGET :F(ASSMT) 22600 OUTPUT = SAVED. STRING 22700 SAVEC. STRING = : (RETURN) 22800 ASSMT OUTPUT = • • TARGET « = • TOKENST ACK •;• 22900 ♦ :(RETURN) 23JC0 -tJECT 23100 * * 23200 * 23300 * SLR(l) EXPRESSION PARSER 2340U * 2 3 500 * * 21600 PARSE. LXPk ST AT EST ACK<1> = 1 23700 SYMBCl STACK<1> = 23d00 CURRENT. STATE = 1 23SJO STACKPTk = 2 2*000 NfEDPEAD ■ «YtS' 53 24100 24200 24300 24400 24500 24600 24700 24800 24900 25000 25100 25^00 25300 25400 25500 25600 25700 25800 25900 26000 26100 26200 26300 26400 26500 26600 267C0 26800 26900 27000 27100 27200 2730C 27400 27500 27600 27700 2780C 27900 28000 28100 28200 28300 28400 28500 28600 28700 28800 28900 29000 29100 29200 29300 29400 29500 29600 29700 29800 29900 30000 * PARSE FIND FCUND * SHIFT I = FIRST. TRANS J = LAST.TRANS EU(I.U) NLEDREAO = 1 DENT ( NEEDRt AD ) 'YES* LUCKAhEAU. TOKEN = SCANO LOOKAHEAC. SYMBOL = SCANTYPE EQILCUK AHEAD. SYMBOL, T RANS . SYMBOL <1 >) I = LT(1,J) I ♦ 1 J = TRANS. ACTIGN EQ(JtO) LT[J,0) S(REDUCE.LRO) S(FIND) S(FOUND) S(FINO)F(SYNTAXERR) S(ACCEPT) S(REDUCE.LRl) REDUCE. REDUCE. Rfc'DUCc GCTO PUSHGOT STATESTACK = J SYMBOLSTACK = LOOKAHEAD. SYMBOL TOKLNSTACK = LOOKAHEAD. TCKEN CURRENT. STATE = J STACKPTR = LT(STACKPTR, DEPTH) STACKPTR ♦ 1 :S(PARSE) OUTPUT = • /******* STACK OVERFLOW ****»*/» : (PARSE) LRO NEEUREAD = 'YES 1 UREDUCE) LR1 NEEDREAD = PROD = -J STACKPTR = STACKPTR - RHS . l.ENGTH CURRENT. STATE = STATESTACK SYMBOLSTACK - LHS. SYMBOL I = FIRST. TRANS J = LAST.TRA,>iS EQ(SYMdLLSTACK,TRANS.SYMBUL) :S(PUSHGOTO) I = LTU,J> I + 1 :S(GOTO)F(SYNTAXERR) J = TRANS. ACTION STATESTACK = J CURRENT. STATE = J SEMANTICSIPROD) STACKPTR = STACKPTR + 1 MPARSE) * ACCEPT SYNTAXERR OUTPUT = * * -EJECT * SCAN — * * SCAN /***** SYNTAX ERROR : (RETURN) *****/t :(FRETURN) LEXICAL ANALYZER. RETURNED VALUE IS TCKEN STRING. ALSO RETURNS SCANTYPE — NUMERICAL TOKEN CODE VARINFO — EXTRA INFC FOR TOKENS SCAN NEXTC GETNL 1 HAR = NBLANKl ) :($( 'CLASS* CLASSNJ) * * * CLASS 1 — LETTERS, BREAK CHARACTERS: BUILD VARIABLES ClASSl SCANTYPE = VARIABLE; VAPINFC = CLASS1A ADK); GETCHAR() GE(CLASSN,l) LE(CLASSN,2) :S(CLASS1A) CLASSlb IDENT(NEXTCHAR, ' •) :F(CLASS1C) GET.\Cr^BLANK( ) CLASSIC Ew(CLASSN,b) :F(OUTI PARENCCUNT = 1 CLASSIC AuOO J GETCHAR( ) 54 30100 PARENCOGNT = E0(CLASSN»8) PAPENCOUNT • ► 1 :S(CLASS1D) 30200 PARENCOUNT = EQ(CLASSN,9) PARENCOUNT - - 1 :F(CLA3S1D) 30300 EUPARENCOUNT.O) :F(CLASS1D) 30400 30500 30600 30700 30800 AOD( ) :(OUT) * CLASS 2 — DIGITS: 8UILD CONSTANTS CLASS2 SCANTYPE = VARIABLE; VAR1NFG = 1 30900 CLASS2A APC(J; GETNONBLANM ) 31000 VARINFO = VARINFO + 1 31100 EQ(CLASSN,21 :S(CLASS2A) 31200 VARINFO = ICENT(NEXTCHAR,' .• ) 100 ;S(CLASS2A) 31300 IDENTINEXTCHAR, 'E' ) :F(GIT) 31400 AODO; GETCHARU 31500 G£(CLASSN,2) L5(CLASSN,4) :F(OUT) 31600 CLASS2B ADOO; GETChAR{) 31700 31800 31900 32000 EG(CLASSN,2) :S(CLASS2B)F(0UT) * CLASShS 3-10 — TERMINALS 3+ 4- 5/ 6* 7= 8( 9) 10; 32100 CLASS3 ACO(l; GETCHARl); SCANTYPE = PLUS 32200 OIFFEKIUNARY) EQ(CLASSN»2) :S(CLASS2)F(0UT) 32300 CLASS4 ALOU; GETCHAR(); SCANTYPE = MINUS 32400 3IFFERIUNARY) £0 (CLASSN, 2 ) :SICLASS2)F(CUT) 32500 CLASS5 ADD(); GETCHAR(); SCANTYPE = SLASH 32600 EU(CLASSN,6) :F(OUT) 32700 CLASS5A REMQVEO; GETNONBLANKl) 32800 CLASS5B EU(CLASSN,6) :F(CLASS5A) 32900 REMOVED; GETNONBLANKO 33000 EQD( ) . • SCANTYPE = EQUALS :(0UT) 33400 CLASS8 AUDI ) ! , SCANTYPE = LPAREN :(OUT) 33500 CLASS9 ACDU , SCANTYPE = RPAKEN :(OUT) 33600 CLASS10 ALU! ) • SCANTYPE = SEMICOLON : (OUT) 33700 CLASS11 ACOU . SCANTYPE = COMMA :, I (A,B,C) LIKE INTERVAL, (S1,S2,S3,S4) FLOAT BINAkY(53), I0P(4) LABEL, OP FIXED BINARY(15); Si = B.IA_LB; S2 = B.IA_UB; S3 = c.ia_lb; GO TO ICP(CP) ; */ */ S4 = c.ia.ub; IOP( 1): A.I A_Lb A.IA_UB RETURN; I0P(2»: A. IA_LB A.IA_UB return; IOP43) : A.IA_LB /* ADDITION */ ROUNDDCWNlSl ♦ S3) ; ROUNDUP (S2 *■ S4) ; /* SUBTRACTION */ ROUNDDQWNISl - S4) ; ROUNDUP (S2 - S3) ; /* MULTIPLlCATlbN */ MINI RGUNDLUWNtSi * S3), RCUNUU0WNCS1 * S4 ) , R0UNCDGWN1S2 * S3), R0UNLCJkN(S2 * S4) 58 12100 12200 12300 12400 12500 12600 12700 12800 12900 13000 13100 13200 13300 13400 13500 13600 13700 13300 13900 14000 14100 14200 14300 14400 14500 14600 14700 14800 14900 15000 15100 15200 15300 15400 15500 15600 15700 15800 15900 16000 16100 16200 16300 16400 16500 16600 167C0 16800 16900 17000 17100 17200 17300 17400 17500 17600 17700 17b00 179 00 16J00 PAXSQ: MINSU: A.IA_U8 = MAX( ROUNDUP (Si * S3), ROUNDUP (SI * S4) , RUUNDUP (S2 * S3), ROUNDUP (S2 * S4) RETURN; I0P(4): /* DIVISION */ A.IA_LB = MIM ROUNDDQwN(Sl / S3), R0UNDD0WN(S1 / S4), RCUNDD0WN(S2 / S3), R0UNDD0WiM(S2 / S4) A.IA.UB = WAX( ROUNDUP (SI / S3), ROUNDUP (SI / S4), ROUNDUP (S2 / S3), ROUNDUP (S2 / S4) RETURN; END INTERVAL; PROC(A) RETURNSIFLOAT BINARY(53)); /* RETURNS MAX VALUE OF A**2 OVER INTERVAL A */ DECLARE 1 A, /* INTERVAL */ 2 LB FLOAT BINARYI21), 2 UB FLOAT BINARY(21) , (T1,T2) FLOAT BINARY(53) ; Tl = PRECISICM A.LB,53)**2; T2 = PRECISICN(A.UB,53)**2; . Tl = NAX(Tl,T2); RETURN(Tl); END MAXSQ; PROCU) RETURNS(FLCAT BINARY(53)); /* RETURNS MIN VALUE OF A**2 OVER INTERVAL A */ DECLARE 1 A, /* INTERVAL */ 2 LB FLOAT BINARY121) , 2 UB FLOAT B1NARY(21), (T1,T2) FLOAT BINARY153) ; IF (A.LB*A.UB <= OEOB) THEN DO; Tl = IEOB/HUGE; = PRECISION! A. LB, 53)**2; = PRECISICN(A.UB,53)**2; = MIN(T1,T2) ; RETUPM(Tl) ; END MINSO; /* INTERVAL CONTAINS ZERO ♦/ RETURN(Tl); END; Tl T2 Tl iDD_RCJNMNG„ERRCR: PROC(A); DECLARE I A LIKE SMEAKED_NUM6 ER , (D,S) FLCAT dINARY(53); /♦ COMPUTE D = MIN(MANTISSAU)) »/ /* S = 1 / (2D)**2 / 3 */ /* = COMPUTATION RCUNOOFF ERROR VARIANCE, */ /• UNITS BETA**(-2T). */ IF A. VARIANCE >= HUGE | A . COMPUT EO_VALUE=OEOB THEN RETURN; = MINIMANT ISSA(A.IA_L8) , MANT I S SA ( A. I A_UB) ) ; S = 1E03 / (C*i)*12E0) ; A. VARIANCE = kUUNQUPl A.VAK IANCE *• C); ENC ADL_kOU>iDIN'-,_tKRCK; 59 18100 18200 18300 18400 18500 18600 1870C 18800 18900 19000 19100 19200 19300 19V00 19500 19600 19700 19600 19900 20000 2010C 20200 2030C 20400 20500 20600 20700 20800 20900 21000 21100 21200 21300 21400 21500 21600 21700 21600 21900 22000 22100 22200 22300 22400 22500 22600 22/00 22800 22900 23000 23100 23200 23300 23400 23500 23600 23700 23d00 23900 24000 MANTISSA: PKOC(X) RE TURNS ( FLOAT ETNARY(53)); DLCLARE (X,Y) FLOAT bINARY(21); Y = AdS(X); 00 WHILE (Y > l.OGOOOOOOUOOOOOOOOOOOCOEOB) ; / = Y * 1.0000000000000000000000E-4B; END; DO WHILE (Y < 1.00000O0OCO0000UO00OO00E-4B) ; Y = Y * 1.0000000COOCOOOOOOOOOOOE+4B; END; RETURN(Y); END MANTISSA; SMEARED_PRINT: PROC(AA.A); DECLARE 1 A LIKE SHE ARED.NUMBER, AA CHAR(*), (S,S1,S2,ABS1,ABS2,ABS3,ABS4) FLOAT BIN(53); S = SQRTIA. VARIANCE); Si = -3.625 * S; /* 99.7% CONFIDENCE LIMITS */ S2 = +3.625 * S; ABS3 = A.COMPUTEO_VfiLUE * (1 + S1*6ETA_MT); ABS4 = A.CCMPUTEO_VALUE * (1 ♦ S2*bETA_MT); IF A.COMPUTED_VALUE < THEN DO; ABS1=ABS3; ABS3=A8S4; A6S4=ABS1; END; PUT FILE(SYSPRINT) SKIP EGIT (AA, 'COMPUTED VALUE: • , A.COMPGTEO.VALUE ) (SKIP,C0L17),A,X(8),A,E(25.15,16)) {•3.625*SIGMA BOUNDS ON RELATIVE ERROR: «,•(', SI , • ,',S2, • ) * BETA**(-T)« ) (SKIP,CGL(12),A,C0L(50),A,2 (E{25,15,16),AJJ ABS4, « )• ) (SKIP, COL 112), A, COL (50), A, 2 (E(25,15,16),A)I (•INTERVAL ANALYSIS BOUNDS:', ■ ( • , A. I A_LB , • ,', A.IA_U8,' )') (SKIP,CGL(12),A,CGL(5G),A,2 (E(16,6,7),X(9),A)); END SMEAREC.PRINT; ROUND: PROC(X) RE TURNS ( FLOAT BINARY(21J); DECLARE (X,HALF) FLOAT BINwRY(5i), FLOAT 8INARY(21), 8IT(c>4) , 3IT(8) DEFINED E PUSITICN(l), 3IT(1) DEFINED E POS IT I 0N( 33) , BIT(64) IMT( •OOOOUOUOOOOOOOOOOOO 000000 000000 01 000000 00000000 00 0000 0000 30000 GO' B) , FEXPT BITlo) DtrlNEC F POSITIGN(l); E = UNSPEC(X) ; RX E EEXPT RBIT F IF /* ROUND */ RBIT THEN DO; FEXPT = EEXPT; UNSPEC(HALF) = F; RX = X ♦ HALF; ENf ; ELSE KX = X ; /* TRUNCATE */ 60 24100 RETURN(RX); 24200 ENO RLUNO; 24300 24400 ROUNDUP: PrtOC(X) RETURNSC FLCAT BINARY(21)1; 24300 DECLARE X FLCAT bINARY(53), 24600 RX FLCAT BINARY(2ll; 24700 IF X > CEOB THEN RETURN ( RCUND (X )) ; 24800 ELSE 00; RX=X; RETURN(RX); END; 24900 END RCUNDUP; 2500C 25100 ROUNDDCJWN: PROC(X) RETURNS ( FLOAT EINARY(2D); 25200 DECLARE X FLCAT BINARY(53J, 25300 RX FLOAT 3INARY(21); 25400 IF X < OEOB THEN RETURN ( R0UNC( X )) ; 25500 ELSE DC; RX=X; FETURN(RX); END; 25600 END RCUNODOWN; 61 7. CONCLUSIONS The important conclusions of this paper are first, that accumulated roundoff errors tend to cluster towards the center of their bounds like normal densities; second, that high-confidence bounds can be proved; third, that automated statistical analysis is not worthwhile. The last conclusion is not obvious. Naturally the statistical model will generate better bounds than the worst-case model of Wilkinson. But beyond that no conclusions are obvious. Interval Analysis treats error densities as uniform, although the Central Limit Theorem tells us they are almost normal. We would expect better results from statistical analysis because it not only combines intervals (read "probability densities on intervals"), it also retains information about the re- sultant error distribution on the new interval. But statistical analysis breaks down because relative errors cannot handle cancellation or numbers near zero, and correlation is almost impossible to monitor practically. One possible exit from the irrepresentability of zero is to drop relative errors altogether and develop a statistical theory of absolute errors. Observations 2 and 3 of Section 2 provide most of the necessary theory if we now treat numbers as random variables instead of the number's relative errors. Zero then becomes representable. How- ever, it is not clear that an analogue of Theorem 5 in Section 4 can be proved for absolute analysis since multiplication and division no longer involve simple convolution, and besides, the problems of correlated errors cannot be eliminated. Statistical error analysis in general seems of questionable use. 62 REFERENCES 1. Feller, W. An Introduction to Probability Theory and Its Applica- tions v. II . NY: John Wiley & Sons, 1966. 2. Gnedenko, B. V. and A. N. Kolmogorov. Limit Distributions for Sums of Independent Random Variables . Reading, Mass.: Addi son- Wesley, 1954. 3. Hamming, R. W. "On the Distribution of Numbers" Bell Sys. Tech. J. 49, 8, 1609 (1970). 4. Henri ci, P. Elements of Numerical Analysis . NY: John Wiley & Sons, 1964. 5. Huskey, H.D. "On the Precision of a Certain Procedure of Numerical Integration." J. Res. Nat. Bur. Stds., 42, 1, 57 (1949). 6. Kaneko, T. and B. Liu. "On Local Roundoff Errors in Floating- Point Arithmetic." JACM 20, 3, 391 (1973). 7. Lever Brothers. "Vanish Advertisement." Shopper's Weekly Magazine, 17 May 1975. 8. McBride, E. B. Obtaining Generating Functions . NY: Springer- Verlag, 1971. 9. Marasa, J. D. and D. W. Matula. "A Simulative Study of Correlated Error Propagation in Various Finite-Precision Arithmetics." IEEE Trans Comp. C-22 , 6, 587 (1973). 10. Papoulis, A. The Fourier Integral and its Applications . NY: McGraw-Hill, 196^ 11. Pennington, R. H. Introductory Computer Methods and Numerical Analysis . Toronto: Macmillan, 1970. 12. Prohorov, Yu. V. and Yu. A. Rozanov. Probability Theory . NY: Springer-Verlag, 1969. 13. Sterbenz, P. H. Floating-Point Computation . Englewood Cliffs, NJ: Prentice-Hall, 1974. 14. Tsao, N. "On the Distribution of Significant Digits and Roundoff Errors." CACM V7, 5, 269 (1974). 15. Uspensky, J. V. Introduction to Mathematical Probability . NY: McGraw-Hill , 1937. 16. Wilkinson, J. H. Rounding Errors in Algebraic Processes . Englewood Cliffs, NJ: Prentice-Hall, 1963. ilBLIOGRAPHIC DATA .HEET 1. Report No. UIUCDCS-R-76-787 3. Recipient's Accession No. . Title and Subt it le THE STATISTICAL THEORY OF RELATIVE ERRORS IN FLOATING-POINT COMPUTATION 5. Report Date March 1976 6. Author(s') Douglass Stott Parker, Jr, 8. Performing Organization Rept. No UIUCDCS-R-76-787 :.iing Organization Name and Address University of Illinois at Urbana-Champaign Department of Computer Science Urbana, Illinois 61801 10. Project/Task/Work Unit N. 11. Contract/Grant No. US NSF DCR73-07980 A02 nsoring Organization Name and Address National Science Foundation Washington, D. C. 13. Type of Report & Period Covered Master's Thesis 14. t mentary Notes 6. Abstracts Experience with i most floating-poi analysis assumes cancellation and alternative to th obtain high-confi paper shows that derived, but poin error analysis in nterval analysis indicates that the error bounds it generates for nt computations are pessimistic. This is necessary since interval some error in every computation, and is oblivious of both error the Central Limit Theorem. It has been questioned whether a possible is generation of bounds is the use of error probability densities to dence intervals, or statistical bounds, on the computed error. This such a probabilistic system for Wilkinsonian relative errors can be ts out numerous crippling drawbacks for this system and statistical general . 7. Key tt'ords and Document Analysis. 17a. Descriptors Automated error analysis Floating-point computation Interval analysis Roundoff error Statistical error bounds 7b. Identifiers 'Open-Ended Terms 7c. COSATI Fie Id /Group ability Statement Release Unlimited 19. Security Class (This Report) UNCLASSIFIED 20. Security ("lass (This Page UNCLASSIFIED 21. No. of Pages 66 22. Price "-"" NTrS-35 (10-70) USCOMM-DC 40329-P71 ;>•■ Vf*** UNIVERSITY OF ILLINOIS URBANA 510 84 IL6R no C002 no 782 787(1976 Design ol Irredundant MOS Networks a p 0112 088402588