Print this page
11210 libm should be cstyle(1ONBLD) clean
   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 


  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 }
   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 


 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 }