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 rhypot(double x, double y)
49 *
50 * Method :
51 * 1. Special cases:
52 * x or y = Inf => 0
53 * x or y = NaN => QNaN
54 * x and y = 0 => Inf + divide-by-zero
55 * 2. Computes rhypot(x,y):
56 * rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
57 * Where:
58 * m = 1/max(|x|,|y|)
59 * xnm = x * m
60 * ynm = y * m
61 *
62 * Compute 1/(xnm * xnm + ynm * ynm) by simulating
63 * muti-precision arithmetic.
64 *
65 * Accuracy:
66 * Maximum error observed: less than 0.869 ulp after 1.000.000.000
67 * results.
68 */
69
70 #define sqrt __sqrt
71
72 extern double sqrt(double);
73
74 extern double fabs(double);
75
76 static const int __vlibm_TBL_rhypot[] = {
77 /* i = [0,127]
78 * TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
79 0x7fe00000, 0x7fdfc07f, 0x7fdf81f8, 0x7fdf4465,
80 0x7fdf07c1, 0x7fdecc07, 0x7fde9131, 0x7fde573a,
81 0x7fde1e1e, 0x7fdde5d6, 0x7fddae60, 0x7fdd77b6,
82 0x7fdd41d4, 0x7fdd0cb5, 0x7fdcd856, 0x7fdca4b3,
83 0x7fdc71c7, 0x7fdc3f8f, 0x7fdc0e07, 0x7fdbdd2b,
84 0x7fdbacf9, 0x7fdb7d6c, 0x7fdb4e81, 0x7fdb2036,
85 0x7fdaf286, 0x7fdac570, 0x7fda98ef, 0x7fda6d01,
86 0x7fda41a4, 0x7fda16d3, 0x7fd9ec8e, 0x7fd9c2d1,
87 0x7fd99999, 0x7fd970e4, 0x7fd948b0, 0x7fd920fb,
88 0x7fd8f9c1, 0x7fd8d301, 0x7fd8acb9, 0x7fd886e5,
89 0x7fd86186, 0x7fd83c97, 0x7fd81818, 0x7fd7f405,
90 0x7fd7d05f, 0x7fd7ad22, 0x7fd78a4c, 0x7fd767dc,
91 0x7fd745d1, 0x7fd72428, 0x7fd702e0, 0x7fd6e1f7,
92 0x7fd6c16c, 0x7fd6a13c, 0x7fd68168, 0x7fd661ec,
93 0x7fd642c8, 0x7fd623fa, 0x7fd60581, 0x7fd5e75b,
411 itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
412 itbl1 -= iexp1;
413 HI(&dd1) = itbl1;
414 LO(&dd1) = 0;
415
416 dd1 = dd1 * (DTWO - dd1 * dres1);
417 dd1 = dd1 * (DTWO - dd1 * dres1);
418 dres1 = dd1 * (DTWO - dd1 * dres1);
419
420 HI(&res1) = HI(&dres1) & 0xffffff00;
421 LO(&res1) = 0;
422 res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
423 res1 = sqrt (res1);
424
425 res1 = scl1 * res1;
426
427 *pz1 = res1;
428 }
429 }
430 }
431
|
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 rhypot(double x, double y)
48 *
49 * Method :
50 * 1. Special cases:
51 * x or y = Inf => 0
52 * x or y = NaN => QNaN
53 * x and y = 0 => Inf + divide-by-zero
54 * 2. Computes rhypot(x,y):
55 * rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
56 * Where:
57 * m = 1/max(|x|,|y|)
58 * xnm = x * m
59 * ynm = y * m
60 *
61 * Compute 1/(xnm * xnm + ynm * ynm) by simulating
62 * muti-precision arithmetic.
63 *
64 * Accuracy:
65 * Maximum error observed: less than 0.869 ulp after 1.000.000.000
66 * results.
67 */
68
69 extern double sqrt(double);
70 extern double fabs(double);
71
72 static const int __vlibm_TBL_rhypot[] = {
73 /* i = [0,127]
74 * TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
75 0x7fe00000, 0x7fdfc07f, 0x7fdf81f8, 0x7fdf4465,
76 0x7fdf07c1, 0x7fdecc07, 0x7fde9131, 0x7fde573a,
77 0x7fde1e1e, 0x7fdde5d6, 0x7fddae60, 0x7fdd77b6,
78 0x7fdd41d4, 0x7fdd0cb5, 0x7fdcd856, 0x7fdca4b3,
79 0x7fdc71c7, 0x7fdc3f8f, 0x7fdc0e07, 0x7fdbdd2b,
80 0x7fdbacf9, 0x7fdb7d6c, 0x7fdb4e81, 0x7fdb2036,
81 0x7fdaf286, 0x7fdac570, 0x7fda98ef, 0x7fda6d01,
82 0x7fda41a4, 0x7fda16d3, 0x7fd9ec8e, 0x7fd9c2d1,
83 0x7fd99999, 0x7fd970e4, 0x7fd948b0, 0x7fd920fb,
84 0x7fd8f9c1, 0x7fd8d301, 0x7fd8acb9, 0x7fd886e5,
85 0x7fd86186, 0x7fd83c97, 0x7fd81818, 0x7fd7f405,
86 0x7fd7d05f, 0x7fd7ad22, 0x7fd78a4c, 0x7fd767dc,
87 0x7fd745d1, 0x7fd72428, 0x7fd702e0, 0x7fd6e1f7,
88 0x7fd6c16c, 0x7fd6a13c, 0x7fd68168, 0x7fd661ec,
89 0x7fd642c8, 0x7fd623fa, 0x7fd60581, 0x7fd5e75b,
407 itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
408 itbl1 -= iexp1;
409 HI(&dd1) = itbl1;
410 LO(&dd1) = 0;
411
412 dd1 = dd1 * (DTWO - dd1 * dres1);
413 dd1 = dd1 * (DTWO - dd1 * dres1);
414 dres1 = dd1 * (DTWO - dd1 * dres1);
415
416 HI(&res1) = HI(&dres1) & 0xffffff00;
417 LO(&res1) = 0;
418 res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
419 res1 = sqrt (res1);
420
421 res1 = scl1 * res1;
422
423 *pz1 = res1;
424 }
425 }
426 }
|