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_synonyms.h"
32 #include "libm_inlines.h"
33
34 #ifdef _LITTLE_ENDIAN
35 #define HI(x) *(1+(int*)x)
36 #define LO(x) *(unsigned*)x
37 #else
38 #define HI(x) *(int*)x
39 #define LO(x) *(1+(unsigned*)x)
40 #endif
41
42 #ifdef __RESTRICT
43 #define restrict _Restrict
44 #else
45 #define restrict
46 #endif
47
48 /* double hypot(double x, double y)
49 *
50 * Method :
51 * 1. Special cases:
52 * x or y is +Inf or -Inf => +Inf
53 * x or y is NaN => QNaN
54 * 2. Computes hypot(x,y):
55 * hypot(x,y) = m * sqrt(xnm * xnm + ynm * ynm)
56 * Where:
57 * m = max(|x|,|y|)
58 * xnm = x * (1/m)
59 * ynm = y * (1/m)
60 *
61 * Compute xnm * xnm + ynm * ynm by simulating
62 * muti-precision arithmetic.
63 *
64 * Accuracy:
65 * Maximum error observed: less than 0.872 ulp after 16.777.216.000
66 * results.
67 */
68
69 #define sqrt __sqrt
70
71 extern double sqrt(double);
72 extern double fabs(double);
73
74 static const unsigned long long LCONST[] = {
75 0x41b0000000000000ULL, /* D2ON28 = 2 ** 28 */
76 0x0010000000000000ULL, /* D2ONM1022 = 2 ** -1022 */
77 0x7fd0000000000000ULL /* D2ONP1022 = 2 ** 1022 */
78 };
79
80 static void
81 __vhypot_n(int n, double * restrict px, int stridex, double * restrict py,
82 int stridey, double * restrict pz, int stridez);
83
84 #pragma no_inline(__vhypot_n)
85
86 #define RETURN(ret) \
87 { \
88 *pz = (ret); \
89 py += stridey; \
90 pz += stridez; \
377 y0 *= scl0;
378
379 x_hi0 = (x0 + D2ON28) - D2ON28;
380 y_hi0 = (y0 + D2ON28) - D2ON28;
381 x_lo0 = x0 - x_hi0;
382 y_lo0 = y0 - y_hi0;
383
384 res0 = (x_hi0 * x_hi0 + y_hi0 * y_hi0);
385 res0 += ((x0 + x_hi0) * x_lo0 + (y0 + y_hi0) * y_lo0);
386
387 res0 = sqrt(res0);
388
389 HI(&scl0) = j0;
390
391 res0 = scl0 * res0;
392 *pz = res0;
393
394 pz += stridez;
395 }
396 }
397
|
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; \
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 }
|