sysdeps/ieee754/dbl-64/sincos32.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /****************************************************************/
  20 /*  MODULE_NAME: sincos32.c                                     */
  21 /*                                                              */
  22 /*  FUNCTIONS: ss32                                             */
  23 /*             cc32                                             */
  24 /*             c32                                              */
  25 /*             sin32                                            */
  26 /*             cos32                                            */
  27 /*             mpsin                                            */
  28 /*             mpcos                                            */
  29 /*             mpranred                                         */
  30 /*             mpsin1                                           */
  31 /*             mpcos1                                           */
  32 /*                                                              */
  33 /* FILES NEEDED: endian.h mpa.h sincos32.h                      */
  34 /*               mpa.c                                          */
  35 /*                                                              */
  36 /* Multi Precision sin() and cos() function with p=32  for sin()*/
  37 /* cos() arcsin() and arccos() routines                         */
  38 /* In addition mpranred() routine  performs range  reduction of */
  39 /* a double number x into multi precision number   y,           */
  40 /* such that y=x-n*pi/2, abs(y)<pi/4,  n=0,+-1,+-2,....         */
  41 /****************************************************************/
  42 #include "endian.h"
  43 #include "mpa.h"
  44 #include "sincos32.h"
  45 #include <math_private.h>
  46 #include <stap-probe.h>
  47
  48 #ifndef SECTION
  49 # define SECTION
  50 #endif
  51
  52 /* Compute Multi-Precision sin() function for given p.  Receive Multi Precision
  53    number x and result stored at y.  */
  54 static void
  55 SECTION
  56 ss32 (mp_no *x, mp_no *y, int p)
  57 {
  58   int i;
  59   double a;
  60   mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
  61   for (i = 1; i <= p; i++)
  62     mpk.d[i] = 0;
  63
  64   __sqr (x, &x2, p);
  65   __cpy (&oofac27, &gor, p);
  66   __cpy (&gor, &sum, p);
  67   for (a = 27.0; a > 1.0; a -= 2.0)
  68     {
  69       mpk.d[1] = a * (a - 1.0);
  70       __mul (&gor, &mpk, &mpt1, p);
  71       __cpy (&mpt1, &gor, p);
  72       __mul (&x2, &sum, &mpt1, p);
  73       __sub (&gor, &mpt1, &sum, p);
  74     }
  75   __mul (x, &sum, y, p);
  76 }
  77
  78 /* Compute Multi-Precision cos() function for given p. Receive Multi Precision
  79    number x and result stored at y.  */
  80 static void
  81 SECTION
  82 cc32 (mp_no *x, mp_no *y, int p)
  83 {
  84   int i;
  85   double a;
  86   mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
  87   for (i = 1; i <= p; i++)
  88     mpk.d[i] = 0;
  89
  90   __sqr (x, &x2, p);
  91   mpk.d[1] = 27.0;
  92   __mul (&oofac27, &mpk, &gor, p);
  93   __cpy (&gor, &sum, p);
  94   for (a = 26.0; a > 2.0; a -= 2.0)
  95     {
  96       mpk.d[1] = a * (a - 1.0);
  97       __mul (&gor, &mpk, &mpt1, p);
  98       __cpy (&mpt1, &gor, p);
  99       __mul (&x2, &sum, &mpt1, p);
 100       __sub (&gor, &mpt1, &sum, p);
 101     }
 102   __mul (&x2, &sum, y, p);
 103 }
 104
 105 /* Compute both sin(x), cos(x) as Multi precision numbers.  */
 106 void
 107 SECTION
 108 __c32 (mp_no *x, mp_no *y, mp_no *z, int p)
 109 {
 110   mp_no u, t, t1, t2, c, s;
 111   int i;
 112   __cpy (x, &u, p);
 113   u.e = u.e - 1;
 114   cc32 (&u, &c, p);
 115   ss32 (&u, &s, p);
 116   for (i = 0; i < 24; i++)
 117     {
 118       __mul (&c, &s, &t, p);
 119       __sub (&s, &t, &t1, p);
 120       __add (&t1, &t1, &s, p);
 121       __sub (&mptwo, &c, &t1, p);
 122       __mul (&t1, &c, &t2, p);
 123       __add (&t2, &t2, &c, p);
 124     }
 125   __sub (&mpone, &c, y, p);
 126   __cpy (&s, z, p);
 127 }
 128
 129 /* Receive double x and two double results of sin(x) and return result which is
 130    more accurate, computing sin(x) with multi precision routine c32.  */
 131 double
 132 SECTION
 133 __sin32 (double x, double res, double res1)
 134 {
 135   int p;
 136   mp_no a, b, c;
 137   p = 32;
 138   __dbl_mp (res, &a, p);
 139   __dbl_mp (0.5 * (res1 - res), &b, p);
 140   __add (&a, &b, &c, p);
 141   if (x > 0.8)
 142     {
 143       __sub (&hp, &c, &a, p);
 144       __c32 (&a, &b, &c, p);
 145     }
 146   else
 147     __c32 (&c, &a, &b, p);      /* b=sin(0.5*(res+res1))  */
 148   __dbl_mp (x, &c, p);          /* c = x  */
 149   __sub (&b, &c, &a, p);
 150   /* if a > 0 return min (res, res1), otherwise return max (res, res1).  */
 151   if ((a.d[0] > 0 && res >= res1) || (a.d[0] <= 0 && res <= res1))
 152     res = res1;
 153   LIBC_PROBE (slowasin, 2, &res, &x);
 154   return res;
 155 }
 156
 157 /* Receive double x and two double results of cos(x) and return result which is
 158    more accurate, computing cos(x) with multi precision routine c32.  */
 159 double
 160 SECTION
 161 __cos32 (double x, double res, double res1)
 162 {
 163   int p;
 164   mp_no a, b, c;
 165   p = 32;
 166   __dbl_mp (res, &a, p);
 167   __dbl_mp (0.5 * (res1 - res), &b, p);
 168   __add (&a, &b, &c, p);
 169   if (x > 2.4)
 170     {
 171       __sub (&pi, &c, &a, p);
 172       __c32 (&a, &b, &c, p);
 173       b.d[0] = -b.d[0];
 174     }
 175   else if (x > 0.8)
 176     {
 177       __sub (&hp, &c, &a, p);
 178       __c32 (&a, &c, &b, p);
 179     }
 180   else
 181     __c32 (&c, &b, &a, p);      /* b=cos(0.5*(res+res1))  */
 182   __dbl_mp (x, &c, p);          /* c = x                  */
 183   __sub (&b, &c, &a, p);
 184   /* if a > 0 return max (res, res1), otherwise return min (res, res1).  */
 185   if ((a.d[0] > 0 && res <= res1) || (a.d[0] <= 0 && res >= res1))
 186     res = res1;
 187   LIBC_PROBE (slowacos, 2, &res, &x);
 188   return res;
 189 }
 190
 191 /* Compute sin() of double-length number (X + DX) as Multi Precision number and
 192    return result as double.  If REDUCE_RANGE is true, X is assumed to be the
 193    original input and DX is ignored.  */
 194 double
 195 SECTION
 196 __mpsin (double x, double dx, bool reduce_range)
 197 {
 198   double y;
 199   mp_no a, b, c, s;
 200   int n;
 201   int p = 32;
 202
 203   if (reduce_range)
 204     {
 205       n = __mpranred (x, &a, p);        /* n is 0, 1, 2 or 3.  */
 206       __c32 (&a, &c, &s, p);
 207     }
 208   else
 209     {
 210       n = -1;
 211       __dbl_mp (x, &b, p);
 212       __dbl_mp (dx, &c, p);
 213       __add (&b, &c, &a, p);
 214       if (x > 0.8)
 215         {
 216           __sub (&hp, &a, &b, p);
 217           __c32 (&b, &s, &c, p);
 218         }
 219       else
 220         __c32 (&a, &c, &s, p);  /* b = sin(x+dx)  */
 221     }
 222
 223   /* Convert result based on which quarter of unit circle y is in.  */
 224   switch (n)
 225     {
 226     case 1:
 227       __mp_dbl (&c, &y, p);
 228       break;
 229
 230     case 3:
 231       __mp_dbl (&c, &y, p);
 232       y = -y;
 233       break;
 234
 235     case 2:
 236       __mp_dbl (&s, &y, p);
 237       y = -y;
 238       break;
 239
 240     /* Quadrant not set, so the result must be sin (X + DX), which is also in
 241        S.  */
 242     case 0:
 243     default:
 244       __mp_dbl (&s, &y, p);
 245     }
 246   LIBC_PROBE (slowsin, 3, &x, &dx, &y);
 247   return y;
 248 }
 249
 250 /* Compute cos() of double-length number (X + DX) as Multi Precision number and
 251    return result as double.  If REDUCE_RANGE is true, X is assumed to be the
 252    original input and DX is ignored.  */
 253 double
 254 SECTION
 255 __mpcos (double x, double dx, bool reduce_range)
 256 {
 257   double y;
 258   mp_no a, b, c, s;
 259   int n;
 260   int p = 32;
 261
 262   if (reduce_range)
 263     {
 264       n = __mpranred (x, &a, p);        /* n is 0, 1, 2 or 3.  */
 265       __c32 (&a, &c, &s, p);
 266     }
 267   else
 268     {
 269       n = -1;
 270       __dbl_mp (x, &b, p);
 271       __dbl_mp (dx, &c, p);
 272       __add (&b, &c, &a, p);
 273       if (x > 0.8)
 274         {
 275           __sub (&hp, &a, &b, p);
 276           __c32 (&b, &s, &c, p);
 277         }
 278       else
 279         __c32 (&a, &c, &s, p);  /* a = cos(x+dx)     */
 280     }
 281
 282   /* Convert result based on which quarter of unit circle y is in.  */
 283   switch (n)
 284     {
 285     case 1:
 286       __mp_dbl (&s, &y, p);
 287       y = -y;
 288       break;
 289
 290     case 3:
 291       __mp_dbl (&s, &y, p);
 292       break;
 293
 294     case 2:
 295       __mp_dbl (&c, &y, p);
 296       y = -y;
 297       break;
 298
 299     /* Quadrant not set, so the result must be cos (X + DX), which is also
 300        stored in C.  */
 301     case 0:
 302     default:
 303       __mp_dbl (&c, &y, p);
 304     }
 305   LIBC_PROBE (slowcos, 3, &x, &dx, &y);
 306   return y;
 307 }
 308
 309 /* Perform range reduction of a double number x into multi precision number y,
 310    such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ...
 311    Return int which indicates in which quarter of circle x is.  */
 312 int
 313 SECTION
 314 __mpranred (double x, mp_no *y, int p)
 315 {
 316   number v;
 317   double t, xn;
 318   int i, k, n;
 319   mp_no a, b, c;
 320
 321   if (ABS (x) < 2.8e14)
 322     {
 323       t = (x * hpinv.d + toint.d);
 324       xn = t - toint.d;
 325       v.d = t;
 326       n = v.i[LOW_HALF] & 3;
 327       __dbl_mp (xn, &a, p);
 328       __mul (&a, &hp, &b, p);
 329       __dbl_mp (x, &c, p);
 330       __sub (&c, &b, y, p);
 331       return n;
 332     }
 333   else
 334     {
 335       /* If x is very big more precision required.  */
 336       __dbl_mp (x, &a, p);
 337       a.d[0] = 1.0;
 338       k = a.e - 5;
 339       if (k < 0)
 340         k = 0;
 341       b.e = -k;
 342       b.d[0] = 1.0;
 343       for (i = 0; i < p; i++)
 344         b.d[i + 1] = toverp[i + k];
 345       __mul (&a, &b, &c, p);
 346       t = c.d[c.e];
 347       for (i = 1; i <= p - c.e; i++)
 348         c.d[i] = c.d[i + c.e];
 349       for (i = p + 1 - c.e; i <= p; i++)
 350         c.d[i] = 0;
 351       c.e = 0;
 352       if (c.d[1] >= HALFRAD)
 353         {
 354           t += 1.0;
 355           __sub (&c, &mpone, &b, p);
 356           __mul (&b, &hp, y, p);
 357         }
 358       else
 359         __mul (&c, &hp, y, p);
 360       n = (int) t;
 361       if (x < 0)
 362         {
 363           y->d[0] = -y->d[0];
 364           n = -n;
 365         }
 366       return (n & 3);
 367     }
 368 }