🚀 go-pugleaf

In article <1992Aug6.085315.10239@mic.ucla.edu> iwelch@agsm.ucla.edu
(Ivo Welch) writes:
> [ referring to the "Numerical Recipes in C" SVDCMP function ]
>Furthermore, perhaps someone can enlighten me as to what the statement
>
>        if (fabs(rv1[l])+anorm == anorm)
>
>tests. anorm is the max of abs(w[]) + abs(rv1[]).

This is actually a nice way of testing whether or not  rv1[l]  is
"negligably small" relative to  w[] .  That is, in general the test

	if (x == (x + e))

(the  fabs()  actually isn't needed) tests whether or not  e  is
negligably small relative to  x .  More precisely, to within at
most a factor of two (which doesn't matter here) it tests for
 |e/x| < eps  where  eps  is the "machine epsilon", i.e. the
smallest positive floating point number such that  1.0+eps != 1.0
in floating point.  (eps = 5.97e-8 for IEEE single, 1.12e-16 for
IEEE double precision.)  This can be useful in testing whether or
not some iteration which updates  x  by adding  e  to it, has
converged.

This is a *nice* way of testing this because it *automatically*
adapts to the floating point format used.  That is, code written
using this style of test should work unchanged with any reasonable
floating point system (even Cray) *without* requiring changes to
"machine constants".

The same trick can also be used to test for things being "almost
negligable":

	if (x == (x + 0.01*e))

tests for  x + e  and  x  differing by at most 100 "units in the
last place" (100 times the value of the least significant bit in
 x ), i.e. for  |e/x| < 100*eps .  Here again the error tolerance
*automatically* adapts to the floating point precision.

(The one caveat in this is that the usual floating point precision
must be used for the computation/storage of  (x + e) , or at least
this quantity must be stored in a usual-precision register or memory
location.  If this computation is done/stored in (say) an extended-
precision register, all bets are off.  The usual hack to force
 (x + e)  to be computed in the "usual precision" is (using Fortran)

		double precision x, e
		double precision temp(2)
		common /extemp/ temp
		...
	c
	c	is  e  negligably small compared to  x ?
		temp(1) = x
		temp(2) = temp(1) + e
		if (temp(1) .eq. temp(2)) then
			... e  is negligably small compared to  x
		else
			... it isn't
		end if

(In C the array would be declared with "extern" storage class, i.e.
it would be declared at "top level", not inside a function, and would
*not* be declared "static".)

Here the use of an array, and its placement in a common block, are
intended to force the compiler to actually compute and store  x + e
in memory, hopefully rounding it to the "usual precision".  This is
rather ugly, but is safe against most optimizing compilers and
extended-precision-registers hardware.

- Jonathan Thornburg
  <jonathan@einstein.ph.utexas.edu> or <jonathan@hermes.chpc.utexas.edu>
  University of Texas at Austin / Physics Dept / Center for Relativity
  and (for a few more months) U of British Columbia / {Astronomy,Physics}