![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to cabs.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / ieee|
Return to cabs.c CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / ieee |
1.1 root 1: /*
2: * Copyright (c) 1985 Regents of the University of California.
3: * All rights reserved.
4: *
5: * Redistribution and use in source and binary forms are permitted
6: * provided that the above copyright notice and this paragraph are
7: * duplicated in all such forms and that any documentation,
8: * advertising materials, and other materials related to such
9: * distribution and use acknowledge that the software was developed
10: * by the University of California, Berkeley. The name of the
11: * University may not be used to endorse or promote products derived
12: * from this software without specific prior written permission.
13: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15: * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16: *
17: * All recipients should regard themselves as participants in an ongoing
18: * research project and hence should feel obligated to report their
19: * experiences (good or bad) with these elementary function codes, using
20: * the sendbug(8) program, to the authors.
21: */
22:
23: #ifndef lint
24: static char sccsid[] = "@(#)cabs.c 5.3 (Berkeley) 6/30/88";
25: #endif /* not lint */
26:
27: /* HYPOT(X,Y)
28: * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY
29: * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
30: * CODED IN C BY K.C. NG, 11/28/84;
31: * REVISED BY K.C. NG, 7/12/85.
32: *
33: * Required system supported functions :
34: * copysign(x,y)
35: * finite(x)
36: * scalb(x,N)
37: * sqrt(x)
38: *
39: * Method :
40: * 1. replace x by |x| and y by |y|, and swap x and
41: * y if y > x (hence x is never smaller than y).
42: * 2. Hypot(x,y) is computed by:
43: * Case I, x/y > 2
44: *
45: * y
46: * hypot = x + -----------------------------
47: * 2
48: * sqrt ( 1 + [x/y] ) + x/y
49: *
50: * Case II, x/y <= 2
51: * y
52: * hypot = x + --------------------------------------------------
53: * 2
54: * [x/y] - 2
55: * (sqrt(2)+1) + (x-y)/y + -----------------------------
56: * 2
57: * sqrt ( 1 + [x/y] ) + sqrt(2)
58: *
59: *
60: *
61: * Special cases:
62: * hypot(x,y) is INF if x or y is +INF or -INF; else
63: * hypot(x,y) is NAN if x or y is NAN.
64: *
65: * Accuracy:
66: * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
67: * in the last place). See Kahan's "Interval Arithmetic Options in the
68: * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
69: * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
70: * code follows in comments.) In a test run with 500,000 random arguments
71: * on a VAX, the maximum observed error was .959 ulps.
72: *
73: * Constants:
74: * The hexadecimal values are the intended ones for the following constants.
75: * The decimal values may be used, provided that the compiler will convert
76: * from decimal to binary accurately enough to produce the hexadecimal values
77: * shown.
78: */
79:
80: #if defined(vax)||defined(tahoe) /* VAX D format */
81: #ifdef vax
82: #define _0x(A,B) 0x/**/A/**/B
83: #else /* vax */
84: #define _0x(A,B) 0x/**/B/**/A
85: #endif /* vax */
86: /* static double */
87: /* r2p1hi = 2.4142135623730950345E0 , Hex 2^ 2 * .9A827999FCEF32 */
88: /* r2p1lo = 1.4349369327986523769E-17 , Hex 2^-55 * .84597D89B3754B */
89: /* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
90: static long r2p1hix[] = { _0x(8279,411a), _0x(ef32,99fc)};
91: static long r2p1lox[] = { _0x(597d,2484), _0x(754b,89b3)};
92: static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)};
93: #define r2p1hi (*(double*)r2p1hix)
94: #define r2p1lo (*(double*)r2p1lox)
95: #define sqrt2 (*(double*)sqrt2x)
96: #else /* defined(vax)||defined(tahoe) */
97: static double
98: r2p1hi = 2.4142135623730949234E0 , /*Hex 2^1 * 1.3504F333F9DE6 */
99: r2p1lo = 1.2537167179050217666E-16 , /*Hex 2^-53 * 1.21165F626CDD5 */
100: sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */
101: #endif /* defined(vax)||defined(tahoe) */
102:
103: double
104: hypot(x,y)
105: double x, y;
106: {
107: static double zero=0, one=1,
108: small=1.0E-18; /* fl(1+small)==1 */
109: static ibig=30; /* fl(1+2**(2*ibig))==1 */
110: double copysign(),scalb(),logb(),sqrt(),t,r;
111: int finite(), exp;
112:
113: if(finite(x))
114: if(finite(y))
115: {
116: x=copysign(x,one);
117: y=copysign(y,one);
118: if(y > x)
119: { t=x; x=y; y=t; }
120: if(x == zero) return(zero);
121: if(y == zero) return(x);
122: exp= logb(x);
123: if(exp-(int)logb(y) > ibig )
124: /* raise inexact flag and return |x| */
125: { one+small; return(x); }
126:
127: /* start computing sqrt(x^2 + y^2) */
128: r=x-y;
129: if(r>y) { /* x/y > 2 */
130: r=x/y;
131: r=r+sqrt(one+r*r); }
132: else { /* 1 <= x/y <= 2 */
133: r/=y; t=r*(r+2.0);
134: r+=t/(sqrt2+sqrt(2.0+t));
135: r+=r2p1lo; r+=r2p1hi; }
136:
137: r=y/r;
138: return(x+r);
139:
140: }
141:
142: else if(y==y) /* y is +-INF */
143: return(copysign(y,one));
144: else
145: return(y); /* y is NaN and x is finite */
146:
147: else if(x==x) /* x is +-INF */
148: return (copysign(x,one));
149: else if(finite(y))
150: return(x); /* x is NaN, y is finite */
151: #if !defined(vax)&&!defined(tahoe)
152: else if(y!=y) return(y); /* x and y is NaN */
153: #endif /* !defined(vax)&&!defined(tahoe) */
154: else return(copysign(y,one)); /* y is INF */
155: }
156:
157: /* CABS(Z)
158: * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY
159: * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
160: * CODED IN C BY K.C. NG, 11/28/84.
161: * REVISED BY K.C. NG, 7/12/85.
162: *
163: * Required kernel function :
164: * hypot(x,y)
165: *
166: * Method :
167: * cabs(z) = hypot(x,y) .
168: */
169:
170: double
171: cabs(z)
172: struct { double x, y;} z;
173: {
174: return hypot(z.x,z.y);
175: }
176:
177: double
178: z_abs(z)
179: struct { double x,y;} *z;
180: {
181: return hypot(z->x,z->y);
182: }
183:
184: /* A faster but less accurate version of cabs(x,y) */
185: #if 0
186: double hypot(x,y)
187: double x, y;
188: {
189: static double zero=0, one=1;
190: small=1.0E-18; /* fl(1+small)==1 */
191: static ibig=30; /* fl(1+2**(2*ibig))==1 */
192: double copysign(),scalb(),logb(),sqrt(),temp;
193: int finite(), exp;
194:
195: if(finite(x))
196: if(finite(y))
197: {
198: x=copysign(x,one);
199: y=copysign(y,one);
200: if(y > x)
201: { temp=x; x=y; y=temp; }
202: if(x == zero) return(zero);
203: if(y == zero) return(x);
204: exp= logb(x);
205: x=scalb(x,-exp);
206: if(exp-(int)logb(y) > ibig )
207: /* raise inexact flag and return |x| */
208: { one+small; return(scalb(x,exp)); }
209: else y=scalb(y,-exp);
210: return(scalb(sqrt(x*x+y*y),exp));
211: }
212:
213: else if(y==y) /* y is +-INF */
214: return(copysign(y,one));
215: else
216: return(y); /* y is NaN and x is finite */
217:
218: else if(x==x) /* x is +-INF */
219: return (copysign(x,one));
220: else if(finite(y))
221: return(x); /* x is NaN, y is finite */
222: else if(y!=y) return(y); /* x and y is NaN */
223: else return(copysign(y,one)); /* y is INF */
224: }
225: #endif
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.