Commit | Line | Data |
---|---|---|
98cc7428 | 1 | /* |
f9fea09f | 2 | * Copyright (c) 1985 Regents of the University of California. |
98cc7428 KB |
3 | * All rights reserved. |
4 | * | |
af359dea C |
5 | * Redistribution and use in source and binary forms, with or without |
6 | * modification, are permitted provided that the following conditions | |
7 | * are met: | |
8 | * 1. Redistributions of source code must retain the above copyright | |
9 | * notice, this list of conditions and the following disclaimer. | |
10 | * 2. Redistributions in binary form must reproduce the above copyright | |
11 | * notice, this list of conditions and the following disclaimer in the | |
12 | * documentation and/or other materials provided with the distribution. | |
13 | * 3. All advertising materials mentioning features or use of this software | |
14 | * must display the following acknowledgement: | |
15 | * This product includes software developed by the University of | |
16 | * California, Berkeley and its contributors. | |
17 | * 4. Neither the name of the University nor the names of its contributors | |
18 | * may be used to endorse or promote products derived from this software | |
19 | * without specific prior written permission. | |
98cc7428 | 20 | * |
af359dea C |
21 | * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND |
22 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
23 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
24 | * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE | |
25 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
26 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
27 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
28 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
29 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
30 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
31 | * SUCH DAMAGE. | |
f9fea09f ZAL |
32 | */ |
33 | ||
34 | #ifndef lint | |
af359dea | 35 | static char sccsid[] = "@(#)cabs.c 5.6 (Berkeley) 10/9/90"; |
98cc7428 | 36 | #endif /* not lint */ |
f9fea09f | 37 | |
f9fea09f ZAL |
38 | /* HYPOT(X,Y) |
39 | * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY | |
40 | * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) | |
41 | * CODED IN C BY K.C. NG, 11/28/84; | |
42 | * REVISED BY K.C. NG, 7/12/85. | |
43 | * | |
44 | * Required system supported functions : | |
45 | * copysign(x,y) | |
46 | * finite(x) | |
47 | * scalb(x,N) | |
48 | * sqrt(x) | |
49 | * | |
50 | * Method : | |
51 | * 1. replace x by |x| and y by |y|, and swap x and | |
52 | * y if y > x (hence x is never smaller than y). | |
53 | * 2. Hypot(x,y) is computed by: | |
54 | * Case I, x/y > 2 | |
55 | * | |
56 | * y | |
57 | * hypot = x + ----------------------------- | |
58 | * 2 | |
59 | * sqrt ( 1 + [x/y] ) + x/y | |
60 | * | |
61 | * Case II, x/y <= 2 | |
62 | * y | |
63 | * hypot = x + -------------------------------------------------- | |
64 | * 2 | |
65 | * [x/y] - 2 | |
66 | * (sqrt(2)+1) + (x-y)/y + ----------------------------- | |
67 | * 2 | |
68 | * sqrt ( 1 + [x/y] ) + sqrt(2) | |
69 | * | |
70 | * | |
71 | * | |
72 | * Special cases: | |
73 | * hypot(x,y) is INF if x or y is +INF or -INF; else | |
74 | * hypot(x,y) is NAN if x or y is NAN. | |
75 | * | |
76 | * Accuracy: | |
77 | * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units | |
78 | * in the last place). See Kahan's "Interval Arithmetic Options in the | |
79 | * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics | |
80 | * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate | |
81 | * code follows in comments.) In a test run with 500,000 random arguments | |
82 | * on a VAX, the maximum observed error was .959 ulps. | |
83 | * | |
84 | * Constants: | |
85 | * The hexadecimal values are the intended ones for the following constants. | |
86 | * The decimal values may be used, provided that the compiler will convert | |
87 | * from decimal to binary accurately enough to produce the hexadecimal values | |
88 | * shown. | |
89 | */ | |
9eda3584 | 90 | #include "mathimpl.h" |
f9fea09f | 91 | |
9eda3584 KB |
92 | vc(r2p1hi, 2.4142135623730950345E0 ,8279,411a,ef32,99fc, 2, .9A827999FCEF32) |
93 | vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B) | |
94 | vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) | |
95 | ||
96 | ic(r2p1hi, 2.4142135623730949234E0 , 1, 1.3504F333F9DE6) | |
97 | ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5) | |
98 | ic(sqrt2, 1.4142135623730951455E0 , 0, 1.6A09E667F3BCD) | |
99 | ||
100 | #ifdef vccast | |
101 | #define r2p1hi vccast(r2p1hi) | |
102 | #define r2p1lo vccast(r2p1lo) | |
103 | #define sqrt2 vccast(sqrt2) | |
104 | #endif | |
f9fea09f | 105 | |
74920388 ZAL |
106 | double |
107 | hypot(x,y) | |
f9fea09f ZAL |
108 | double x, y; |
109 | { | |
9eda3584 | 110 | static const double zero=0, one=1, |
f9fea09f | 111 | small=1.0E-18; /* fl(1+small)==1 */ |
9eda3584 KB |
112 | static const ibig=30; /* fl(1+2**(2*ibig))==1 */ |
113 | double t,r; | |
114 | int exp; | |
f9fea09f ZAL |
115 | |
116 | if(finite(x)) | |
117 | if(finite(y)) | |
118 | { | |
119 | x=copysign(x,one); | |
120 | y=copysign(y,one); | |
121 | if(y > x) | |
122 | { t=x; x=y; y=t; } | |
123 | if(x == zero) return(zero); | |
124 | if(y == zero) return(x); | |
125 | exp= logb(x); | |
126 | if(exp-(int)logb(y) > ibig ) | |
127 | /* raise inexact flag and return |x| */ | |
128 | { one+small; return(x); } | |
129 | ||
130 | /* start computing sqrt(x^2 + y^2) */ | |
131 | r=x-y; | |
132 | if(r>y) { /* x/y > 2 */ | |
133 | r=x/y; | |
134 | r=r+sqrt(one+r*r); } | |
135 | else { /* 1 <= x/y <= 2 */ | |
136 | r/=y; t=r*(r+2.0); | |
137 | r+=t/(sqrt2+sqrt(2.0+t)); | |
138 | r+=r2p1lo; r+=r2p1hi; } | |
139 | ||
140 | r=y/r; | |
141 | return(x+r); | |
142 | ||
143 | } | |
144 | ||
145 | else if(y==y) /* y is +-INF */ | |
146 | return(copysign(y,one)); | |
147 | else | |
148 | return(y); /* y is NaN and x is finite */ | |
149 | ||
150 | else if(x==x) /* x is +-INF */ | |
151 | return (copysign(x,one)); | |
152 | else if(finite(y)) | |
153 | return(x); /* x is NaN, y is finite */ | |
859dc438 | 154 | #if !defined(vax)&&!defined(tahoe) |
f9fea09f | 155 | else if(y!=y) return(y); /* x and y is NaN */ |
859dc438 | 156 | #endif /* !defined(vax)&&!defined(tahoe) */ |
f9fea09f ZAL |
157 | else return(copysign(y,one)); /* y is INF */ |
158 | } | |
159 | ||
74920388 ZAL |
160 | /* CABS(Z) |
161 | * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY | |
162 | * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) | |
163 | * CODED IN C BY K.C. NG, 11/28/84. | |
164 | * REVISED BY K.C. NG, 7/12/85. | |
165 | * | |
166 | * Required kernel function : | |
167 | * hypot(x,y) | |
168 | * | |
169 | * Method : | |
170 | * cabs(z) = hypot(x,y) . | |
171 | */ | |
172 | ||
173 | double | |
174 | cabs(z) | |
175 | struct { double x, y;} z; | |
176 | { | |
177 | return hypot(z.x,z.y); | |
178 | } | |
179 | ||
180 | double | |
181 | z_abs(z) | |
182 | struct { double x,y;} *z; | |
183 | { | |
184 | return hypot(z->x,z->y); | |
185 | } | |
186 | ||
f9fea09f ZAL |
187 | /* A faster but less accurate version of cabs(x,y) */ |
188 | #if 0 | |
189 | double hypot(x,y) | |
190 | double x, y; | |
191 | { | |
9eda3584 | 192 | static const double zero=0, one=1; |
f9fea09f | 193 | small=1.0E-18; /* fl(1+small)==1 */ |
9eda3584 KB |
194 | static const ibig=30; /* fl(1+2**(2*ibig))==1 */ |
195 | double temp; | |
196 | int exp; | |
f9fea09f ZAL |
197 | |
198 | if(finite(x)) | |
199 | if(finite(y)) | |
200 | { | |
201 | x=copysign(x,one); | |
202 | y=copysign(y,one); | |
203 | if(y > x) | |
204 | { temp=x; x=y; y=temp; } | |
205 | if(x == zero) return(zero); | |
206 | if(y == zero) return(x); | |
207 | exp= logb(x); | |
208 | x=scalb(x,-exp); | |
209 | if(exp-(int)logb(y) > ibig ) | |
210 | /* raise inexact flag and return |x| */ | |
211 | { one+small; return(scalb(x,exp)); } | |
212 | else y=scalb(y,-exp); | |
213 | return(scalb(sqrt(x*x+y*y),exp)); | |
214 | } | |
215 | ||
216 | else if(y==y) /* y is +-INF */ | |
217 | return(copysign(y,one)); | |
218 | else | |
219 | return(y); /* y is NaN and x is finite */ | |
220 | ||
221 | else if(x==x) /* x is +-INF */ | |
222 | return (copysign(x,one)); | |
223 | else if(finite(y)) | |
224 | return(x); /* x is NaN, y is finite */ | |
225 | else if(y!=y) return(y); /* x and y is NaN */ | |
226 | else return(copysign(y,one)); /* y is INF */ | |
227 | } | |
228 | #endif |