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
Text File  |  1991-06-06  |  4KB  |  103 lines
  1. The following letter is for consideration of publication in
  2. TECH Specialist.
  3. ----------------------------------------------------------------
  4. Dear Mr. Ward:
  5.  
  6. I have several comments on "The Mathematical Programmer, Part 3"
  7. by Scott Ladd (June 91).  
  8.  
  9. The tonearest() function is the most accurate and complete
  10. integer rounding routine in C that I have seen.  Readers should
  11. remember to insert #include <math.h> at the top.
  12.  
  13. The concept of using an allowable tolerance (epsilon) for
  14. terminating a converging iteration is good, but the magnitude of
  15. the numbers being compared must be considered.  The ANSI-
  16. specified value recommended in the article, FLT_EPSILON,
  17. represents the smallest possible change to 1.0F, not the smallest
  18. possible difference between two float values as stated.
  19.  
  20. For solution values over a wide range, the proper value for
  21. epsilon should normally vary with the magnitude of the numbers
  22. being compared to achieve a given number of significant digits. 
  23. The relative difference can be determined by dividing the
  24. difference by the old or new value of x before comparing, e.g.:
  25.  
  26.   fabs((xnew-xold)/xnew) <= tolerance.
  27.  
  28. For some functions an absolute epsilon should be used for values
  29. of x near zero, resulting in two terminating conditions:
  30.  
  31. #include <math.h>
  32. #define ABSOLUTE_EPSILON      1e-10
  33. #define RELATIVE_EPSILON      1e-6   /* 6 digits */
  34. float xnew, xold;
  35.  
  36. xnew = 1.0;  /* some initial estimate */
  37. do { 
  38.   xold = xnew;
  39.   xnew = solution_refinement (xold);
  40. } while ( fabs(xnew - xold) > ABSOLUTE_EPSILON &&
  41.      fabs (xnew-xold) > fabs(xnew) * RELATIVE_EPSILON );
  42.  
  43. I cannot agree with the author's recommendation of limited-
  44. precision rounding during intermediate calculations.  Each
  45. instance of such rounding adds more error to the final result.  A
  46. more accurate estimate is given by carrying full precision
  47. throughout the calculations, then rounding to the desired
  48. precision on presentation (usually by a conversion format
  49. specification).  The limited-precision rounding can indicate an
  50. estimate of accuracy, but not enhance it.
  51.  
  52. The randomd() function, written to return a random double value,
  53. converts from an int to a double _after_ dividing, resulting in a
  54. return value of 0 for all calls (except when rand() is RAND_MAX)!
  55.  
  56. Assuming that a number is to be generated in the interval from
  57. 0.0 up to but excluding 1.0 (the normally desired case), the
  58. function should read
  59.  
  60. #include <stdlib.h>
  61. double randomd (void) {
  62.   return (rand() / (RAND_MAX + 1.0)) ;
  63. }
  64.  
  65. Note the use of RAND_MAX, defined to be maximum value returned by
  66. rand().  
  67.  
  68. Readers should be aware than the technique of setting a pointer
  69. to the address of an array minus one element to use for one-based
  70. indexing results in technically undefined behavior in ANSI C,
  71. although it is probably safe for most MS-DOS compilers.  Another
  72. way, which is defined under ANSI C, is to define a macro for the
  73. indexing:
  74.  
  75. #define fap(i)  fa[(i)-1]
  76.  
  77. The programmer then uses parentheses rather than brackets for the
  78. subscripts (this may be an advantage for translating Fortran
  79. code!).
  80.  
  81. The IEEE 754 short real format for floating point has 23
  82. signficand bits plus an implied leading 1 bit, resulting in 24
  83. bits of precision.  This is equivalent to 24 * log(2),
  84. approximately 7.2, decimal digits, rather than 6 as stated in the
  85. article.  Microsoft's value of 7 for FLT_DIG is thus correct. 
  86. This can be demonstrated by reading values into a float (on a
  87. machine using IEEE short real format), then printing them.  For
  88. each 7-digit test case that I tried (not all 9 million!) I got
  89. the exact number back.
  90.  
  91. Sincerely,
  92.  
  93.  
  94. Thad Smith III
  95. T3 Systems
  96. 2001 N. 75th St.
  97. Boulder, CO 80301
  98.  
  99.  
  100. For editor's use only:
  101. (303)449-3322 (voice)
  102. Fidonet: 1:104/89.3
  103.