]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/powerpc/fpu/e_rem_pio2f.c
1 /* e_rem_pio2f.c -- float version of e_rem_pio2.c
2 Copyright (C) 2011-2018 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011
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.
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.
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/>. */
22 #include <math_private.h>
23 #include "s_float_bitwise.h"
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
);
28 /* __ieee754_rem_pio2f(x,y)
30 * return the remainder of x rem pi/2 in y[0]+y[1]
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
45 static const float zero
= 0.0000000000e+00;
46 static const float two8
= 2.5600000000e+02;
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;
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;
64 __ieee754_rem_pio2f (float x
, float *y
)
66 float ax
, z
, n
, r
, w
, t
, e0
;
70 ax
= __builtin_fabsf (x
);
82 if (!__float_and_test28 (ax
, pio2_24b
))
85 y
[1] = (z
- y
[0]) - pio2_1t
;
91 y
[1] = (z
- y
[0]) - pio2_2t
;
98 if (!__float_and_test28 (ax
, pio2_24b
))
101 y
[1] = (z
- y
[0]) + pio2_1t
;
107 y
[1] = (z
- y
[0]) + pio2_2t
;
114 n
= __floorf (ax
* invpio2
+ half
);
117 w
= n
* pio2_1t
; /* 1st round good to 40 bit */
118 if (i
< 32 && !__float_and_test24 (ax
, npio2_hw
[i
- 1]))
125 j
= __float_and8 (ax
);
127 i
= __float_and8 (y
[0]);
128 if (j
/ i
> 256.0 || j
/ i
< 3.9062500e-3)
129 { /* 2nd iterations needed, good to 57 */
133 w
= n
* pio2_2t
- ((t
- 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 */
141 w
= n
* pio2_3t
- ((t
- r
) - w
);
146 y
[1] = (r
- y
[0]) - w
;
159 /* all other (large) arguments */
160 if (isnanf (x
) || isinff (x
))
166 /* set z = scalbn(|x|,ilogb(x)-7) */
167 e0
= __float_and8 (ax
/ 128.0);
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
);
177 while (tx
[nx
- 1] == zero
)
180 i
= __fp_kernel_rem_pio2f (tx
, y
, e0
, nx
);