Back to index

glibc  2.9
k_sincosl.c
Go to the documentation of this file.
00001 /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>.
00002    Copyright (C) 1999,2004,2006 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
00005 
00006    The GNU C Library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Lesser General Public
00008    License as published by the Free Software Foundation; either
00009    version 2.1 of the License, or (at your option) any later version.
00010 
00011    The GNU C Library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Lesser General Public License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public
00017    License along with the GNU C Library; if not, write to the Free
00018    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00019    02111-1307 USA.  */
00020 
00021 #include "math.h"
00022 #include "math_private.h"
00023 
00024 static const long double c[] = {
00025 #define ONE c[0]
00026  1.00000000000000000000000000000000000E+00L, /* 3fff0000000000000000000000000000 */
00027 
00028 /* cos x ~ ONE + x^2 ( SCOS1 + SCOS2 * x^2 + ... + SCOS4 * x^6 + SCOS5 * x^8 )
00029    x in <0,1/256>  */
00030 #define SCOS1 c[1]
00031 #define SCOS2 c[2]
00032 #define SCOS3 c[3]
00033 #define SCOS4 c[4]
00034 #define SCOS5 c[5]
00035 -5.00000000000000000000000000000000000E-01L, /* bffe0000000000000000000000000000 */
00036  4.16666666666666666666666666556146073E-02L, /* 3ffa5555555555555555555555395023 */
00037 -1.38888888888888888888309442601939728E-03L, /* bff56c16c16c16c16c16a566e42c0375 */
00038  2.48015873015862382987049502531095061E-05L, /* 3fefa01a01a019ee02dcf7da2d6d5444 */
00039 -2.75573112601362126593516899592158083E-07L, /* bfe927e4f5dce637cb0b54908754bde0 */
00040 
00041 /* cos x ~ ONE + x^2 ( COS1 + COS2 * x^2 + ... + COS7 * x^12 + COS8 * x^14 )
00042    x in <0,0.1484375>  */
00043 #define COS1 c[6]
00044 #define COS2 c[7]
00045 #define COS3 c[8]
00046 #define COS4 c[9]
00047 #define COS5 c[10]
00048 #define COS6 c[11]
00049 #define COS7 c[12]
00050 #define COS8 c[13]
00051 -4.99999999999999999999999999999999759E-01L, /* bffdfffffffffffffffffffffffffffb */
00052  4.16666666666666666666666666651287795E-02L, /* 3ffa5555555555555555555555516f30 */
00053 -1.38888888888888888888888742314300284E-03L, /* bff56c16c16c16c16c16c16a463dfd0d */
00054  2.48015873015873015867694002851118210E-05L, /* 3fefa01a01a01a01a0195cebe6f3d3a5 */
00055 -2.75573192239858811636614709689300351E-07L, /* bfe927e4fb7789f5aa8142a22044b51f */
00056  2.08767569877762248667431926878073669E-09L, /* 3fe21eed8eff881d1e9262d7adff4373 */
00057 -1.14707451049343817400420280514614892E-11L, /* bfda9397496922a9601ed3d4ca48944b */
00058  4.77810092804389587579843296923533297E-14L, /* 3fd2ae5f8197cbcdcaf7c3fb4523414c */
00059 
00060 /* sin x ~ ONE * x + x^3 ( SSIN1 + SSIN2 * x^2 + ... + SSIN4 * x^6 + SSIN5 * x^8 )
00061    x in <0,1/256>  */
00062 #define SSIN1 c[14]
00063 #define SSIN2 c[15]
00064 #define SSIN3 c[16]
00065 #define SSIN4 c[17]
00066 #define SSIN5 c[18]
00067 -1.66666666666666666666666666666666659E-01L, /* bffc5555555555555555555555555555 */
00068  8.33333333333333333333333333146298442E-03L, /* 3ff81111111111111111111110fe195d */
00069 -1.98412698412698412697726277416810661E-04L, /* bff2a01a01a01a01a019e7121e080d88 */
00070  2.75573192239848624174178393552189149E-06L, /* 3fec71de3a556c640c6aaa51aa02ab41 */
00071 -2.50521016467996193495359189395805639E-08L, /* bfe5ae644ee90c47dc71839de75b2787 */
00072 
00073 /* sin x ~ ONE * x + x^3 ( SIN1 + SIN2 * x^2 + ... + SIN7 * x^12 + SIN8 * x^14 )
00074    x in <0,0.1484375>  */
00075 #define SIN1 c[19]
00076 #define SIN2 c[20]
00077 #define SIN3 c[21]
00078 #define SIN4 c[22]
00079 #define SIN5 c[23]
00080 #define SIN6 c[24]
00081 #define SIN7 c[25]
00082 #define SIN8 c[26]
00083 -1.66666666666666666666666666666666538e-01L, /* bffc5555555555555555555555555550 */
00084  8.33333333333333333333333333307532934e-03L, /* 3ff811111111111111111111110e7340 */
00085 -1.98412698412698412698412534478712057e-04L, /* bff2a01a01a01a01a01a019e7a626296 */
00086  2.75573192239858906520896496653095890e-06L, /* 3fec71de3a556c7338fa38527474b8f5 */
00087 -2.50521083854417116999224301266655662e-08L, /* bfe5ae64567f544e16c7de65c2ea551f */
00088  1.60590438367608957516841576404938118e-10L, /* 3fde6124613a811480538a9a41957115 */
00089 -7.64716343504264506714019494041582610e-13L, /* bfd6ae7f3d5aef30c7bc660b060ef365 */
00090  2.81068754939739570236322404393398135e-15L, /* 3fce9510115aabf87aceb2022a9a9180 */
00091 };
00092 
00093 #define SINCOSL_COS_HI 0
00094 #define SINCOSL_COS_LO 1
00095 #define SINCOSL_SIN_HI 2
00096 #define SINCOSL_SIN_LO 3
00097 extern const long double __sincosl_table[];
00098 
00099 void
00100 __kernel_sincosl(long double x, long double y, long double *sinx, long double *cosx, int iy)
00101 {
00102   long double h, l, z, sin_l, cos_l_m1;
00103   int64_t ix;
00104   u_int32_t tix, hix, index;
00105   GET_LDOUBLE_MSW64 (ix, x);
00106   tix = ((u_int64_t)ix) >> 32;
00107   tix &= ~0x80000000;                     /* tix = |x|'s high 32 bits */
00108   if (tix < 0x3fc30000)                   /* |x| < 0.1484375 */
00109     {
00110       /* Argument is small enough to approximate it by a Chebyshev
00111         polynomial of degree 16(17).  */
00112       if (tix < 0x3c600000)        /* |x| < 2^-57 */
00113        if (!((int)x))                     /* generate inexact */
00114          {
00115            *sinx = x;
00116            *cosx = ONE;
00117            return;
00118          }
00119       z = x * x;
00120       *sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
00121                      z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
00122       *cosx = ONE + (z*(COS1+z*(COS2+z*(COS3+z*(COS4+
00123                    z*(COS5+z*(COS6+z*(COS7+z*COS8))))))));
00124     }
00125   else
00126     {
00127       /* So that we don't have to use too large polynomial,  we find
00128         l and h such that x = l + h,  where fabsl(l) <= 1.0/256 with 83
00129         possible values for h.  We look up cosl(h) and sinl(h) in
00130         pre-computed tables,  compute cosl(l) and sinl(l) using a
00131         Chebyshev polynomial of degree 10(11) and compute
00132         sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and
00133         cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l).  */
00134       int six = tix;
00135       tix = ((six - 0x3ff00000) >> 4) + 0x3fff0000;
00136       index = 0x3ffe - (tix >> 16);
00137       hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
00138       x = fabsl (x);
00139       switch (index)
00140        {
00141        case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
00142        case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
00143        default:
00144        case 2: index = (hix - 0x3ffc3000) >> 10; break;
00145        }
00146       hix = (hix << 4) & 0x3fffffff;
00147 /*
00148     The following should work for double but generates the wrong index.
00149     For now the code above converts double to ieee extended to compute
00150     the index back to double for the h value. 
00151     
00152 
00153       index = 0x3fe - (tix >> 20);
00154       hix = (tix + (0x2000 << index)) & (0xffffc000 << index);
00155       x = fabsl (x);
00156       switch (index)
00157        {
00158        case 0: index = ((45 << 14) + hix - 0x3fe00000) >> 12; break;
00159        case 1: index = ((13 << 15) + hix - 0x3fd00000) >> 13; break;
00160        default:
00161        case 2: index = (hix - 0x3fc30000) >> 14; break;
00162        }
00163 */
00164       SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
00165       if (iy)
00166        l = y - (h - x);
00167       else
00168        l = x - h;
00169       z = l * l;
00170       sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5)))));
00171       cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5))));
00172       z = __sincosl_table [index + SINCOSL_SIN_HI]
00173          + (__sincosl_table [index + SINCOSL_SIN_LO]
00174             + (__sincosl_table [index + SINCOSL_SIN_HI] * cos_l_m1)
00175             + (__sincosl_table [index + SINCOSL_COS_HI] * sin_l));
00176       *sinx = (ix < 0) ? -z : z;
00177       *cosx = __sincosl_table [index + SINCOSL_COS_HI]
00178              + (__sincosl_table [index + SINCOSL_COS_LO]
00179                - (__sincosl_table [index + SINCOSL_SIN_HI] * sin_l
00180                   - __sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1));
00181     }
00182 }