125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2425c28e83SPiotr Jasiukajtis */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2725c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis */
2925c28e83SPiotr Jasiukajtis
3025c28e83SPiotr Jasiukajtis /* INDENT OFF */
3125c28e83SPiotr Jasiukajtis /*
3225c28e83SPiotr Jasiukajtis * void sincospi(double x, double *s, double *c)
3325c28e83SPiotr Jasiukajtis * *s = sin(pi*x); *c = cos(pi*x);
3425c28e83SPiotr Jasiukajtis *
3525c28e83SPiotr Jasiukajtis * Algorithm, 10/17/2002, K.C. Ng
3625c28e83SPiotr Jasiukajtis * ------------------------------
3725c28e83SPiotr Jasiukajtis * Let y = |4x|, z = floor(y), and n = (int)(z mod 8.0) (displayed in binary).
3825c28e83SPiotr Jasiukajtis * 1. If y == z, then x is a multiple of pi/4. Return the following values:
3925c28e83SPiotr Jasiukajtis * ---------------------------------------------------
4025c28e83SPiotr Jasiukajtis * n x mod 2 sin(x*pi) cos(x*pi) tan(x*pi)
4125c28e83SPiotr Jasiukajtis * ---------------------------------------------------
4225c28e83SPiotr Jasiukajtis * 000 0.00 +0 ___ +1 ___ +0
4325c28e83SPiotr Jasiukajtis * 001 0.25 +\/0.5 +\/0.5 +1
4425c28e83SPiotr Jasiukajtis * 010 0.50 +1 ___ +0 ___ +inf
4525c28e83SPiotr Jasiukajtis * 011 0.75 +\/0.5 -\/0.5 -1
4625c28e83SPiotr Jasiukajtis * 100 1.00 -0 ___ -1 ___ +0
4725c28e83SPiotr Jasiukajtis * 101 1.25 -\/0.5 -\/0.5 +1
4825c28e83SPiotr Jasiukajtis * 110 1.50 -1 ___ -0 ___ +inf
4925c28e83SPiotr Jasiukajtis * 111 1.75 -\/0.5 +\/0.5 -1
5025c28e83SPiotr Jasiukajtis * ---------------------------------------------------
5125c28e83SPiotr Jasiukajtis * 2. Otherwise,
5225c28e83SPiotr Jasiukajtis * ---------------------------------------------------
5325c28e83SPiotr Jasiukajtis * n t sin(x*pi) cos(x*pi) tan(x*pi)
5425c28e83SPiotr Jasiukajtis * ---------------------------------------------------
5525c28e83SPiotr Jasiukajtis * 000 (y-z)/4 sinpi(t) cospi(t) tanpi(t)
5625c28e83SPiotr Jasiukajtis * 001 (z+1-y)/4 cospi(t) sinpi(t) 1/tanpi(t)
5725c28e83SPiotr Jasiukajtis * 010 (y-z)/4 cospi(t) -sinpi(t) -1/tanpi(t)
5825c28e83SPiotr Jasiukajtis * 011 (z+1-y)/4 sinpi(t) -cospi(t) -tanpi(t)
5925c28e83SPiotr Jasiukajtis * 100 (y-z)/4 -sinpi(t) -cospi(t) tanpi(t)
6025c28e83SPiotr Jasiukajtis * 101 (z+1-y)/4 -cospi(t) -sinpi(t) 1/tanpi(t)
6125c28e83SPiotr Jasiukajtis * 110 (y-z)/4 -cospi(t) sinpi(t) -1/tanpi(t)
6225c28e83SPiotr Jasiukajtis * 111 (z+1-y)/4 -sinpi(t) cospi(t) -tanpi(t)
6325c28e83SPiotr Jasiukajtis * ---------------------------------------------------
6425c28e83SPiotr Jasiukajtis *
6525c28e83SPiotr Jasiukajtis * NOTE. This program compute sinpi/cospi(t<0.25) by __k_sin/cos(pi*t, 0.0).
6625c28e83SPiotr Jasiukajtis * This will return a result with error slightly more than one ulp (but less
6725c28e83SPiotr Jasiukajtis * than 2 ulp). If one wants accurate result, one may break up pi*t in
6825c28e83SPiotr Jasiukajtis * high (tpi_h) and low (tpi_l) parts and call __k_sin/cos(tip_h, tip_lo)
6925c28e83SPiotr Jasiukajtis * instead.
7025c28e83SPiotr Jasiukajtis */
7125c28e83SPiotr Jasiukajtis
7225c28e83SPiotr Jasiukajtis #include "libm.h"
7325c28e83SPiotr Jasiukajtis #include "libm_protos.h"
7425c28e83SPiotr Jasiukajtis #include "libm_macros.h"
7525c28e83SPiotr Jasiukajtis #include <math.h>
7625c28e83SPiotr Jasiukajtis #if defined(__SUNPRO_C)
7725c28e83SPiotr Jasiukajtis #include <sunmath.h>
7825c28e83SPiotr Jasiukajtis #endif
7925c28e83SPiotr Jasiukajtis
8025c28e83SPiotr Jasiukajtis static const double
8125c28e83SPiotr Jasiukajtis pi = 3.14159265358979323846, /* 400921FB,54442D18 */
8225c28e83SPiotr Jasiukajtis sqrth_h = 0.70710678118654757273731092936941422522068023681640625,
8325c28e83SPiotr Jasiukajtis sqrth_l = -4.8336466567264565185935844299127932213411660131004e-17;
8425c28e83SPiotr Jasiukajtis /* INDENT ON */
8525c28e83SPiotr Jasiukajtis
8625c28e83SPiotr Jasiukajtis void
sincospi(double x,double * s,double * c)87*709db5a3SToomas Soome sincospi(double x, double *s, double *c)
88*709db5a3SToomas Soome {
8925c28e83SPiotr Jasiukajtis double y, z, t;
9025c28e83SPiotr Jasiukajtis int n, ix, k;
91*709db5a3SToomas Soome int hx = ((int *)&x)[HIWORD];
92*709db5a3SToomas Soome unsigned h, lx = ((unsigned *)&x)[LOWORD];
9325c28e83SPiotr Jasiukajtis
9425c28e83SPiotr Jasiukajtis ix = hx & ~0x80000000;
9525c28e83SPiotr Jasiukajtis n = (ix >> 20) - 0x3ff;
9625c28e83SPiotr Jasiukajtis if (n >= 51) { /* |x| >= 2**51 */
97*709db5a3SToomas Soome if (n >= 1024) {
9825c28e83SPiotr Jasiukajtis #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
9925c28e83SPiotr Jasiukajtis *s = *c = ix >= 0x7ff80000 ? x : x - x;
10025c28e83SPiotr Jasiukajtis /* assumes sparc-like QNaN */
10125c28e83SPiotr Jasiukajtis #else
10225c28e83SPiotr Jasiukajtis *s = *c = x - x;
10325c28e83SPiotr Jasiukajtis #endif
104*709db5a3SToomas Soome } else {
105*709db5a3SToomas Soome if (n >= 53) {
10625c28e83SPiotr Jasiukajtis *s = 0.0;
10725c28e83SPiotr Jasiukajtis *c = 1.0;
108*709db5a3SToomas Soome } else if (n == 52) {
10925c28e83SPiotr Jasiukajtis if ((lx & 1) == 0) {
11025c28e83SPiotr Jasiukajtis *s = 0.0;
11125c28e83SPiotr Jasiukajtis *c = 1.0;
112*709db5a3SToomas Soome } else {
11325c28e83SPiotr Jasiukajtis *s = -0.0;
11425c28e83SPiotr Jasiukajtis *c = -1.0;
11525c28e83SPiotr Jasiukajtis }
116*709db5a3SToomas Soome } else { /* n == 51 */
11725c28e83SPiotr Jasiukajtis if ((lx & 1) == 0) {
11825c28e83SPiotr Jasiukajtis *s = 0.0;
11925c28e83SPiotr Jasiukajtis *c = 1.0;
120*709db5a3SToomas Soome } else {
12125c28e83SPiotr Jasiukajtis *s = 1.0;
12225c28e83SPiotr Jasiukajtis *c = 0.0;
12325c28e83SPiotr Jasiukajtis }
12425c28e83SPiotr Jasiukajtis if ((lx & 2) != 0) {
12525c28e83SPiotr Jasiukajtis *s = -*s;
12625c28e83SPiotr Jasiukajtis *c = -*c;
12725c28e83SPiotr Jasiukajtis }
12825c28e83SPiotr Jasiukajtis }
12925c28e83SPiotr Jasiukajtis }
130*709db5a3SToomas Soome } else if (n < -2) /* |x| < 0.25 */
13125c28e83SPiotr Jasiukajtis *s = __k_sincos(pi * fabs(x), 0.0, c);
13225c28e83SPiotr Jasiukajtis else {
13325c28e83SPiotr Jasiukajtis /* y = |4x|, z = floor(y), and n = (int)(z mod 8.0) */
13425c28e83SPiotr Jasiukajtis if (ix < 0x41C00000) { /* |x| < 2**29 */
13525c28e83SPiotr Jasiukajtis y = 4.0 * fabs(x);
136*709db5a3SToomas Soome n = (int)y; /* exact */
137*709db5a3SToomas Soome z = (double)n;
13825c28e83SPiotr Jasiukajtis k = z == y;
13925c28e83SPiotr Jasiukajtis t = (y - z) * 0.25;
140*709db5a3SToomas Soome } else { /* 2**29 <= |x| < 2**51 */
14125c28e83SPiotr Jasiukajtis y = fabs(x);
14225c28e83SPiotr Jasiukajtis k = 50 - n;
14325c28e83SPiotr Jasiukajtis n = lx >> k;
14425c28e83SPiotr Jasiukajtis h = n << k;
145*709db5a3SToomas Soome ((unsigned *)&z)[LOWORD] = h;
146*709db5a3SToomas Soome ((int *)&z)[HIWORD] = ix;
14725c28e83SPiotr Jasiukajtis k = h == lx;
14825c28e83SPiotr Jasiukajtis t = y - z;
14925c28e83SPiotr Jasiukajtis }
15025c28e83SPiotr Jasiukajtis if (k) { /* x = N/4 */
151*709db5a3SToomas Soome if ((n & 1) != 0) {
15225c28e83SPiotr Jasiukajtis *s = *c = sqrth_h + sqrth_l;
153*709db5a3SToomas Soome } else {
15425c28e83SPiotr Jasiukajtis if ((n & 2) == 0) {
15525c28e83SPiotr Jasiukajtis *s = 0.0;
15625c28e83SPiotr Jasiukajtis *c = 1.0;
157*709db5a3SToomas Soome } else {
15825c28e83SPiotr Jasiukajtis *s = 1.0;
15925c28e83SPiotr Jasiukajtis *c = 0.0;
16025c28e83SPiotr Jasiukajtis }
161*709db5a3SToomas Soome }
162*709db5a3SToomas Soome if ((n & 4) != 0)
163*709db5a3SToomas Soome *s = -*s;
164*709db5a3SToomas Soome if (((n + 1) & 4) != 0)
165*709db5a3SToomas Soome *c = -*c;
166*709db5a3SToomas Soome } else {
16725c28e83SPiotr Jasiukajtis if ((n & 1) != 0)
16825c28e83SPiotr Jasiukajtis t = 0.25 - t;
16925c28e83SPiotr Jasiukajtis if (((n + (n & 1)) & 2) == 0)
17025c28e83SPiotr Jasiukajtis *s = __k_sincos(pi * t, 0.0, c);
17125c28e83SPiotr Jasiukajtis else
17225c28e83SPiotr Jasiukajtis *c = __k_sincos(pi * t, 0.0, s);
173*709db5a3SToomas Soome if ((n & 4) != 0)
174*709db5a3SToomas Soome *s = -*s;
175*709db5a3SToomas Soome if (((n + 2) & 4) != 0)
176*709db5a3SToomas Soome *c = -*c;
17725c28e83SPiotr Jasiukajtis }
17825c28e83SPiotr Jasiukajtis }
17925c28e83SPiotr Jasiukajtis if (hx < 0)
18025c28e83SPiotr Jasiukajtis *s = -*s;
18125c28e83SPiotr Jasiukajtis }
182