sysdeps/powerpc/fpu/e_rem_pio2f.c

   1 /* e_rem_pio2f.c -- float version of e_rem_pio2.c
   2    Copyright (C) 2011-2014 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Library General Public License as
   8    published by the Free Software Foundation; either version 2 of the
   9    License, or (at your option) any later version.
  10
  11    The GNU C Library 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 GNU
  14    Library General Public License for more details.
  15
  16    You should have received a copy of the GNU Library General Public
  17    License along with the GNU C Library; see the file COPYING.LIB.  If
  18    not, see <http://www.gnu.org/licenses/>.  */
  19
  20 #include <math.h>
  21
  22 #include <math_private.h>
  23 #include "s_float_bitwise.h"
  24
  25 /* defined in sysdeps/powerpc/fpu/k_rem_pio2f.c */
  26 int __fp_kernel_rem_pio2f (float *x, float *y, float e0, int32_t nx);
  27
  28 /* __ieee754_rem_pio2f(x,y)
  29  *
  30  * return the remainder of x rem pi/2 in y[0]+y[1]
  31  */
  32
  33 static const float npio2_hw[] = {
  34   1.57077026e+00, 3.14154053e+00, 4.71228027e+00, 6.28308105e+00,
  35   7.85388184e+00, 9.42456055e+00, 1.09953613e+01, 1.25661621e+01,
  36   1.41369629e+01, 1.57077637e+01, 1.72783203e+01, 1.88491211e+01,
  37   2.04199219e+01, 2.19907227e+01, 2.35615234e+01, 2.51323242e+01,
  38   2.67031250e+01, 2.82739258e+01, 2.98447266e+01, 3.14155273e+01,
  39   3.29863281e+01, 3.45566406e+01, 3.61279297e+01, 3.76982422e+01,
  40   3.92695312e+01, 4.08398438e+01, 4.24111328e+01, 4.39814453e+01,
  41   4.55527344e+01, 4.71230469e+01, 4.86943359e+01, 5.02646484e+01
  42 };
  43
  44
  45 static const float zero  = 0.0000000000e+00;
  46 static const float two8  = 2.5600000000e+02;
  47
  48 static const float half    = 5.0000000000e-01;
  49 static const float invpio2 = 6.3661980629e-01;
  50 static const float pio2_1  = 1.5707855225e+00;
  51 static const float pio2_1t = 1.0804334124e-05;
  52 static const float pio2_2  = 1.0804273188e-05;
  53 static const float pio2_2t = 6.0770999344e-11;
  54 static const float pio2_3  = 6.0770943833e-11;
  55 static const float pio2_3t = 6.1232342629e-17;
  56
  57 static const float pio4     = 7.8539801e-01;
  58 static const float pio3_4   = 2.3561945e+00;
  59 static const float pio2_24b = 1.5707951e+00;
  60 static const float pio2_2e7 = 2.0106054e+02;
  61
  62
  63 int32_t
  64 __ieee754_rem_pio2f (float x, float *y)
  65 {
  66   float ax, z, n, r, w, t, e0;
  67   float tx[3];
  68   int32_t i, nx;
  69
  70   ax = __builtin_fabsf (x);
  71   if (ax <= pio4)
  72     {
  73       y[0] = x;
  74       y[1] = 0;
  75       return 0;
  76     }
  77   if (ax < pio3_4)
  78     {
  79       if (x > 0)
  80         {
  81           z = x - pio2_1;
  82           if (!__float_and_test28 (ax, pio2_24b))
  83             {
  84               y[0] = z - pio2_1t;
  85               y[1] = (z - y[0]) - pio2_1t;
  86             }
  87           else
  88             {
  89               z -= pio2_2;
  90               y[0] = z - pio2_2t;
  91               y[1] = (z - y[0]) - pio2_2t;
  92             }
  93           return 1;
  94         }
  95       else
  96         {
  97           z = x + pio2_1;
  98           if (!__float_and_test28 (ax, pio2_24b))
  99             {
 100               y[0] = z + pio2_1t;
 101               y[1] = (z - y[0]) + pio2_1t;
 102             }
 103           else
 104             {
 105               z += pio2_2;
 106               y[0] = z + pio2_2t;
 107               y[1] = (z - y[0]) + pio2_2t;
 108             }
 109           return -1;
 110         }
 111     }
 112   if (ax <= pio2_2e7)
 113     {
 114       n = __floorf (ax * invpio2 + half);
 115       i = (int32_t) n;
 116       r = ax - n * pio2_1;
 117       w = n * pio2_1t;          /* 1st round good to 40 bit */
 118       if (i < 32 && !__float_and_test24 (ax, npio2_hw[i - 1]))
 119         {
 120           y[0] = r - w;
 121         }
 122       else
 123         {
 124           float i, j;
 125           j = __float_and8 (ax);
 126           y[0] = r - w;
 127           i = __float_and8 (y[0]);
 128           if (j / i > 256.0 || j / i < 3.9062500e-3)
 129             {                   /* 2nd iterations needed, good to 57 */
 130               t = r;
 131               w = n * pio2_2;
 132               r = t - w;
 133               w = n * pio2_2t - ((t - r) - w);
 134               y[0] = r - w;
 135               i = __float_and8 (y[0]);
 136               if (j / i > 33554432 || j / i < 2.9802322e-8)
 137                 {               /* 3rd iteration needed, 74 bits acc */
 138                   t = r;
 139                   w = n * pio2_3;
 140                   r = t - w;
 141                   w = n * pio2_3t - ((t - r) - w);
 142                   y[0] = r - w;
 143                 }
 144             }
 145         }
 146       y[1] = (r - y[0]) - w;
 147       if (x < 0)
 148         {
 149           y[0] = -y[0];
 150           y[1] = -y[1];
 151           return -i;
 152         }
 153       else
 154         {
 155           return i;
 156         }
 157     }
 158
 159   /* all other (large) arguments */
 160   if (isnanf (x) || isinff (x))
 161     {
 162       y[0] = y[1] = x - x;
 163       return 0;
 164     }
 165
 166   /* set z = scalbn(|x|,ilogb(x)-7) */
 167   e0 = __float_and8 (ax / 128.0);
 168   z = ax / e0;
 169
 170   tx[0] = __floorf (z);
 171   z = (z - tx[0]) * two8;
 172   tx[1] = __floorf (z);
 173   z = (z - tx[1]) * two8;
 174   tx[2] = __floorf (z);
 175
 176   nx = 3;
 177   while (tx[nx - 1] == zero)
 178     nx--;
 179
 180   i = __fp_kernel_rem_pio2f (tx, y, e0, nx);
 181   if (x < 0)
 182     {
 183       y[0] = -y[0];
 184       y[1] = -y[1];
 185       return -i;
 186     }
 187   return i;
 188 }