Adjust thresholds in Bessel function implementations (bug 14469).
[platform/upstream/glibc.git] / math / e_scalb_template.c
1 /* Multiply by integral power of radix.
2
3    Copyright (C) 2011-2020 Free Software Foundation, Inc.
4
5    This file is part of the GNU C Library.
6
7    The GNU C Library is free software; you can redistribute it and/or
8    modify it under the terms of the GNU Lesser General Public
9    License as published by the Free Software Foundation; either
10    version 2.1 of the License, or (at your option) any later version.
11
12    The GNU C Library is distributed in the hope that it will be useful,
13    but WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15    Lesser General Public License for more details.
16
17    You should have received a copy of the GNU Lesser General Public
18    License along with the GNU C Library; if not, see
19    <https://www.gnu.org/licenses/>.  */
20
21 #include <math.h>
22 #include <math_private.h>
23 #include <libm-alias-finite.h>
24
25 static FLOAT
26 __attribute__ ((noinline))
27 invalid_fn (FLOAT x, FLOAT fn)
28 {
29   if (M_SUF (rint) (fn) != fn)
30     return (fn - fn) / (fn - fn);
31   else if (fn > M_LIT (65000.0))
32     return M_SUF (__scalbn) (x, 65000);
33   else
34     return M_SUF (__scalbn) (x,-65000);
35 }
36
37
38 FLOAT
39 M_DECL_FUNC (__ieee754_scalb) (FLOAT x, FLOAT fn)
40 {
41   if (__glibc_unlikely (isnan (x)))
42     return x * fn;
43   if (__glibc_unlikely (!isfinite (fn)))
44     {
45       if (isnan (fn) || fn > M_LIT (0.0))
46         return x * fn;
47       if (x == M_LIT (0.0))
48         return x;
49       return x / -fn;
50     }
51   if (__glibc_unlikely (M_FABS (fn) >= M_LIT (0x1p31)
52                         || (FLOAT) (int) fn != fn))
53     return invalid_fn (x, fn);
54
55   return M_SCALBN (x, (int) fn);
56 }
57 declare_mgen_finite_alias (__ieee754_scalb, __scalb)