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 __cosf = cosf
  32 
  33 /*
  34  * See sincosf.c
  35  */
  36 
  37 #include "libm.h"
  38 
  39 extern const int _TBL_ipio2_inf[];
  40 extern int __rem_pio2m(double *, double *, int, int, int, const int *);
  41 
  42 #if defined(__i386) && !defined(__amd64)
  43 extern int __swapRP(int);
  44 #endif
  45 
  46 static const double C[] = {
  47         1.85735322054308378716204874632872525989806770558e-0003,
  48         -1.95035094218403635082921458859320791358115801259e-0004,
  49         5.38400550766074785970952495168558701485841707252e+0002,
  50         -3.31975110777873728964197739157371509422022905947e+0001,
  51         1.09349482127188401868272000389539985058873853699e-0003,
  52         -5.03324285989964979398034700054920226866107675091e-0004,
  53         2.43792880266971107750418061559602239831538067410e-0005,
  54         9.14499072605666582228127405245558035523741471271e+0002,
  55         -3.63151270591815439197122504991683846785293207730e+0001,
  56         0.636619772367581343075535,     /* 2^ -1  * 1.45F306DC9C883 */
  57         0.5,
  58         1.570796326734125614166,        /* 2^  0  * 1.921FB54400000 */
  59         6.077100506506192601475e-11,    /* 2^-34  * 1.0B4611A626331 */
  60 };
  61 
  62 #define S0              C[0]
  63 #define S1              C[1]
  64 #define S2              C[2]
  65 #define S3              C[3]
  66 #define C0              C[4]
  67 #define C1              C[5]
  68 #define C2              C[6]
  69 #define C3              C[7]
  70 #define C4              C[8]
  71 #define invpio2         C[9]
  72 #define half            C[10]
  73 #define pio2_1          C[11]
  74 #define pio2_t          C[12]
  75 
  76 float
  77 cosf(float x)
  78 {
  79         double y, z, w;
  80         float f;
  81         int n, ix, hx, hy;
  82         volatile int i __unused;
  83 
  84         hx = *((int *)&x);
  85         ix = hx & 0x7fffffff;
  86 
  87         y = (double)x;
  88 
  89         if (ix <= 0x4016cbe4) {                      /* |x| < 3*pi/4 */
  90                 if (ix <= 0x3f490fdb) {              /* |x| < pi/4 */
  91                         if (ix <= 0x39800000) {      /* |x| <= 2**-12 */
  92                                 i = (int)y;
  93 #ifdef lint
  94                                 i = i;
  95 #endif
  96                                 return (1.0f);
  97                         }
  98 
  99                         z = y * y;
 100                         return ((float)(((C0 + z * C1) + (z * z) * C2) * (C3 +
 101                             z * (C4 + z))));
 102                 } else if (hx > 0) {
 103                         y = (y - pio2_1) - pio2_t;
 104                         z = y * y;
 105                         return ((float)-((y * (S0 + z * S1)) * (S2 + z * (S3 +
 106                             z))));
 107                 } else {
 108                         y = (y + pio2_1) + pio2_t;
 109                         z = y * y;
 110                         return ((float)((y * (S0 + z * S1)) * (S2 + z * (S3 +
 111                             z))));
 112                 }
 113         } else if (ix <= 0x49c90fdb) {       /* |x| < 2^19*pi */
 114 #if defined(__i386) && !defined(__amd64)
 115                 int rp;
 116 
 117                 rp = __swapRP(fp_extended);
 118 #endif
 119                 w = y * invpio2;
 120 
 121                 if (hx < 0)
 122                         n = (int)(w - half);
 123                 else
 124                         n = (int)(w + half);
 125 
 126                 y = (y - n * pio2_1) - n * pio2_t;
 127                 n++;
 128 #if defined(__i386) && !defined(__amd64)
 129                 if (rp != fp_extended)
 130                         (void) __swapRP(rp);
 131 #endif
 132         } else {
 133                 if (ix >= 0x7f800000)
 134                         return (x / x);         /* cos(Inf or NaN) is NaN */
 135 
 136                 hy = ((int *)&y)[HIWORD];
 137                 n = ((hy >> 20) & 0x7ff) - 1046;
 138                 ((int *)&w)[HIWORD] = (hy & 0xfffff) | 0x41600000;
 139                 ((int *)&w)[LOWORD] = ((int *)&y)[LOWORD];
 140                 n = __rem_pio2m(&w, &y, n, 1, 0, _TBL_ipio2_inf) + 1;
 141         }
 142 
 143         if (n & 1) {
 144                 /* compute cos y */
 145                 z = y * y;
 146                 f = (float)(((C0 + z * C1) + (z * z) * C2) * (C3 + z * (C4 +
 147                     z)));
 148         } else {
 149                 /* compute sin y */
 150                 z = y * y;
 151                 f = (float)((y * (S0 + z * S1)) * (S2 + z * (S3 + z)));
 152         }
 153 
 154         return ((n & 2) ? -f : f);
 155 }