| 1 | /* |
|---|
| 2 | * Single-precision vector log function. |
|---|
| 3 | * |
|---|
| 4 | * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
|---|
| 5 | * See https://llvm.org/LICENSE.txt for license information. |
|---|
| 6 | * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
|---|
| 7 | */ |
|---|
| 8 | |
|---|
| 9 | #include "mathlib.h" |
|---|
| 10 | #include "v_math.h" |
|---|
| 11 | #if V_SUPPORTED |
|---|
| 12 | |
|---|
| 13 | static const float Poly[] = { |
|---|
| 14 | /* 3.34 ulp error */ |
|---|
| 15 | -0x1.3e737cp-3f, 0x1.5a9aa2p-3f, -0x1.4f9934p-3f, 0x1.961348p-3f, |
|---|
| 16 | -0x1.00187cp-2f, 0x1.555d7cp-2f, -0x1.ffffc8p-2f, |
|---|
| 17 | }; |
|---|
| 18 | #define P7 v_f32 (Poly[0]) |
|---|
| 19 | #define P6 v_f32 (Poly[1]) |
|---|
| 20 | #define P5 v_f32 (Poly[2]) |
|---|
| 21 | #define P4 v_f32 (Poly[3]) |
|---|
| 22 | #define P3 v_f32 (Poly[4]) |
|---|
| 23 | #define P2 v_f32 (Poly[5]) |
|---|
| 24 | #define P1 v_f32 (Poly[6]) |
|---|
| 25 | |
|---|
| 26 | #define Ln2 v_f32 (0x1.62e43p-1f) /* 0x3f317218 */ |
|---|
| 27 | #define Min v_u32 (0x00800000) |
|---|
| 28 | #define Max v_u32 (0x7f800000) |
|---|
| 29 | #define Mask v_u32 (0x007fffff) |
|---|
| 30 | #define Off v_u32 (0x3f2aaaab) /* 0.666667 */ |
|---|
| 31 | |
|---|
| 32 | VPCS_ATTR |
|---|
| 33 | __attribute__ ((noinline)) static v_f32_t |
|---|
| 34 | specialcase (v_f32_t x, v_f32_t y, v_u32_t cmp) |
|---|
| 35 | { |
|---|
| 36 | /* Fall back to scalar code. */ |
|---|
| 37 | return v_call_f32 (f: logf, x, y, p: cmp); |
|---|
| 38 | } |
|---|
| 39 | |
|---|
| 40 | VPCS_ATTR |
|---|
| 41 | v_f32_t |
|---|
| 42 | V_NAME(logf) (v_f32_t x) |
|---|
| 43 | { |
|---|
| 44 | v_f32_t n, p, q, r, r2, y; |
|---|
| 45 | v_u32_t u, cmp; |
|---|
| 46 | |
|---|
| 47 | u = v_as_u32_f32 (x); |
|---|
| 48 | cmp = v_cond_u32 (x: u - Min >= Max - Min); |
|---|
| 49 | |
|---|
| 50 | /* x = 2^n * (1+r), where 2/3 < 1+r < 4/3 */ |
|---|
| 51 | u -= Off; |
|---|
| 52 | n = v_to_f32_s32 (x: v_as_s32_u32 (x: u) >> 23); /* signextend */ |
|---|
| 53 | u &= Mask; |
|---|
| 54 | u += Off; |
|---|
| 55 | r = v_as_f32_u32 (x: u) - v_f32 (x: 1.0f); |
|---|
| 56 | |
|---|
| 57 | /* y = log(1+r) + n*ln2. */ |
|---|
| 58 | r2 = r * r; |
|---|
| 59 | /* n*ln2 + r + r2*(P1 + r*P2 + r2*(P3 + r*P4 + r2*(P5 + r*P6 + r2*P7))). */ |
|---|
| 60 | p = v_fma_f32 (P6, y: r, P5); |
|---|
| 61 | q = v_fma_f32 (P4, y: r, P3); |
|---|
| 62 | y = v_fma_f32 (P2, y: r, P1); |
|---|
| 63 | p = v_fma_f32 (P7, y: r2, z: p); |
|---|
| 64 | q = v_fma_f32 (x: p, y: r2, z: q); |
|---|
| 65 | y = v_fma_f32 (x: q, y: r2, z: y); |
|---|
| 66 | p = v_fma_f32 (Ln2, y: n, z: r); |
|---|
| 67 | y = v_fma_f32 (x: y, y: r2, z: p); |
|---|
| 68 | |
|---|
| 69 | if (unlikely (v_any_u32 (cmp))) |
|---|
| 70 | return specialcase (x, y, cmp); |
|---|
| 71 | return y; |
|---|
| 72 | } |
|---|
| 73 | VPCS_ALIAS |
|---|
| 74 | #endif |
|---|
| 75 | |
|---|
