1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #include <sys/isa_defs.h>
31 #include "libm_inlines.h"
32
33 #ifdef _LITTLE_ENDIAN
34 #define HI(x) *(1+(int*)x)
35 #define LO(x) *(unsigned*)x
36 #else
37 #define HI(x) *(int*)x
38 #define LO(x) *(1+(unsigned*)x)
39 #endif
40
41 #ifdef __RESTRICT
42 #define restrict _Restrict
43 #else
44 #define restrict
45 #endif
46
47 /* double hypot(double x, double y)
48 *
49 * Method :
50 * 1. Special cases:
51 * x or y is +Inf or -Inf => +Inf
52 * x or y is NaN => QNaN
53 * 2. Computes hypot(x,y):
54 * hypot(x,y) = m * sqrt(xnm * xnm + ynm * ynm)
55 * Where:
56 * m = max(|x|,|y|)
57 * xnm = x * (1/m)
58 * ynm = y * (1/m)
59 *
60 * Compute xnm * xnm + ynm * ynm by simulating
61 * muti-precision arithmetic.
62 *
63 * Accuracy:
64 * Maximum error observed: less than 0.872 ulp after 16.777.216.000
65 * results.
66 */
67
68 extern double sqrt(double);
69 extern double fabs(double);
70
71 static const unsigned long long LCONST[] = {
72 0x41b0000000000000ULL, /* D2ON28 = 2 ** 28 */
73 0x0010000000000000ULL, /* D2ONM1022 = 2 ** -1022 */
74 0x7fd0000000000000ULL /* D2ONP1022 = 2 ** 1022 */
75 };
76
77 static void
78 __vhypot_n(int n, double * restrict px, int stridex, double * restrict py,
79 int stridey, double * restrict pz, int stridez);
80
81 #pragma no_inline(__vhypot_n)
82
83 #define RETURN(ret) \
84 { \
85 *pz = (ret); \
86 py += stridey; \
87 pz += stridez; \
88 if (n_n == 0) \
89 { \
90 hx0 = HI(px); \
91 hy0 = HI(py); \
92 spx = px; spy = py; spz = pz; \
93 continue; \
94 } \
95 n--; \
96 break; \
97 }
98
99 void
100 __vhypot(int n, double * restrict px, int stridex, double * restrict py,
101 int stridey, double * restrict pz, int stridez)
102 {
103 int hx0, hx1, hy0, j0, diff;
104 double x_hi, x_lo, y_hi, y_lo;
105 double scl = 0;
106 double x, y, res;
107 double *spx, *spy, *spz;
108 int n_n;
109 double D2ON28 = ((double*)LCONST)[0]; /* 2 ** 28 */
110 double D2ONM1022 = ((double*)LCONST)[1]; /* 2 **-1022 */
111 double D2ONP1022 = ((double*)LCONST)[2]; /* 2 ** 1022 */
112
113 while (n > 1)
114 {
115 n_n = 0;
116 spx = px;
117 spy = py;
118 spz = pz;
119 hx0 = HI(px);
120 hy0 = HI(py);
121 for (; n > 1 ; n--)
122 {
123 px += stridex;
124 hx0 &= 0x7fffffff;
125 hy0 &= 0x7fffffff;
126
127 if (hx0 >= 0x7fe00000) /* |X| >= 2**1023 or Inf or NaN */
128 {
129 diff = hy0 - hx0;
130 j0 = diff >> 31;
131 j0 = hy0 - (diff & j0);
132 j0 &= 0x7ff00000;
133 x = *(px - stridex);
134 y = *py;
135 x = fabs(x);
136 y = fabs(y);
137 if (j0 >= 0x7ff00000) /* |X| or |Y| = Inf or NaN */
138 {
139 int lx = LO((px - stridex));
140 int ly = LO(py);
141 if (hx0 == 0x7ff00000 && lx == 0) res = x == y ? y : x;
142 else if (hy0 == 0x7ff00000 && ly == 0) res = x == y ? x : y;
143 else res = x + y;
144 RETURN (res)
145 }
146 else
147 {
148 j0 = diff >> 31;
149 if (((diff ^ j0) - j0) < 0x03600000) /* max(|X|,|Y|)/min(|X|,|Y|) < 2**54 */
150 {
151 x *= D2ONM1022;
152 y *= D2ONM1022;
153
154 x_hi = (x + D2ON28) - D2ON28;
155 x_lo = x - x_hi;
156 y_hi = (y + D2ON28) - D2ON28;
157 y_lo = y - y_hi;
158 res = (x_hi * x_hi + y_hi * y_hi);
159 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
160
161 res = sqrt (res);
162
163 res = D2ONP1022 * res;
164 RETURN (res)
165 }
166 else RETURN (x + y)
167 }
168 }
169 if (hy0 >= 0x7fe00000) /* |Y| >= 2**1023 or Inf or NaN */
170 {
171 diff = hy0 - hx0;
172 j0 = diff >> 31;
173 j0 = hy0 - (diff & j0);
174 j0 &= 0x7ff00000;
175 x = *(px - stridex);
176 y = *py;
177 x = fabs(x);
178 y = fabs(y);
179 if (j0 >= 0x7ff00000) /* |X| or |Y| = Inf or NaN */
180 {
181 int lx = LO((px - stridex));
182 int ly = LO(py);
183 if (hx0 == 0x7ff00000 && lx == 0) res = x == y ? y : x;
184 else if (hy0 == 0x7ff00000 && ly == 0) res = x == y ? x : y;
185 else res = x + y;
186 RETURN (res)
187 }
188 else
189 {
190 j0 = diff >> 31;
191 if (((diff ^ j0) - j0) < 0x03600000) /* max(|X|,|Y|)/min(|X|,|Y|) < 2**54 */
192 {
193 x *= D2ONM1022;
194 y *= D2ONM1022;
195
196 x_hi = (x + D2ON28) - D2ON28;
197 x_lo = x - x_hi;
198 y_hi = (y + D2ON28) - D2ON28;
199 y_lo = y - y_hi;
200 res = (x_hi * x_hi + y_hi * y_hi);
201 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
202
203 res = sqrt (res);
204
205 res = D2ONP1022 * res;
206 RETURN (res)
207 }
208 else RETURN (x + y)
209 }
210 }
211
212 hx1 = HI(px);
213
214 if (hx0 < 0x00100000 && hy0 < 0x00100000) /* X and Y are subnormal */
215 {
216 x = *(px - stridex);
217 y = *py;
218
219 x *= D2ONP1022;
220 y *= D2ONP1022;
221
222 x_hi = (x + D2ON28) - D2ON28;
223 x_lo = x - x_hi;
224 y_hi = (y + D2ON28) - D2ON28;
225 y_lo = y - y_hi;
226 res = (x_hi * x_hi + y_hi * y_hi);
227 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
228
229 res = sqrt(res);
230
231 res = D2ONM1022 * res;
232 RETURN (res)
233 }
234
235 hx0 = hx1;
236 py += stridey;
237 pz += stridez;
238 n_n++;
239 hy0 = HI(py);
240 }
241 if (n_n > 0)
242 __vhypot_n (n_n, spx, stridex, spy, stridey, spz, stridez);
243 }
244
245 if (n > 0)
246 {
247 x = *px;
248 y = *py;
249 hx0 = HI(px);
250 hy0 = HI(py);
251
252 hx0 &= 0x7fffffff;
253 hy0 &= 0x7fffffff;
254
255 diff = hy0 - hx0;
256 j0 = diff >> 31;
257 j0 = hy0 - (diff & j0);
258 j0 &= 0x7ff00000;
259
260 if (j0 >= 0x7fe00000) /* max(|X|,|Y|) >= 2**1023 or X or Y = Inf or NaN */
261 {
262 x = fabs(x);
263 y = fabs(y);
264 if (j0 >= 0x7ff00000) /* |X| or |Y| = Inf or NaN */
265 {
266 int lx = LO(px);
267 int ly = LO(py);
268 if (hx0 == 0x7ff00000 && lx == 0) res = x == y ? y : x;
269 else if (hy0 == 0x7ff00000 && ly == 0) res = x == y ? x : y;
270 else res = x + y;
271 *pz = res;
272 return;
273 }
274 else
275 {
276 j0 = diff >> 31;
277 if (((diff ^ j0) - j0) < 0x03600000) /* max(|X|,|Y|)/min(|X|,|Y|) < 2**54 */
278 {
279 x *= D2ONM1022;
280 y *= D2ONM1022;
281
282 x_hi = (x + D2ON28) - D2ON28;
283 x_lo = x - x_hi;
284 y_hi = (y + D2ON28) - D2ON28;
285 y_lo = y - y_hi;
286 res = (x_hi * x_hi + y_hi * y_hi);
287 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
288
289 res = sqrt (res);
290
291 res = D2ONP1022 * res;
292 *pz = res;
293 return;
294 }
295 else
296 {
297 *pz = x + y;
298 return;
299 }
300 }
301 }
302
303 if (j0 < 0x00100000) /* X and Y are subnormal */
304 {
305 x *= D2ONP1022;
306 y *= D2ONP1022;
307
308 x_hi = (x + D2ON28) - D2ON28;
309 x_lo = x - x_hi;
310 y_hi = (y + D2ON28) - D2ON28;
311 y_lo = y - y_hi;
312 res = (x_hi * x_hi + y_hi * y_hi);
313 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
314
315 res = sqrt(res);
316
317 res = D2ONM1022 * res;
318 *pz = res;
319 return;
320 }
321
322 HI(&scl) = (0x7fe00000 - j0);
323
324 x *= scl;
325 y *= scl;
326
327 x_hi = (x + D2ON28) - D2ON28;
328 y_hi = (y + D2ON28) - D2ON28;
329 x_lo = x - x_hi;
330 y_lo = y - y_hi;
331
332 res = (x_hi * x_hi + y_hi * y_hi);
333 res += ((x + x_hi) * x_lo + (y + y_hi) * y_lo);
334
335 res = sqrt(res);
336
337 HI(&scl) = j0;
338
339 res = scl * res;
340 *pz = res;
341 }
342 }
343
344 static void
345 __vhypot_n(int n, double * restrict px, int stridex, double * restrict py,
346 int stridey, double * restrict pz, int stridez)
347 {
348 int hx0, hy0, j0, diff0;
349 double x_hi0, x_lo0, y_hi0, y_lo0, scl0 = 0;
350 double x0, y0, res0;
351 double D2ON28 = ((double*)LCONST)[0]; /* 2 ** 28 */
352
353 for(; n > 0 ; n--)
354 {
355 x0 = *px;
356 y0 = *py;
357 hx0 = HI(px);
358 hy0 = HI(py);
359
360 hx0 &= 0x7fffffff;
361 hy0 &= 0x7fffffff;
362
363 diff0 = hy0 - hx0;
364 j0 = diff0 >> 31;
365 j0 = hy0 - (diff0 & j0);
366 j0 &= 0x7ff00000;
367
368 px += stridex;
369 py += stridey;
370
371 HI(&scl0) = (0x7fe00000 - j0);
372
373 x0 *= scl0;
374 y0 *= scl0;
375
376 x_hi0 = (x0 + D2ON28) - D2ON28;
377 y_hi0 = (y0 + D2ON28) - D2ON28;
378 x_lo0 = x0 - x_hi0;
379 y_lo0 = y0 - y_hi0;
380
381 res0 = (x_hi0 * x_hi0 + y_hi0 * y_hi0);
382 res0 += ((x0 + x_hi0) * x_lo0 + (y0 + y_hi0) * y_lo0);
383
384 res0 = sqrt(res0);
385
386 HI(&scl0) = j0;
387
388 res0 = scl0 * res0;
389 *pz = res0;
390
391 pz += stridez;
392 }
393 }