xref: /freebsd/contrib/arm-optimized-routines/math/aarch64/sve/log2f.c (revision f3087bef11543b42e0d69b708f367097a4118d24)
1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner  * Single-precision vector/SVE log2 function.
3*f3087befSAndrew Turner  *
4*f3087befSAndrew Turner  * Copyright (c) 2022-2024, Arm Limited.
5*f3087befSAndrew Turner  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*f3087befSAndrew Turner  */
7*f3087befSAndrew Turner 
8*f3087befSAndrew Turner #include "sv_math.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner 
12*f3087befSAndrew Turner static const struct data
13*f3087befSAndrew Turner {
14*f3087befSAndrew Turner   float poly_02468[5];
15*f3087befSAndrew Turner   float poly_1357[4];
16*f3087befSAndrew Turner   uint32_t off, lower;
17*f3087befSAndrew Turner } data = {
18*f3087befSAndrew Turner   .poly_1357 = {
19*f3087befSAndrew Turner     /* Coefficients copied from the AdvSIMD routine, then rearranged so that coeffs
20*f3087befSAndrew Turner        1, 3, 5 and 7 can be loaded as a single quad-word, hence used with _lane
21*f3087befSAndrew Turner        variant of MLA intrinsic.  */
22*f3087befSAndrew Turner     -0x1.715458p-1f, -0x1.7171a4p-2f, -0x1.e5143ep-3f, -0x1.c675bp-3f
23*f3087befSAndrew Turner   },
24*f3087befSAndrew Turner   .poly_02468 = { 0x1.715476p0f, 0x1.ec701cp-2f, 0x1.27a0b8p-2f,
25*f3087befSAndrew Turner 		  0x1.9d8ecap-3f, 0x1.9e495p-3f },
26*f3087befSAndrew Turner   .off = 0x3f2aaaab,
27*f3087befSAndrew Turner   /* Lower bound is the smallest positive normal float 0x00800000. For
28*f3087befSAndrew Turner      optimised register use subnormals are detected after offset has been
29*f3087befSAndrew Turner      subtracted, so lower bound is 0x0080000 - offset (which wraps around).  */
30*f3087befSAndrew Turner   .lower = 0x00800000 - 0x3f2aaaab
31*f3087befSAndrew Turner };
32*f3087befSAndrew Turner 
33*f3087befSAndrew Turner #define Thresh (0x7f000000) /* asuint32(inf) - 0x00800000.  */
34*f3087befSAndrew Turner #define MantissaMask (0x007fffff)
35*f3087befSAndrew Turner 
36*f3087befSAndrew Turner static svfloat32_t NOINLINE
special_case(svuint32_t u_off,svfloat32_t p,svfloat32_t r2,svfloat32_t y,svbool_t cmp)37*f3087befSAndrew Turner special_case (svuint32_t u_off, svfloat32_t p, svfloat32_t r2, svfloat32_t y,
38*f3087befSAndrew Turner 	      svbool_t cmp)
39*f3087befSAndrew Turner {
40*f3087befSAndrew Turner   return sv_call_f32 (
41*f3087befSAndrew Turner       log2f, svreinterpret_f32 (svadd_x (svptrue_b32 (), u_off, data.off)),
42*f3087befSAndrew Turner       svmla_x (svptrue_b32 (), p, r2, y), cmp);
43*f3087befSAndrew Turner }
44*f3087befSAndrew Turner 
45*f3087befSAndrew Turner /* Optimised implementation of SVE log2f, using the same algorithm
46*f3087befSAndrew Turner    and polynomial as AdvSIMD log2f.
47*f3087befSAndrew Turner    Maximum error is 2.48 ULPs:
48*f3087befSAndrew Turner    SV_NAME_F1 (log2)(0x1.558174p+0) got 0x1.a9be84p-2
49*f3087befSAndrew Turner 				   want 0x1.a9be8p-2.  */
SV_NAME_F1(log2)50*f3087befSAndrew Turner svfloat32_t SV_NAME_F1 (log2) (svfloat32_t x, const svbool_t pg)
51*f3087befSAndrew Turner {
52*f3087befSAndrew Turner   const struct data *d = ptr_barrier (&data);
53*f3087befSAndrew Turner 
54*f3087befSAndrew Turner   svuint32_t u_off = svreinterpret_u32 (x);
55*f3087befSAndrew Turner 
56*f3087befSAndrew Turner   u_off = svsub_x (pg, u_off, d->off);
57*f3087befSAndrew Turner   svbool_t special = svcmpge (pg, svsub_x (pg, u_off, d->lower), Thresh);
58*f3087befSAndrew Turner 
59*f3087befSAndrew Turner   /* x = 2^n * (1+r), where 2/3 < 1+r < 4/3.  */
60*f3087befSAndrew Turner   svfloat32_t n = svcvt_f32_x (
61*f3087befSAndrew Turner       pg, svasr_x (pg, svreinterpret_s32 (u_off), 23)); /* Sign-extend.  */
62*f3087befSAndrew Turner   svuint32_t u = svand_x (pg, u_off, MantissaMask);
63*f3087befSAndrew Turner   u = svadd_x (pg, u, d->off);
64*f3087befSAndrew Turner   svfloat32_t r = svsub_x (pg, svreinterpret_f32 (u), 1.0f);
65*f3087befSAndrew Turner 
66*f3087befSAndrew Turner   /* y = log2(1+r) + n.  */
67*f3087befSAndrew Turner   svfloat32_t r2 = svmul_x (svptrue_b32 (), r, r);
68*f3087befSAndrew Turner 
69*f3087befSAndrew Turner   /* Evaluate polynomial using pairwise Horner scheme.  */
70*f3087befSAndrew Turner   svfloat32_t p_1357 = svld1rq (svptrue_b32 (), &d->poly_1357[0]);
71*f3087befSAndrew Turner   svfloat32_t q_01 = svmla_lane (sv_f32 (d->poly_02468[0]), r, p_1357, 0);
72*f3087befSAndrew Turner   svfloat32_t q_23 = svmla_lane (sv_f32 (d->poly_02468[1]), r, p_1357, 1);
73*f3087befSAndrew Turner   svfloat32_t q_45 = svmla_lane (sv_f32 (d->poly_02468[2]), r, p_1357, 2);
74*f3087befSAndrew Turner   svfloat32_t q_67 = svmla_lane (sv_f32 (d->poly_02468[3]), r, p_1357, 3);
75*f3087befSAndrew Turner   svfloat32_t y = svmla_x (pg, q_67, r2, sv_f32 (d->poly_02468[4]));
76*f3087befSAndrew Turner   y = svmla_x (pg, q_45, r2, y);
77*f3087befSAndrew Turner   y = svmla_x (pg, q_23, r2, y);
78*f3087befSAndrew Turner   y = svmla_x (pg, q_01, r2, y);
79*f3087befSAndrew Turner 
80*f3087befSAndrew Turner   if (unlikely (svptest_any (pg, special)))
81*f3087befSAndrew Turner     return special_case (u_off, n, r, y, special);
82*f3087befSAndrew Turner   return svmla_x (svptrue_b32 (), n, r, y);
83*f3087befSAndrew Turner }
84*f3087befSAndrew Turner 
85*f3087befSAndrew Turner TEST_SIG (SV, F, 1, log2, 0.01, 11.1)
86*f3087befSAndrew Turner TEST_ULP (SV_NAME_F1 (log2), 1.99)
87*f3087befSAndrew Turner TEST_DISABLE_FENV (SV_NAME_F1 (log2))
88*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), -0.0, -0x1p126, 4000)
89*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), 0.0, 0x1p-126, 4000)
90*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), 0x1p-126, 0x1p-23, 50000)
91*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), 0x1p-23, 1.0, 50000)
92*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), 1.0, 100, 50000)
93*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_F1 (log2), 100, inf, 50000)
94*f3087befSAndrew Turner CLOSE_SVE_ATTR
95