4.3BSD - /usr/src/usr.lib/libm/IEEE/cabs.c

   1: /*
   2:  * Copyright (c) 1985 Regents of the University of California.
   3:  *
   4:  * Use and reproduction of this software are granted  in  accordance  with
   5:  * the terms and conditions specified in  the  Berkeley  Software  License
   6:  * Agreement (in particular, this entails acknowledgement of the programs'
   7:  * source, and inclusion of this notice) with the additional understanding
   8:  * that  all  recipients  should regard themselves as participants  in  an
   9:  * ongoing  research  project and hence should  feel  obligated  to report
  10:  * their  experiences (good or bad) with these elementary function  codes,
  11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
  12:  */
  13: 
  14: #ifndef lint
  15: static char sccsid[] = "@(#)cabs.c	1.2 (Berkeley) 8/21/85";
  16: #endif not lint
  17: 
  18: /* CABS(Z)
  19:  * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER  Z = X + iY
  20:  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  21:  * CODED IN C BY K.C. NG, 11/28/84.
  22:  * REVISED BY K.C. NG, 7/12/85.
  23:  *
  24:  * Required kernel function :
  25:  *	hypot(x,y)
  26:  *
  27:  * Method :
  28:  *	cabs(z) = hypot(x,y) .
  29:  */
  30: 
  31: double cabs(z)
  32: struct { double x, y;} z;
  33: {
  34:     double hypot();
  35:     return(hypot(z.x,z.y));
  36: }
  37: 
  38: 
  39: /* HYPOT(X,Y)
  40:  * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
  41:  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  42:  * CODED IN C BY K.C. NG, 11/28/84;
  43:  * REVISED BY K.C. NG, 7/12/85.
  44:  *
  45:  * Required system supported functions :
  46:  *	copysign(x,y)
  47:  *	finite(x)
  48:  *	scalb(x,N)
  49:  *	sqrt(x)
  50:  *
  51:  * Method :
  52:  *	1. replace x by |x| and y by |y|, and swap x and
  53:  *	   y if y > x (hence x is never smaller than y).
  54:  *	2. Hypot(x,y) is computed by:
  55:  *	   Case I, x/y > 2
  56:  *
  57:  *				       y
  58:  *		hypot = x + -----------------------------
  59:  *			 		    2
  60:  *			    sqrt ( 1 + [x/y]  )  +  x/y
  61:  *
  62:  *	   Case II, x/y <= 2
  63:  *				                   y
  64:  *		hypot = x + --------------------------------------------------
  65:  *				          		     2
  66:  *				     			[x/y]   -  2
  67:  *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
  68:  *			 		    			  2
  69:  *			    			  sqrt ( 1 + [x/y]  )  + sqrt(2)
  70:  *
  71:  *
  72:  *
  73:  * Special cases:
  74:  *	hypot(x,y) is INF if x or y is +INF or -INF; else
  75:  *	hypot(x,y) is NAN if x or y is NAN.
  76:  *
  77:  * Accuracy:
  78:  * 	hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
  79:  *	in the last place). See Kahan's "Interval Arithmetic Options in the
  80:  *	Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
  81:  *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
  82:  *	code follows in	comments.) In a test run with 500,000 random arguments
  83:  *	on a VAX, the maximum observed error was .959 ulps.
  84:  *
  85:  * Constants:
  86:  * The hexadecimal values are the intended ones for the following constants.
  87:  * The decimal values may be used, provided that the compiler will convert
  88:  * from decimal to binary accurately enough to produce the hexadecimal values
  89:  * shown.
  90:  */
  91: 
  92: #ifdef VAX  /* VAX D format */
  93: /* static double */
  94: /* r2p1hi =  2.4142135623730950345E0     , Hex  2^  2   *  .9A827999FCEF32 */
  95: /* r2p1lo =  1.4349369327986523769E-17   , Hex  2^-55   *  .84597D89B3754B */
  96: /* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
  97: static long    r2p1hix[] = { 0x8279411a, 0xef3299fc};
  98: static long    r2p1lox[] = { 0x597d2484, 0x754b89b3};
  99: static long     sqrt2x[] = { 0x04f340b5, 0xde6533f9};
 100: #define   r2p1hi    (*(double*)r2p1hix)
 101: #define   r2p1lo    (*(double*)r2p1lox)
 102: #define    sqrt2    (*(double*)sqrt2x)
 103: #else       /* IEEE double format */
 104: static double
 105: r2p1hi =  2.4142135623730949234E0     , /*Hex  2^1     *  1.3504F333F9DE6 */
 106: r2p1lo =  1.2537167179050217666E-16   , /*Hex  2^-53   *  1.21165F626CDD5 */
 107: sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
 108: #endif
 109: 
 110: double hypot(x,y)
 111: double x, y;
 112: {
 113:     static double zero=0, one=1,
 114:               small=1.0E-18;    /* fl(1+small)==1 */
 115:     static ibig=30; /* fl(1+2**(2*ibig))==1 */
 116:     double copysign(),scalb(),logb(),sqrt(),t,r;
 117:     int finite(), exp;
 118: 
 119:     if(finite(x))
 120:         if(finite(y))
 121:         {
 122:         x=copysign(x,one);
 123:         y=copysign(y,one);
 124:         if(y > x)
 125:             { t=x; x=y; y=t; }
 126:         if(x == zero) return(zero);
 127:         if(y == zero) return(x);
 128:         exp= logb(x);
 129:         if(exp-(int)logb(y) > ibig )
 130:             /* raise inexact flag and return |x| */
 131:            { one+small; return(x); }
 132: 
 133:         /* start computing sqrt(x^2 + y^2) */
 134:         r=x-y;
 135:         if(r>y) {   /* x/y > 2 */
 136:             r=x/y;
 137:             r=r+sqrt(one+r*r); }
 138:         else {      /* 1 <= x/y <= 2 */
 139:             r/=y; t=r*(r+2.0);
 140:             r+=t/(sqrt2+sqrt(2.0+t));
 141:             r+=r2p1lo; r+=r2p1hi; }
 142: 
 143:         r=y/r;
 144:         return(x+r);
 145: 
 146:         }
 147: 
 148:         else if(y==y)          /* y is +-INF */
 149:              return(copysign(y,one));
 150:         else
 151:              return(y);    /* y is NaN and x is finite */
 152: 
 153:     else if(x==x)          /* x is +-INF */
 154:              return (copysign(x,one));
 155:     else if(finite(y))
 156:              return(x);        /* x is NaN, y is finite */
 157:     else if(y!=y) return(y);  /* x and y is NaN */
 158:     else return(copysign(y,one));   /* y is INF */
 159: }
 160: 
 161: /* A faster but less accurate version of cabs(x,y) */
 162: #if 0
 163: double hypot(x,y)
 164: double x, y;
 165: {
 166:     static double zero=0, one=1;
 167:               small=1.0E-18;    /* fl(1+small)==1 */
 168:     static ibig=30; /* fl(1+2**(2*ibig))==1 */
 169:     double copysign(),scalb(),logb(),sqrt(),temp;
 170:     int finite(), exp;
 171: 
 172:     if(finite(x))
 173:         if(finite(y))
 174:         {
 175:         x=copysign(x,one);
 176:         y=copysign(y,one);
 177:         if(y > x)
 178:             { temp=x; x=y; y=temp; }
 179:         if(x == zero) return(zero);
 180:         if(y == zero) return(x);
 181:         exp= logb(x);
 182:         x=scalb(x,-exp);
 183:         if(exp-(int)logb(y) > ibig )
 184:             /* raise inexact flag and return |x| */
 185:            { one+small; return(scalb(x,exp)); }
 186:         else y=scalb(y,-exp);
 187:         return(scalb(sqrt(x*x+y*y),exp));
 188:         }
 189: 
 190:         else if(y==y)          /* y is +-INF */
 191:              return(copysign(y,one));
 192:         else
 193:              return(y);    /* y is NaN and x is finite */
 194: 
 195:     else if(x==x)          /* x is +-INF */
 196:              return (copysign(x,one));
 197:     else if(finite(y))
 198:              return(x);        /* x is NaN, y is finite */
 199:     else if(y!=y) return(y);    /* x and y is NaN */
 200:     else return(copysign(y,one));   /* y is INF */
 201: }
 202: #endif

Defined functions

Defined variables

Defined macros