xref: /freebsd/contrib/arm-optimized-routines/math/log2f.c (revision 072a4ba82a01476eaee33781ccd241033eefcf0b)
131914882SAlex Richardson /*
231914882SAlex Richardson  * Single-precision log2 function.
331914882SAlex Richardson  *
431914882SAlex Richardson  * Copyright (c) 2017-2018, Arm Limited.
5*072a4ba8SAndrew Turner  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
631914882SAlex Richardson  */
731914882SAlex Richardson 
831914882SAlex Richardson #include <math.h>
931914882SAlex Richardson #include <stdint.h>
1031914882SAlex Richardson #include "math_config.h"
1131914882SAlex Richardson 
1231914882SAlex Richardson /*
1331914882SAlex Richardson LOG2F_TABLE_BITS = 4
1431914882SAlex Richardson LOG2F_POLY_ORDER = 4
1531914882SAlex Richardson 
1631914882SAlex Richardson ULP error: 0.752 (nearest rounding.)
1731914882SAlex Richardson Relative error: 1.9 * 2^-26 (before rounding.)
1831914882SAlex Richardson */
1931914882SAlex Richardson 
2031914882SAlex Richardson #define N (1 << LOG2F_TABLE_BITS)
2131914882SAlex Richardson #define T __log2f_data.tab
2231914882SAlex Richardson #define A __log2f_data.poly
2331914882SAlex Richardson #define OFF 0x3f330000
2431914882SAlex Richardson 
2531914882SAlex Richardson float
2631914882SAlex Richardson log2f (float x)
2731914882SAlex Richardson {
2831914882SAlex Richardson   /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
2931914882SAlex Richardson   double_t z, r, r2, p, y, y0, invc, logc;
3031914882SAlex Richardson   uint32_t ix, iz, top, tmp;
3131914882SAlex Richardson   int k, i;
3231914882SAlex Richardson 
3331914882SAlex Richardson   ix = asuint (x);
3431914882SAlex Richardson #if WANT_ROUNDING
3531914882SAlex Richardson   /* Fix sign of zero with downward rounding when x==1.  */
3631914882SAlex Richardson   if (unlikely (ix == 0x3f800000))
3731914882SAlex Richardson     return 0;
3831914882SAlex Richardson #endif
3931914882SAlex Richardson   if (unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000))
4031914882SAlex Richardson     {
4131914882SAlex Richardson       /* x < 0x1p-126 or inf or nan.  */
4231914882SAlex Richardson       if (ix * 2 == 0)
4331914882SAlex Richardson 	return __math_divzerof (1);
4431914882SAlex Richardson       if (ix == 0x7f800000) /* log2(inf) == inf.  */
4531914882SAlex Richardson 	return x;
4631914882SAlex Richardson       if ((ix & 0x80000000) || ix * 2 >= 0xff000000)
4731914882SAlex Richardson 	return __math_invalidf (x);
4831914882SAlex Richardson       /* x is subnormal, normalize it.  */
4931914882SAlex Richardson       ix = asuint (x * 0x1p23f);
5031914882SAlex Richardson       ix -= 23 << 23;
5131914882SAlex Richardson     }
5231914882SAlex Richardson 
5331914882SAlex Richardson   /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
5431914882SAlex Richardson      The range is split into N subintervals.
5531914882SAlex Richardson      The ith subinterval contains z and c is near its center.  */
5631914882SAlex Richardson   tmp = ix - OFF;
5731914882SAlex Richardson   i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N;
5831914882SAlex Richardson   top = tmp & 0xff800000;
5931914882SAlex Richardson   iz = ix - top;
6031914882SAlex Richardson   k = (int32_t) tmp >> 23; /* arithmetic shift */
6131914882SAlex Richardson   invc = T[i].invc;
6231914882SAlex Richardson   logc = T[i].logc;
6331914882SAlex Richardson   z = (double_t) asfloat (iz);
6431914882SAlex Richardson 
6531914882SAlex Richardson   /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */
6631914882SAlex Richardson   r = z * invc - 1;
6731914882SAlex Richardson   y0 = logc + (double_t) k;
6831914882SAlex Richardson 
6931914882SAlex Richardson   /* Pipelined polynomial evaluation to approximate log1p(r)/ln2.  */
7031914882SAlex Richardson   r2 = r * r;
7131914882SAlex Richardson   y = A[1] * r + A[2];
7231914882SAlex Richardson   y = A[0] * r2 + y;
7331914882SAlex Richardson   p = A[3] * r + y0;
7431914882SAlex Richardson   y = y * r2 + p;
7531914882SAlex Richardson   return eval_as_float (y);
7631914882SAlex Richardson }
7731914882SAlex Richardson #if USE_GLIBC_ABI
7831914882SAlex Richardson strong_alias (log2f, __log2f_finite)
7931914882SAlex Richardson hidden_alias (log2f, __ieee754_log2f)
8031914882SAlex Richardson #endif
81