]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/powerpc/power7/fpu/s_logb.c
Update copyright dates with scripts/update-copyrights.
[thirdparty/glibc.git] / sysdeps / powerpc / power7 / fpu / s_logb.c
CommitLineData
777b1eea 1/* logb(). PowerPC/POWER7 version.
b168057a 2 Copyright (C) 2012-2015 Free Software Foundation, Inc.
777b1eea
AZ
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <http://www.gnu.org/licenses/>. */
18
45470df3 19#include <math_ldbl_opt.h>
777b1eea
AZ
20
21/* This implementation avoids FP to INT conversions by using VSX
22 bitwise instructions over FP values. */
23
24static const double two1div52 = 2.220446049250313e-16; /* 1/2**52 */
25static const double two10m1 = -1023.0; /* 2**10 -1 */
26
27/* FP mask to extract the exponent. */
28static const union {
29 unsigned long long mask;
30 double d;
31} mask = { 0x7ff0000000000000ULL };
32
33double
34__logb (double x)
35{
36 double ret;
37
38 if (__builtin_expect (x == 0.0, 0))
39 /* Raise FE_DIVBYZERO and return -HUGE_VAL[LF]. */
40 return -1.0 / __builtin_fabs (x);
41
42 /* ret = x & 0x7ff0000000000000; */
43 asm (
44 "xxland %x0,%x1,%x2\n"
45 "fcfid %0,%0"
46 : "=f" (ret)
47 : "f" (x), "f" (mask.d));
48 /* ret = (ret >> 52) - 1023.0; */
49 ret = (ret * two1div52) + two10m1;
50 if (__builtin_expect (ret > -two10m1, 0))
51 /* Multiplication is used to set logb (+-INF) = INF. */
52 return (x * x);
53 else if (__builtin_expect (ret == two10m1, 0))
54 {
55 /* POSIX specifies that denormal numbers are treated as
56 though they were normalized. */
57 int32_t lx, ix;
25dbcb27 58 int ma;
777b1eea 59
25dbcb27
AS
60 EXTRACT_WORDS (ix, lx, x);
61 if (ix == 0)
62 ma = __builtin_clz (lx) + 32;
63 else
64 ma = __builtin_clz (ix);
65 return (double) (-1023 - (ma - 12));
777b1eea
AZ
66 }
67 /* Test to avoid logb_downward (0.0) == -0.0. */
68 return ret == -0.0 ? 0.0 : ret;
69}
777b1eea 70weak_alias (__logb, logb)
777b1eea
AZ
71#ifdef NO_LONG_DOUBLE
72strong_alias (__logb, __logbl)
73weak_alias (__logb, logbl)
74#endif
45470df3
AZ
75
76#if LONG_DOUBLE_COMPAT (libm, GLIBC_2_0)
77compat_symbol (libm, logb, logbl, GLIBC_2_0);
78#endif