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
32 #ifdef _LITTLE_ENDIAN
33 #define HI(x) *(1+(int*)x)
34 #define LO(x) *(unsigned*)x
35 #else
36 #define HI(x) *(int*)x
37 #define LO(x) *(1+(unsigned*)x)
38 #endif
39
40 #ifdef __RESTRICT
41 #define restrict _Restrict
42 #else
43 #define restrict
44 #endif
45
46 extern const double __vlibm_TBL_sincos_hi[], __vlibm_TBL_sincos_lo[];
47 extern int __vlibm_rem_pio2m( double *, double *, int, int, int );
48
49 static const double
50 zero = 0.0,
51 one = 1.0,
52 two24 = 16777216.0,
53 pp1 = -1.666666666605760465276263943134982554676e-0001,
54 pp2 = 8.333261209690963126718376566146180944442e-0003,
55 p1 = -1.666666666666629669805215138920301589656e-0001,
56 p2 = 8.333333332390951295683993455280336376663e-0003,
57 p3 = -1.984126237997976692791551778230098403960e-0004,
58 p4 = 2.753403624854277237649987622848330351110e-0006,
59 qq1 = -4.999999999977710986407023955908711557870e-0001,
60 qq2 = 4.166654863857219350645055881018842089580e-0002,
61 q1 = -4.999999999999931701464060878888294524481e-0001,
62 q2 = 4.166666666394861917535640593963708222319e-0002,
63 q3 = -1.388888552656142867832756687736851681462e-0003,
64 q4 = 2.478519423681460796618128289454530524759e-0005;
65
66 void
67 __vlibm_vsincos_big( int n, double * restrict x, int stridex,
68 double * restrict ss, int stridess,
69 double * restrict cc, int stridecc, int thresh )
70 {
71 for ( ; n--; x += stridex, ss += stridess, cc += stridecc )
72 {
73 double ts, tc, tx, tt[3], ty[2], t, w, z, c, s;
74 unsigned hx, xsb;
75 int e0, nx, j;
76
77 hx = HI(x);
78 xsb = hx & 0x80000000;
79 hx &= ~0x80000000;
80 if ( hx <= thresh || hx >= 0x7ff00000 )
81 continue;
82
83 /*
84 * Argument reduction part.
85 */
86 e0 = ( hx >> 20 ) - 1046;
87 HI(&tx) = 0x41600000 | ( hx & 0xfffff );
88 LO(&tx) = LO(x);
89 tt[0] = (double)( (int) tx );
90 tx = ( tx - tt[0] ) * two24;
91 if ( tx != zero )
92 {
93 nx = 2;
94 tt[1] = (double)( (int) tx );
95 tt[2] = ( tx - tt[1] ) * two24;
96 if ( tt[2] != zero )
97 nx = 3;
98 }
99 else
100 {
101 nx = 1;
102 tt[1] = tt[2] = zero;
103 }
104 nx = __vlibm_rem_pio2m( tt, ty, e0, nx, 2 );
105 if ( xsb )
106 {
107 nx = -nx;
108 ty[0] = -ty[0];
109 ty[1] = -ty[1];
110 }
111
112 /* now nx and ty[*] are the quadrant and reduced arg */
113 hx = HI(&ty[0]);
114 xsb = 0;
115 if ( hx & 0x80000000 )
116 {
117 ty[0] = -ty[0];
118 ty[1] = -ty[1];
119 hx &= ~0x80000000;
120 xsb = 1;
121 }
122 if ( hx < 0x3fc40000 )
123 {
124 z = ty[0] * ty[0];
125 t = z * ( q1 + z * ( q2 + z * ( q3 + z * q4 ) ) );
126 c = one + t;
127 t = z * ( p1 + z * ( p2 + z * ( p3 + z * p4 ) ) );
128 s = ty[0] + ( ty[1] + ty[0] * t );
129 }
130 else {
131 j = ( hx + 0x4000 ) & 0x7fff8000;
132 HI(&t) = j;
133 LO(&t) = 0;
134 ty[0] = ( ty[0] - t ) + ty[1];
135 z = ty[0] * ty[0];
136 t = z * ( qq1 + z * qq2 );
137 w = ty[0] * ( one + z * ( pp1 + z * pp2 ) );
138 j = ( ( j - 0x3fc40000 ) >> 13 ) & ~3;
139
140 c = __vlibm_TBL_sincos_hi[j+1];
141 tc = __vlibm_TBL_sincos_lo[j+1] - ( __vlibm_TBL_sincos_hi[j] * w - c * t );
142 c += tc;
143
144 s = __vlibm_TBL_sincos_hi[j];
145 ts = ( __vlibm_TBL_sincos_hi[j+1] * w + s * t ) + __vlibm_TBL_sincos_lo[j];
146 s += ts;
147 }
148 if ( xsb ) {
149 s = -s;
150 }
151
152 switch ( nx & 3 ) {
153 case 0:
154 *ss = s;
155 *cc = c;
156 break;
157
158 case 1:
159 *ss = c;
160 *cc = -s;
161 break;
162
163 case 2:
164 *ss = -s;
165 *cc = -c;
166 break;
167
168 case 3:
169 *ss = -c;
170 *cc = s;
171 break;
172 }
173 }
174 }