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