1*f504f610SAugustin Cavalier /*
2*f504f610SAugustin Cavalier * Single-precision log function.
3*f504f610SAugustin Cavalier *
4*f504f610SAugustin Cavalier * Copyright (c) 2017-2018, Arm Limited.
5*f504f610SAugustin Cavalier * SPDX-License-Identifier: MIT
6*f504f610SAugustin Cavalier */
7*f504f610SAugustin Cavalier
8*f504f610SAugustin Cavalier #include <math.h>
9*f504f610SAugustin Cavalier #include <stdint.h>
10*f504f610SAugustin Cavalier #include "libm.h"
11*f504f610SAugustin Cavalier #include "logf_data.h"
12*f504f610SAugustin Cavalier
13*f504f610SAugustin Cavalier /*
14*f504f610SAugustin Cavalier LOGF_TABLE_BITS = 4
15*f504f610SAugustin Cavalier LOGF_POLY_ORDER = 4
16*f504f610SAugustin Cavalier
17*f504f610SAugustin Cavalier ULP error: 0.818 (nearest rounding.)
18*f504f610SAugustin Cavalier Relative error: 1.957 * 2^-26 (before rounding.)
19*f504f610SAugustin Cavalier */
20*f504f610SAugustin Cavalier
21*f504f610SAugustin Cavalier #define T __logf_data.tab
22*f504f610SAugustin Cavalier #define A __logf_data.poly
23*f504f610SAugustin Cavalier #define Ln2 __logf_data.ln2
24*f504f610SAugustin Cavalier #define N (1 << LOGF_TABLE_BITS)
25*f504f610SAugustin Cavalier #define OFF 0x3f330000
26*f504f610SAugustin Cavalier
logf(float x)27*f504f610SAugustin Cavalier float logf(float x)
28*f504f610SAugustin Cavalier {
29*f504f610SAugustin Cavalier double_t z, r, r2, y, y0, invc, logc;
30*f504f610SAugustin Cavalier uint32_t ix, iz, tmp;
31*f504f610SAugustin Cavalier int k, i;
32*f504f610SAugustin Cavalier
33*f504f610SAugustin Cavalier ix = asuint(x);
34*f504f610SAugustin Cavalier /* Fix sign of zero with downward rounding when x==1. */
35*f504f610SAugustin Cavalier if (WANT_ROUNDING && predict_false(ix == 0x3f800000))
36*f504f610SAugustin Cavalier return 0;
37*f504f610SAugustin Cavalier if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
38*f504f610SAugustin Cavalier /* x < 0x1p-126 or inf or nan. */
39*f504f610SAugustin Cavalier if (ix * 2 == 0)
40*f504f610SAugustin Cavalier return __math_divzerof(1);
41*f504f610SAugustin Cavalier if (ix == 0x7f800000) /* log(inf) == inf. */
42*f504f610SAugustin Cavalier return x;
43*f504f610SAugustin Cavalier if ((ix & 0x80000000) || ix * 2 >= 0xff000000)
44*f504f610SAugustin Cavalier return __math_invalidf(x);
45*f504f610SAugustin Cavalier /* x is subnormal, normalize it. */
46*f504f610SAugustin Cavalier ix = asuint(x * 0x1p23f);
47*f504f610SAugustin Cavalier ix -= 23 << 23;
48*f504f610SAugustin Cavalier }
49*f504f610SAugustin Cavalier
50*f504f610SAugustin Cavalier /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
51*f504f610SAugustin Cavalier The range is split into N subintervals.
52*f504f610SAugustin Cavalier The ith subinterval contains z and c is near its center. */
53*f504f610SAugustin Cavalier tmp = ix - OFF;
54*f504f610SAugustin Cavalier i = (tmp >> (23 - LOGF_TABLE_BITS)) % N;
55*f504f610SAugustin Cavalier k = (int32_t)tmp >> 23; /* arithmetic shift */
56*f504f610SAugustin Cavalier iz = ix - (tmp & 0x1ff << 23);
57*f504f610SAugustin Cavalier invc = T[i].invc;
58*f504f610SAugustin Cavalier logc = T[i].logc;
59*f504f610SAugustin Cavalier z = (double_t)asfloat(iz);
60*f504f610SAugustin Cavalier
61*f504f610SAugustin Cavalier /* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */
62*f504f610SAugustin Cavalier r = z * invc - 1;
63*f504f610SAugustin Cavalier y0 = logc + (double_t)k * Ln2;
64*f504f610SAugustin Cavalier
65*f504f610SAugustin Cavalier /* Pipelined polynomial evaluation to approximate log1p(r). */
66*f504f610SAugustin Cavalier r2 = r * r;
67*f504f610SAugustin Cavalier y = A[1] * r + A[2];
68*f504f610SAugustin Cavalier y = A[0] * r2 + y;
69*f504f610SAugustin Cavalier y = y * r2 + (y0 + r);
70*f504f610SAugustin Cavalier return eval_as_float(y);
71*f504f610SAugustin Cavalier }
72