usr/src/lib/libm/common/R/sincosf.c

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