]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/flt-32/e_logf.c
1 /* Single-precision log function.
2 Copyright (C) 2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
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.
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.
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/>. */
21 #include <shlib-compat.h>
22 #include "math_config.h"
28 ULP error: 0.818 (nearest rounding.)
29 Relative error: 1.957 * 2^-26 (before rounding.)
32 #define T __logf_data.tab
33 #define A __logf_data.poly
34 #define Ln2 __logf_data.ln2
35 #define N (1 << LOGF_TABLE_BITS)
36 #define OFF 0x3f330000
41 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
42 double_t z
, r
, r2
, y
, y0
, invc
, logc
;
48 /* Fix sign of zero with downward rounding when x==1. */
49 if (__glibc_unlikely (ix
== 0x3f800000))
52 if (__glibc_unlikely (ix
- 0x00800000 >= 0x7f800000 - 0x00800000))
54 /* x < 0x1p-126 or inf or nan. */
56 return __math_divzerof (1);
57 if (ix
== 0x7f800000) /* log(inf) == inf. */
59 if ((ix
& 0x80000000) || ix
* 2 >= 0xff000000)
60 return __math_invalidf (x
);
61 /* x is subnormal, normalize it. */
62 ix
= asuint (x
* 0x1p
23f
);
66 /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
67 The range is split into N subintervals.
68 The ith subinterval contains z and c is near its center. */
70 i
= (tmp
>> (23 - LOGF_TABLE_BITS
)) % N
;
71 k
= (int32_t) tmp
>> 23; /* arithmetic shift */
72 iz
= ix
- (tmp
& 0x1ff << 23);
75 z
= (double_t
) asfloat (iz
);
77 /* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */
79 y0
= logc
+ (double_t
) k
* Ln2
;
81 /* Pipelined polynomial evaluation to approximate log1p(r). */
85 y
= y
* r2
+ (y0
+ r
);
89 strong_alias (__logf
, __ieee754_logf
)
90 strong_alias (__logf
, __logf_finite
)
91 versioned_symbol (libm
, __logf
, logf
, GLIBC_2_27
);