home *** CD-ROM | disk | FTP | other *** search
/ SPACE 2 / SPACE - Library 2 - Volume 1.iso / program / 167 / pascal / expbugs.doc < prev    next >
Encoding:
Text File  |  1987-08-19  |  3.9 KB  |  126 lines
  1. 03/22/86    [ using V1.11]
  2.  
  3.                 Inaccuracies in the Ln and Exp functions
  4.                 and some ways around them.
  5.  
  6. While writing some financial calculation routines, I needed a power function
  7. to raise a number to any power.  Since Pascal does not have a built-in power
  8. function, I wrote a very simple one:
  9.  
  10.                       --------------------------------
  11.  
  12.         FUNTION Power(base, exponent : Real ) : Real ;
  13.  
  14.         {* This function merely raises the base to a power specified by the
  15.          * exponent using the relation
  16.          *                                y
  17.          *                            Ln X  = y * Ln X
  18.          *
  19.          *                                y    y * Ln X
  20.          *                               X  = e                         }
  21.  
  22.  
  23.                 BEGIN
  24.                      Power := Exp( exponent * Ln(base) ) ;
  25.                 END;
  26.  
  27.                       ---------------------------------
  28.  
  29.  
  30. However, by comparing the answers I was getting to standard financial tables, I
  31. noticed that the answers I was getting were only accurate to 5 or 6 digits.
  32.  
  33.         I then wrote a simple routine to display ln(x) and exp(x) for various
  34. values of x and found that the routines were accurate for values of x near one,
  35. but became increasingly less accurate for larger values.
  36.  
  37.         I then wrote the following two functions which can be used in place of
  38. the built in exp and ln functions.  Neither is extremely elegant or fast for
  39. all arguments, but they give the right answers.
  40.  
  41. -----------------------------------------------------
  42.  
  43. FUNCTION  Exp_new ( exponent : Real ) : Real ;
  44.  
  45. {* The exp and ln functions furnished with Personal Pascal are not accurate
  46.  * for large ranges of arguments.  This function calculates exp(x) using a
  47.  * Taylor Series. If the exponent is less than zero, then the absolute value
  48.  * of the exponent is used in the calculation and the result inverted. }
  49.  
  50. CONST   Error  =  5.0E-12 ;
  51.  
  52. VAR     Total, Term, Dummy                           : Real        ;
  53.         Count                                        : Integer     ;
  54.  
  55. Begin
  56.  
  57.  
  58.     IF (exponent < 0) THEN
  59.     Dummy := -1 * exponent
  60.     Else  Dummy := exponent  ;
  61.  
  62.         Total := 0 ;
  63.         Count := 0 ;
  64.         Term  := 1 ;
  65.  
  66.         Repeat
  67.                 Count := Count + 1 ;
  68.                 Total := Total + Term ;
  69.                 Term  := (Term * Dummy) / Count ;
  70.         Until ( Abs(Term) < Total * Error ) ;
  71.  
  72.     IF ( exponent >= 0 ) THEN  Exp_new := Total
  73.     Else Exp_new := 1/Total ;
  74.  
  75. End ;
  76.  
  77. ------------------------------------------------------------
  78.  
  79.  
  80. FUNCTION  Ln_new ( arg : Real ) : Real ;
  81.  
  82. {* The exp and ln functions furnished with Personal Pascal are not accurate
  83.  * for large ranges of arguments.  This function calculates ln(x) using
  84.  * Newton's iterative method to solve the equation:
  85.  *                       x
  86.  *                      e  - arg = 0
  87.  *
  88.  * Notice that this function has to call exp_new.  }
  89.  
  90. CONST   Error = 5.0E-12 ;
  91.  
  92. VAR     X1, X2 : Real ;
  93.  
  94. BEGIN
  95.  
  96. { * The first IF-THEN takes care of the special case when arg=1 }
  97.         IF (arg = 1) THEN
  98.         Ln_new := 0
  99.         ELSE
  100.         BEGIN
  101.  
  102. { * The next IF-THEN chooses a more favorable starting point for iterations for
  103.   * large values of arg. [ You could improve on this section if you want more
  104.   * speed for large values of arg.]
  105.  
  106.                 IF (arg > 10) THEN
  107.                 X2 := arg/10
  108.                 ELSE
  109.                 X2 := arg ;
  110.  
  111. { * Now solve the equation. }
  112.                 Repeat
  113.                         X1 := X2 ;
  114.                         X2 := X1 - 1 + arg * Exp_new( -1 * X1 ) ;
  115.                 Until ( ABS((X2-X1)/X1) <= Error) ;
  116.                 Ln_new := X2 ;
  117.         END ;
  118.  
  119.  
  120. END ;
  121.  
  122. ------------------------------------------------
  123. MIKE GUIDRY
  124. CIS: 71046,144
  125. I also hang out on the Seattle and New Jersey ST BBS's.
  126.