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 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __sincosf = sincosf
  32 
  33 
  34 /*
  35  * For |x| < pi/4, let z = x * x, and approximate sin(x) by
  36  *
  37  *      S(x) = x(S0 + S1*z)(S2 + S3*z + z*z)
  38  * where
  39  *      S0 =   1.85735322054308378716204874632872525989806770558e-0003,
  40  *      S1 =  -1.95035094218403635082921458859320791358115801259e-0004,
  41  *      S2 =   5.38400550766074785970952495168558701485841707252e+0002,
  42  *      S3 =  -3.31975110777873728964197739157371509422022905947e+0001,
  43  *
  44  * with error bounded by |(sin(x) - S(x))/x| < 2**(-28.2), and
  45  * cos(x) by
  46  *
  47  *      C(x) = (C0 + C1*z + C2*z*z) * (C3 + C4*z + z*z)
  48  * where
  49  *      C0 =   1.09349482127188401868272000389539985058873853699e-0003
  50  *      C1 =  -5.03324285989964979398034700054920226866107675091e-0004
  51  *      C2 =   2.43792880266971107750418061559602239831538067410e-0005
  52  *      C3 =   9.14499072605666582228127405245558035523741471271e+0002
  53  *      C4 =  -3.63151270591815439197122504991683846785293207730e+0001
  54  *
  55  * with error bounded by |cos(x) - C(x)| < 2**(-34.2).
  56  */
  57 
  58 #include "libm.h"
  59 
  60 extern const int _TBL_ipio2_inf[];
  61 extern int __rem_pio2m(double *, double *, int, int, int, const int *);
  62 
  63 #if defined(__i386) && !defined(__amd64)
  64 extern int __swapRP(int);
  65 #endif
  66 
  67 static const double C[] = {
  68         1.85735322054308378716204874632872525989806770558e-0003,
  69         -1.95035094218403635082921458859320791358115801259e-0004,
  70         5.38400550766074785970952495168558701485841707252e+0002,
  71         -3.31975110777873728964197739157371509422022905947e+0001,
  72         1.09349482127188401868272000389539985058873853699e-0003,
  73         -5.03324285989964979398034700054920226866107675091e-0004,
  74         2.43792880266971107750418061559602239831538067410e-0005,
  75         9.14499072605666582228127405245558035523741471271e+0002,
  76         -3.63151270591815439197122504991683846785293207730e+0001,
  77         0.636619772367581343075535,     /* 2^ -1  * 1.45F306DC9C883 */
  78         0.5,
  79         1.570796326734125614166,        /* 2^  0  * 1.921FB54400000 */
  80         6.077100506506192601475e-11,    /* 2^-34  * 1.0B4611A626331 */
  81 };
  82 
  83 #define S0              C[0]
  84 #define S1              C[1]
  85 #define S2              C[2]
  86 #define S3              C[3]
  87 #define C0              C[4]
  88 #define C1              C[5]
  89 #define C2              C[6]
  90 #define C3              C[7]
  91 #define C4              C[8]
  92 #define invpio2         C[9]
  93 #define half            C[10]
  94 #define pio2_1          C[11]
  95 #define pio2_t          C[12]
  96 
  97 void
  98 sincosf(float x, float *s, float *c)
  99 {
 100         double y, z, w;
 101         float f, g;
 102         int n, ix, hx, hy;
 103         volatile int i __unused;
 104 
 105         hx = *((int *)&x);
 106         ix = hx & 0x7fffffff;
 107 
 108         y = (double)x;
 109 
 110         if (ix <= 0x4016cbe4) {                      /* |x| < 3*pi/4 */
 111                 if (ix <= 0x3f490fdb) {              /* |x| < pi/4 */
 112                         if (ix <= 0x39800000) {      /* |x| <= 2**-12 */
 113                                 i = (int)y;
 114 #ifdef lint
 115                                 i = i;
 116 #endif
 117                                 *s = x;
 118                                 *c = 1.0f;
 119                                 return;
 120                         }
 121 
 122                         z = y * y;
 123                         *s = (float)((y * (S0 + z * S1)) * (S2 + z * (S3 + z)));
 124                         *c = (float)(((C0 + z * C1) + (z * z) * C2) * (C3 + z *
 125                             (C4 + z)));
 126                 } else if (hx > 0) {
 127                         y = (y - pio2_1) - pio2_t;
 128                         z = y * y;
 129                         *s = (float)(((C0 + z * C1) + (z * z) * C2) * (C3 + z *
 130                             (C4 + z)));
 131                         *c = (float)-((y * (S0 + z * S1)) * (S2 + z * (S3 +
 132                             z)));
 133                 } else {
 134                         y = (y + pio2_1) + pio2_t;
 135                         z = y * y;
 136                         *s = (float)-(((C0 + z * C1) + (z * z) * C2) * (C3 + z *
 137                             (C4 + z)));
 138                         *c = (float)((y * (S0 + z * S1)) * (S2 + z * (S3 + z)));
 139                 }
 140 
 141                 return;
 142         } else if (ix <= 0x49c90fdb) {       /* |x| < 2^19*pi */
 143 #if defined(__i386) && !defined(__amd64)
 144                 int rp;
 145 
 146                 rp = __swapRP(fp_extended);
 147 #endif
 148                 w = y * invpio2;
 149 
 150                 if (hx < 0)
 151                         n = (int)(w - half);
 152                 else
 153                         n = (int)(w + half);
 154 
 155                 y = (y - n * pio2_1) - n * pio2_t;
 156 #if defined(__i386) && !defined(__amd64)
 157                 if (rp != fp_extended)
 158                         (void) __swapRP(rp);
 159 #endif
 160         } else {
 161                 if (ix >= 0x7f800000) {
 162                         *s = *c = x / x;
 163                         return;
 164                 }
 165 
 166                 hy = ((int *)&y)[HIWORD];
 167                 n = ((hy >> 20) & 0x7ff) - 1046;
 168                 ((int *)&w)[HIWORD] = (hy & 0xfffff) | 0x41600000;
 169                 ((int *)&w)[LOWORD] = ((int *)&y)[LOWORD];
 170                 n = __rem_pio2m(&w, &y, n, 1, 0, _TBL_ipio2_inf);
 171 
 172                 if (hy < 0) {
 173                         y = -y;
 174                         n = -n;
 175                 }
 176         }
 177 
 178         z = y * y;
 179         f = (float)((y * (S0 + z * S1)) * (S2 + z * (S3 + z)));
 180         g = (float)(((C0 + z * C1) + (z * z) * C2) * (C3 + z * (C4 + z)));
 181 
 182         if (n & 2) {
 183                 f = -f;
 184                 g = -g;
 185         }
 186 
 187         if (n & 1) {
 188                 *s = g;
 189                 *c = -f;
 190         } else {
 191                 *s = f;
 192                 *c = g;
 193         }
 194 }