home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The C Users' Group Library 1994 August
/
wc-cdrom-cusersgrouplibrary-1994-08.iso
/
listings
/
v_02_12
/
371.att
/
TSMITH.LTR
Wrap
Text File
|
1991-06-06
|
4KB
|
103 lines
The following letter is for consideration of publication in
TECH Specialist.
----------------------------------------------------------------
Dear Mr. Ward:
I have several comments on "The Mathematical Programmer, Part 3"
by Scott Ladd (June 91).
The tonearest() function is the most accurate and complete
integer rounding routine in C that I have seen. Readers should
remember to insert #include <math.h> at the top.
The concept of using an allowable tolerance (epsilon) for
terminating a converging iteration is good, but the magnitude of
the numbers being compared must be considered. The ANSI-
specified value recommended in the article, FLT_EPSILON,
represents the smallest possible change to 1.0F, not the smallest
possible difference between two float values as stated.
For solution values over a wide range, the proper value for
epsilon should normally vary with the magnitude of the numbers
being compared to achieve a given number of significant digits.
The relative difference can be determined by dividing the
difference by the old or new value of x before comparing, e.g.:
fabs((xnew-xold)/xnew) <= tolerance.
For some functions an absolute epsilon should be used for values
of x near zero, resulting in two terminating conditions:
#include <math.h>
#define ABSOLUTE_EPSILON 1e-10
#define RELATIVE_EPSILON 1e-6 /* 6 digits */
float xnew, xold;
xnew = 1.0; /* some initial estimate */
do {
xold = xnew;
xnew = solution_refinement (xold);
} while ( fabs(xnew - xold) > ABSOLUTE_EPSILON &&
fabs (xnew-xold) > fabs(xnew) * RELATIVE_EPSILON );
I cannot agree with the author's recommendation of limited-
precision rounding during intermediate calculations. Each
instance of such rounding adds more error to the final result. A
more accurate estimate is given by carrying full precision
throughout the calculations, then rounding to the desired
precision on presentation (usually by a conversion format
specification). The limited-precision rounding can indicate an
estimate of accuracy, but not enhance it.
The randomd() function, written to return a random double value,
converts from an int to a double _after_ dividing, resulting in a
return value of 0 for all calls (except when rand() is RAND_MAX)!
Assuming that a number is to be generated in the interval from
0.0 up to but excluding 1.0 (the normally desired case), the
function should read
#include <stdlib.h>
double randomd (void) {
return (rand() / (RAND_MAX + 1.0)) ;
}
Note the use of RAND_MAX, defined to be maximum value returned by
rand().
Readers should be aware than the technique of setting a pointer
to the address of an array minus one element to use for one-based
indexing results in technically undefined behavior in ANSI C,
although it is probably safe for most MS-DOS compilers. Another
way, which is defined under ANSI C, is to define a macro for the
indexing:
#define fap(i) fa[(i)-1]
The programmer then uses parentheses rather than brackets for the
subscripts (this may be an advantage for translating Fortran
code!).
The IEEE 754 short real format for floating point has 23
signficand bits plus an implied leading 1 bit, resulting in 24
bits of precision. This is equivalent to 24 * log(2),
approximately 7.2, decimal digits, rather than 6 as stated in the
article. Microsoft's value of 7 for FLT_DIG is thus correct.
This can be demonstrated by reading values into a float (on a
machine using IEEE short real format), then printing them. For
each 7-digit test case that I tried (not all 9 million!) I got
the exact number back.
Sincerely,
Thad Smith III
T3 Systems
2001 N. 75th St.
Boulder, CO 80301
For editor's use only:
(303)449-3322 (voice)
Fidonet: 1:104/89.3