home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / X / mit / lib / X / XcmsMath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-15  |  3.2 KB  |  134 lines
  1. /*
  2.  * square and cube roots by Newton's method
  3.  *
  4.  * $XConsortium: XcmsMath.c,v 1.8 91/08/14 15:13:41 rws Exp $
  5.  *
  6.  * Copyright 1990 Massachusetts Institute of Technology
  7.  *
  8.  * Permission to use, copy, modify, distribute, and sell this software and its
  9.  * documentation for any purpose is hereby granted without fee, provided that
  10.  * the above copyright notice appear in all copies and that both that
  11.  * copyright notice and this permission notice appear in supporting
  12.  * documentation, and that the name of M.I.T. not be used in advertising or
  13.  * publicity pertaining to distribution of the software without specific,
  14.  * written prior permission.  M.I.T. makes no representations about the
  15.  * suitability of this software for any purpose.  It is provided "as is"
  16.  * without express or implied warranty.
  17.  *
  18.  * M.I.T. DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL
  19.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL M.I.T.
  20.  * BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  21.  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
  22.  * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN 
  23.  * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  24.  *
  25.  * Stephen Gildea, MIT X Consortium, January 1991
  26.  */
  27.  
  28. #include "Xlibint.h"
  29. #include "Xcmsint.h"
  30.  
  31. #if !defined(X_NOT_STDC_ENV) && (__STDC__ || !(defined(sun) || (defined(sony) && !defined(SYSTYPE_SYSV))))
  32. #include <float.h>
  33. #endif
  34. #ifndef DBL_EPSILON
  35. #define DBL_EPSILON 1e-6
  36. #endif
  37.  
  38. #ifdef _X_ROOT_STATS
  39. int cbrt_loopcount;
  40. int sqrt_loopcount;
  41. #endif
  42.  
  43. /* Newton's Method:  x_n+1 = x_n - ( f(x_n) / f'(x_n) ) */
  44.  
  45.  
  46. /* for cube roots, x^3 - a = 0,  x_new = x - 1/3 (x - a/x^2) */
  47.  
  48. double
  49. _XcmsCubeRoot(a)
  50.     double a;
  51. {
  52.     register double abs_a, cur_guess, delta;
  53.  
  54. #ifdef DEBUG
  55.     printf("_XcmsCubeRoot passed in %g\n", a);
  56. #endif
  57. #ifdef _X_ROOT_STATS
  58.     cbrt_loopcount = 0;
  59. #endif
  60.     if (a == 0.)
  61.     return 0.;
  62.  
  63.     abs_a = a<0. ? -a : a;    /* convert to positive to speed loop tests */
  64.  
  65.     /* arbitrary first guess */
  66.     if (abs_a > 1.)
  67.     cur_guess = abs_a/8.;
  68.     else
  69.     cur_guess = abs_a*8.;
  70.  
  71.     do {
  72. #ifdef _X_ROOT_STATS
  73.     cbrt_loopcount++;
  74. #endif
  75.     delta = (cur_guess - abs_a/(cur_guess*cur_guess))/3.;
  76.     cur_guess -= delta;
  77.     if (delta < 0.) delta = -delta;
  78.     } while (delta >= cur_guess*DBL_EPSILON);
  79.  
  80.     if (a < 0.)
  81.     cur_guess = -cur_guess;
  82.  
  83. #ifdef DEBUG
  84.     printf("_XcmsCubeRoot returning %g\n", cur_guess);
  85. #endif
  86.     return cur_guess;
  87. }
  88.     
  89.  
  90.  
  91. /* for square roots, x^2 - a = 0,  x_new = x - 1/2 (x - a/x) */
  92.  
  93. double
  94. _XcmsSquareRoot(a)
  95.     double a;
  96. {
  97.     register double cur_guess, delta;
  98.  
  99. #ifdef DEBUG
  100.     printf("_XcmsSquareRoot passed in %g\n", a);
  101. #endif
  102. #ifdef _X_ROOT_STATS
  103.     sqrt_loopcount = 0;
  104. #endif
  105.     if (a == 0.)
  106.     return 0.;
  107.  
  108.     if (a < 0.) {
  109.     /* errno = EDOM; */
  110.     return 0.;
  111.     }
  112.  
  113.     /* arbitrary first guess */
  114.     if (a > 1.)
  115.     cur_guess = a/4.;
  116.     else
  117.     cur_guess = a*4.;
  118.  
  119.     do {
  120. #ifdef _X_ROOT_STATS
  121.     sqrt_loopcount++;
  122. #endif
  123.     delta = (cur_guess - a/cur_guess)/2.;
  124.     cur_guess -= delta;
  125.     if (delta < 0.) delta = -delta;
  126.     } while (delta >= cur_guess*DBL_EPSILON);
  127.  
  128. #ifdef DEBUG
  129.     printf("_XcmsSquareRoot returning %g\n", cur_guess);
  130. #endif
  131.     return cur_guess;
  132. }
  133.     
  134.