2.9BSD - /usr/src/lib/m/sqrt.c

   1: /*	sqrt - floating-point square root */
   2: /*	Optimized for PDP-11 */
   3: /*	by Bruce R. Julian, U.S. Geological Survey, Menlo Park, CA    22 Nov 1982*/
   4: 
   5: /*	calls frexp */
   6: 
   7: #include <whoami.h>
   8: #include <errno.h>
   9: #define SQRT2 1.4142135623730950488
  10: 
  11: int errno;
  12: double frexp();
  13: 
  14: double
  15: sqrt(arg)
  16: double arg;
  17: {
  18:     double x, temp;
  19:     int exp;
  20: 
  21:     /* Test for negative or zero argument */
  22:     if(arg <= 0.) {
  23:         if(arg < 0.)
  24:             errno = EDOM;
  25:         return(0.);
  26:     }
  27: 
  28:     /* Split into fraction and exponent and estimate sqrt(fraction) */
  29:     x = frexp(arg,&exp);
  30: 
  31: #ifndef PDP11
  32:     /* Insure fraction is in range (.5,1] */
  33:     /* (Necessary only with non-binary floating point) */
  34:     while (x < 0.5) {
  35:         x *= 2.;
  36:         exp--;
  37:     }
  38: #endif
  39: 
  40:     temp = 0.59016*x + 0.41731;     /* Minimax relative error on [.5, 1] */
  41: 
  42:     /* Make exponent even */
  43:     /* NOTE: This won't work on 1's complement machines */
  44:     if(exp & 1) {               /* Exponent odd? */
  45:         temp *= SQRT2;
  46:         exp--;
  47:     }
  48: 
  49:     /* Get exponent into range [-60, 60] */
  50:     while(exp > 60) {
  51:         temp *= (1L<<30);
  52:         exp -= 60;
  53:     }
  54:     while(exp < -60) {
  55:         temp /= (1L<<30);
  56:         exp += 60;
  57:     }
  58: 
  59:     /* Multiply temp by sqrt(2**exp) */
  60:     if(exp >= 0)
  61:         temp *= 1L << (exp/2);
  62:     else
  63:         temp /= 1L << (-exp/2);
  64: 
  65:     /* Newton's method (3 iterations) */
  66:     temp = 0.5*(temp + arg/temp);
  67:     temp = temp + arg/temp;
  68:     return(0.25*temp + arg/temp);
  69: }

Defined functions

Defined variables

Defined macros