1 /* Double-precision vector (Advanced SIMD) log function.
3 Copyright (C) 2023 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the 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 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <https://www.gnu.org/licenses/>. */
22 static const struct data
25 uint32x4_t special_bound
;
28 uint64x2_t sign_exp_mask
;
30 /* Worst-case error: 1.17 + 0.5 ulp.
31 Rel error: 0x1.6272e588p-56 in [ -0x1.fc1p-9 0x1.009p-8 ]. */
32 .poly
= { V2 (-0x1.ffffffffffff7p
-2), V2 (0x1.55555555170d4p
-2),
33 V2 (-0x1.0000000399c27p
-2), V2 (0x1.999b2e90e94cap
-3),
34 V2 (-0x1.554e550bd501ep
-3) },
35 .ln2
= V2 (0x1.62e42fefa39efp
-1),
36 .min_norm
= V2 (0x0010000000000000),
37 .special_bound
= V4 (0x7fe00000), /* asuint64(inf) - min_norm. */
38 .sign_exp_mask
= V2 (0xfff0000000000000)
41 #define A(i) d->poly[i]
42 #define N (1 << V_LOG_TABLE_BITS)
43 #define IndexMask (N - 1)
44 #define Off v_u64 (0x3fe6900900000000)
52 static inline struct entry
55 /* Since N is a power of 2, n % N = n & (N - 1). */
57 uint64_t i0
= (i
[0] >> (52 - V_LOG_TABLE_BITS
)) & IndexMask
;
58 uint64_t i1
= (i
[1] >> (52 - V_LOG_TABLE_BITS
)) & IndexMask
;
59 float64x2_t e0
= vld1q_f64 (&__v_log_data
.table
[i0
].invc
);
60 float64x2_t e1
= vld1q_f64 (&__v_log_data
.table
[i1
].invc
);
61 e
.invc
= vuzp1q_f64 (e0
, e1
);
62 e
.logc
= vuzp2q_f64 (e0
, e1
);
66 static float64x2_t VPCS_ATTR NOINLINE
67 special_case (float64x2_t x
, float64x2_t y
, float64x2_t hi
, float64x2_t r2
,
70 return v_call_f64 (log
, x
, vfmaq_f64 (hi
, y
, r2
), vmovl_u32 (cmp
));
73 float64x2_t VPCS_ATTR
V_NAME_D1 (log
) (float64x2_t x
)
75 const struct data
*d
= ptr_barrier (&data
);
76 float64x2_t z
, r
, r2
, p
, y
, kd
, hi
;
77 uint64x2_t ix
, iz
, tmp
;
82 ix
= vreinterpretq_u64_f64 (x
);
83 cmp
= vcge_u32 (vsubhn_u64 (ix
, d
->min_norm
),
84 vget_low_u32 (d
->special_bound
));
86 /* x = 2^k z; where z is in range [Off,2*Off) and exact.
87 The range is split into N subintervals.
88 The ith subinterval contains z and c is near its center. */
89 tmp
= vsubq_u64 (ix
, Off
);
90 k
= vshrq_n_s64 (vreinterpretq_s64_u64 (tmp
), 52); /* arithmetic shift. */
91 iz
= vsubq_u64 (ix
, vandq_u64 (tmp
, d
->sign_exp_mask
));
92 z
= vreinterpretq_f64_u64 (iz
);
95 /* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */
96 r
= vfmaq_f64 (v_f64 (-1.0), z
, e
.invc
);
97 kd
= vcvtq_f64_s64 (k
);
99 /* hi = r + log(c) + k*Ln2. */
100 hi
= vfmaq_f64 (vaddq_f64 (e
.logc
, r
), kd
, d
->ln2
);
101 /* y = r2*(A0 + r*A1 + r2*(A2 + r*A3 + r2*A4)) + hi. */
102 r2
= vmulq_f64 (r
, r
);
103 y
= vfmaq_f64 (A (2), A (3), r
);
104 p
= vfmaq_f64 (A (0), A (1), r
);
105 y
= vfmaq_f64 (y
, A (4), r2
);
106 y
= vfmaq_f64 (p
, y
, r2
);
108 if (__glibc_unlikely (v_any_u32h (cmp
)))
109 return special_case (x
, y
, hi
, r2
, cmp
);
110 return vfmaq_f64 (hi
, y
, r2
);