[PATCH 12/17] math: Use log10f from CORE-MATH

[PATCH 12/17] math: Use log10f from CORE-MATH

[caiyinyu]

[-- Attachment #1: Type: text/plain, Size: 2896 bytes --]


HI,

The commit 9247f53219 "math: Use log10f from CORE-MATH" caused the 
following tests to fail for loongarch and riscv when input values are 0/-0.
math/test-float-log10
math/test-float32-log10

The tests passed on x86_64, aarch64, and other architectures because 
they had access to the log10f function through the 
math/w_log10f_compat.c wrapper code,

which was available in glibc before 2.27 for x86_64 and probably earlier 
for other architectures.

However, riscv and loongarch architectures entered the community later, 
with riscv in 2.27 and loongarch in 2.36.

Is this a bug? Or do riscv and loongarch need to implement their own 
related functions?

riscv debug site:
(gdb) disassemble /m __math_invalidf
Dump of assembler code for function __math_invalidf:
79    {
    0x0000003ff7fac6da <+0>:    fmv.s    fa4,fa0
    0x0000003ff7fac6ec <+18>:    addi    sp,sp,-16
    0x0000003ff7fac6ee <+20>:    sd    ra,8(sp)

80      float y = (x - x) / (x - x);
    0x0000003ff7fac6de <+4>:    fsub.s    fa5,fa0,fa0
    0x0000003ff7fac6e2 <+8>:    fdiv.s fa0,fa5,fa5              
<------   here 0/0 , invalid operation set

81      return isnan (x) ? y : with_errnof (y, EDOM);
=> 0x0000003ff7fac6e6 <+12>:    feq.s    a5,fa4,fa4
    0x0000003ff7fac6ea <+16>:    beqz    a5,0x3ff7fac6fe 
<__math_invalidf+36>
    0x0000003ff7fac6f0 <+22>:    li    a0,33
    0x0000003ff7fac6f4 <+26>:    jal    ra,0x3ff7fac628 <with_errnof>

82    }
    0x0000003ff7fac6f8 <+30>:    ld    ra,8(sp)
    0x0000003ff7fac6fa <+32>:    addi    sp,sp,16
    0x0000003ff7fac6fc <+34>:    ret
    0x0000003ff7fac6fe <+36>:    ret

(gdb) bt
#0  __math_invalidf (x=<optimized out>) at 
../sysdeps/ieee754/flt-32/math_errf.c:81
#1  0x0000003ff7fa6b14 in as_special (x=<optimized out>) at 
../sysdeps/ieee754/flt-32/e_log10f.c:45
#2  0x0000003ff7fa6c40 in __ieee754_log10f (x=<optimized out>) at 
../sysdeps/ieee754/flt-32/e_log10f.c:121
#3  0x0000003ff7fa2c6e in __log10f (x=<optimized out>) at 
./w_log10_template.c:41
#4  0x00000000000123c0 in log10_test () at 
/home/caiyinyu/glibc/build/math/libm-test-log10.c:908
#5  0x000000000001265e in do_test () at 
/home/caiyinyu/glibc/build/math/libm-test-log10.c:914
#6  0x0000000000012670 in main (argc=<optimized out>, argv=<optimized 
out>) at ./libm-test-driver.c:1135
#7  0x0000003ff7e68322 in __libc_start_call_main (main=<optimized out>, 
argc=<optimized out>, argv=<optimized out>) at 
../sysdeps/nptl/libc_start_call_main.h:58
#8  0x0000003ff7e683aa in __libc_start_main_impl (main=<optimized out>, 
argc=<optimized out>, argv=<optimized out>, init=<optimized out>, 
fini=<optimized out>, rtld_fini=<optimized out>, stack_end=<optimized 
out>) at libc-start.c:360
#9  0x00000000000122c0 in _start () at ../sysdeps/riscv/start.S:67

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[Xi Ruoyao]

On Thu, 2024-11-07 at 17:21 +0800, caiyinyu wrote:
>  HI, 
>  
>  The commit 9247f53219 "math: Use log10f from CORE-MATH" caused the
> following tests to fail for loongarch and riscv when input values are
> 0/-0. 
>  math/test-float-log10 
>  math/test-float32-log10 
>  
>  The tests passed on x86_64, aarch64, and other architectures because
> they had access to the log10f function through the
> math/w_log10f_compat.c wrapper code, 
>  
>  which was available in glibc before 2.27 for x86_64 and probably
> earlier for other architectures. 
>  
>  However, riscv and loongarch architectures entered the community
> later, with riscv in 2.27 and loongarch in 2.36. 
>  
>  Is this a bug? Or do riscv and loongarch need to implement their own
> related functions? 

This seems a bug.  It looks to me we need

diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c
index 058ce31f09..61fb84c7aa 100644
--- a/sysdeps/ieee754/flt-32/e_log10f.c
+++ b/sysdeps/ieee754/flt-32/e_log10f.c
@@ -38,7 +38,7 @@ as_special (float x)
   uint32_t ax = ux << 1;
   if (ax == 0u)
     { /* -0.0 */
-      __math_divzerof (1);
+      return __math_divzerof (1);
     }
   if (ax > 0xff000000u)
     return x + x; /* nan */

(not tested)

-- 
Xi Ruoyao <xry111@xry111.site>
School of Aerospace Science and Technology, Xidian University

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[Xi Ruoyao]

On Thu, 2024-11-07 at 17:33 +0800, Xi Ruoyao wrote:
> On Thu, 2024-11-07 at 17:21 +0800, caiyinyu wrote:
> >  HI, 
> >  
> >  The commit 9247f53219 "math: Use log10f from CORE-MATH" caused the
> > following tests to fail for loongarch and riscv when input values are
> > 0/-0. 
> >  math/test-float-log10 
> >  math/test-float32-log10 
> >  
> >  The tests passed on x86_64, aarch64, and other architectures because
> > they had access to the log10f function through the
> > math/w_log10f_compat.c wrapper code, 
> >  
> >  which was available in glibc before 2.27 for x86_64 and probably
> > earlier for other architectures. 
> >  
> >  However, riscv and loongarch architectures entered the community
> > later, with riscv in 2.27 and loongarch in 2.36. 
> >  
> >  Is this a bug? Or do riscv and loongarch need to implement their own
> > related functions? 
> 
> This seems a bug.  It looks to me we need
> 
> diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c
> index 058ce31f09..61fb84c7aa 100644
> --- a/sysdeps/ieee754/flt-32/e_log10f.c
> +++ b/sysdeps/ieee754/flt-32/e_log10f.c
> @@ -38,7 +38,7 @@ as_special (float x)
>    uint32_t ax = ux << 1;
>    if (ax == 0u)
>      { /* -0.0 */
> -      __math_divzerof (1);
> +      return __math_divzerof (1);
>      }
>    if (ax > 0xff000000u)
>      return x + x; /* nan */
> 
> (not tested)

In the original CORE-MATH there is a return anyway:
https://gitlab.inria.fr/core-math/core-math/-/blob/61d7befe655ec7dea6461cd5658c607c1db48462/src/binary32/log10/log10f.c#L49

-- 
Xi Ruoyao <xry111@xry111.site>
School of Aerospace Science and Technology, Xidian University

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[Adhemerval Zanella Netto]

On 07/11/24 06:39, Xi Ruoyao wrote:
> On Thu, 2024-11-07 at 17:33 +0800, Xi Ruoyao wrote:
>> On Thu, 2024-11-07 at 17:21 +0800, caiyinyu wrote:
>>>  HI, 
>>>  
>>>  The commit 9247f53219 "math: Use log10f from CORE-MATH" caused the
>>> following tests to fail for loongarch and riscv when input values are
>>> 0/-0. 
>>>  math/test-float-log10 
>>>  math/test-float32-log10 
>>>  
>>>  The tests passed on x86_64, aarch64, and other architectures because
>>> they had access to the log10f function through the
>>> math/w_log10f_compat.c wrapper code, 
>>>  
>>>  which was available in glibc before 2.27 for x86_64 and probably
>>> earlier for other architectures. 
>>>  
>>>  However, riscv and loongarch architectures entered the community
>>> later, with riscv in 2.27 and loongarch in 2.36. 
>>>  
>>>  Is this a bug? Or do riscv and loongarch need to implement their own
>>> related functions? 
>>
>> This seems a bug.  It looks to me we need
>>
>> diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c
>> index 058ce31f09..61fb84c7aa 100644
>> --- a/sysdeps/ieee754/flt-32/e_log10f.c
>> +++ b/sysdeps/ieee754/flt-32/e_log10f.c
>> @@ -38,7 +38,7 @@ as_special (float x)
>>    uint32_t ax = ux << 1;
>>    if (ax == 0u)
>>      { /* -0.0 */
>> -      __math_divzerof (1);
>> +      return __math_divzerof (1);
>>      }
>>    if (ax > 0xff000000u)
>>      return x + x; /* nan */
>>
>> (not tested)
> 
> In the original CORE-MATH there is a return anyway:
> https://gitlab.inria.fr/core-math/core-math/-/blob/61d7befe655ec7dea6461cd5658c607c1db48462/src/binary32/log10/log10f.c#L49
> 

It was my mistake, thanks for catching this. I will fix this and I try to
extend the ABI range I usually tests (I tested on x86_64, i686, aarch64,
armhf, and powerpc64le).

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[Adhemerval Zanella Netto]

On 07/11/24 08:57, Wilco Dijkstra wrote:
> Hi Adhemerval,
> 
> I guess this shows we also need to remove the wrappers for these functions so
> that they are equal on all targets.
> 
> Cheers,
> Wilco

Yes, I avoided to this now because it would require newer symbols and lot of
internal machinery. But this is on my backlog.

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[Wilco Dijkstra]

Hi Adhemerval,

I guess this shows we also need to remove the wrappers for these functions so
that they are equal on all targets.

Cheers,
Wilco

Re: [PATCH 12/17] math: Use log10f from CORE-MATH

[DJ Delorie]

LGTM
Reviewed-by: DJ Delorie <dj@redhat.com>

Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
> diff --git a/sysdeps/i386/fpu/e_log10f.S b/sysdeps/i386/fpu/e_log10f.S
> deleted file mode 100644

Ok.

> diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c

> -/* e_log10f.c -- float version of e_log10.c.
> - */
> -
> -/*
> - * ====================================================
> - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
> - *
> - * Developed at SunPro, a Sun Microsystems, Inc. business.
> - * Permission to use, copy, modify, and distribute this
> - * software is freely granted, provided that this notice
> - * is preserved.
> - * ====================================================
> - */
> +/* Correctly-rounded radix-10 logarithm function for binary32 value.
> +
> +Copyright (c) 2022-2023 Alexei Sibidanov.
> +
> +This file is part of the CORE-MATH project
> +project (file src/binary32/expm1/expm1f.c, revision bc385c2).
> +
> +Permission is hereby granted, free of charge, to any person obtaining a copy
> +of this software and associated documentation files (the "Software"), to deal
> +in the Software without restriction, including without limitation the rights
> +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> +copies of the Software, and to permit persons to whom the Software is
> +furnished to do so, subject to the following conditions:
> +
> +The above copyright notice and this permission notice shall be included in all
> +copies or substantial portions of the Software.
> +
> +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> +SOFTWARE.
> +*/

Ok.

>  #include <math.h>
> -#include <math_private.h>
> -#include <fix-int-fp-convert-zero.h>
> +#include <stdint.h>
>  #include <libm-alias-finite.h>
> +#include "math_config.h"

Ok.

> -static const float
> -two25      =  3.3554432000e+07, /* 0x4c000000 */
> -ivln10     =  4.3429449201e-01, /* 0x3ede5bd9 */
> -log10_2hi  =  3.0102920532e-01, /* 0x3e9a2080 */
> -log10_2lo  =  7.9034151668e-07; /* 0x355427db */

> +static __attribute__ ((noinline)) float
> +as_special (float x)
> +{
> +  uint32_t ux = asuint (x);
> +  if (ux == 0x7f800000u)
> +    return x; /* +inf */
> +  uint32_t ax = ux << 1;
> +  if (ax == 0u)
> +    { /* -0.0 */
> +      __math_divzerof (1);
> +    }
> +  if (ax > 0xff000000u)
> +    return x + x; /* nan */
> +  return __math_invalidf (x);
> +}

Ok.

>  float
> -__ieee754_log10f(float x)
> +__ieee754_log10f (float x)

Ok.

>  {
> -	float y,z;
> -	int32_t i,k,hx;
> -
> -	GET_FLOAT_WORD(hx,x);
> -
> -	k=0;
> -	if (hx < 0x00800000) {			/* x < 2**-126  */
> -	    if (__builtin_expect((hx&0x7fffffff)==0, 0))
> -	      return -two25/fabsf (x);	/* log(+-0)=-inf  */
> -	    if (__builtin_expect(hx<0, 0))
> -		return (x-x)/(x-x);	/* log(-#) = NaN */
> -	    k -= 25; x *= two25; /* subnormal number, scale up x */
> -	    GET_FLOAT_WORD(hx,x);

> +  static const double tr[] =
> +    {
> +      0x1p+0,         0x1.f81f82p-1,  0x1.f07c1fp-1,  0x1.e9131acp-1,
> +      0x1.e1e1e1ep-1, 0x1.dae6077p-1, 0x1.d41d41dp-1, 0x1.cd85689p-1,
> +      0x1.c71c71cp-1, 0x1.c0e0704p-1, 0x1.bacf915p-1, 0x1.b4e81b5p-1,
> +      0x1.af286bdp-1, 0x1.a98ef6p-1,  0x1.a41a41ap-1, 0x1.9ec8e95p-1,
> +      0x1.999999ap-1, 0x1.948b0fdp-1, 0x1.8f9c19p-1,  0x1.8acb90fp-1,
> +      0x1.8618618p-1, 0x1.8181818p-1, 0x1.7d05f41p-1, 0x1.78a4c81p-1,
> +      0x1.745d174p-1, 0x1.702e05cp-1, 0x1.6c16c17p-1, 0x1.6816817p-1,
> +      0x1.642c859p-1, 0x1.605816p-1,  0x1.5c9882cp-1, 0x1.58ed231p-1,
> +      0x1.5555555p-1, 0x1.51d07ebp-1, 0x1.4e5e0a7p-1, 0x1.4afd6ap-1,
> +      0x1.47ae148p-1, 0x1.446f865p-1, 0x1.4141414p-1, 0x1.3e22cbdp-1,
> +      0x1.3b13b14p-1, 0x1.3813814p-1, 0x1.3521cfbp-1, 0x1.323e34ap-1,
> +      0x1.2f684bep-1, 0x1.2c9fb4ep-1, 0x1.29e412ap-1, 0x1.27350b9p-1,
> +      0x1.2492492p-1, 0x1.21fb781p-1, 0x1.1f7047ep-1, 0x1.1cf06aep-1,
> +      0x1.1a7b961p-1, 0x1.1811812p-1, 0x1.15b1e5fp-1, 0x1.135c811p-1,
> +      0x1.1111111p-1, 0x1.0ecf56cp-1, 0x1.0c9715p-1,  0x1.0a6810ap-1,
> +      0x1.0842108p-1, 0x1.0624dd3p-1, 0x1.041041p-1,  0x1.0204081p-1,
> +      0.5
> +    };
> +  static const double tl[] =
> +    {
> +      -0x1.d45fd6237ebe3p-47, 0x1.b947689311b6ep-8, 0x1.b5e909c96d7d5p-7,
> +      0x1.45f4f59ed2165p-6,   0x1.af5f92cbd8f1ep-6, 0x1.0ba01a606de8cp-5,
> +      0x1.3ed119b9a2b7bp-5,   0x1.714834298eec2p-5, 0x1.a30a9d98357fbp-5,
> +      0x1.d41d512670813p-5,   0x1.02428c0f65519p-4, 0x1.1a23444eecc3ep-4,
> +      0x1.31b30543f4cb4p-4,   0x1.48f3ed39bfd04p-4, 0x1.5fe8049a0e423p-4,
> +      0x1.769140a6aa008p-4,   0x1.8cf1836c98cb3p-4, 0x1.a30a9d55541a1p-4,
> +      0x1.b8de4d1ee823ep-4,   0x1.ce6e4202ca2e6p-4, 0x1.e3bc1accace07p-4,
> +      0x1.f8c9683b5abd4p-4,   0x1.06cbd68ca9a6ep-3, 0x1.11142f19df73p-3,
> +      0x1.1b3e71fa7a97fp-3,   0x1.254b4d37a46e3p-3, 0x1.2f3b6912cbf07p-3,
> +      0x1.390f683115886p-3,   0x1.42c7e7fffc5a8p-3, 0x1.4c65808c78d3cp-3,
> +      0x1.55e8c50751c55p-3,   0x1.5f52445dec3d8p-3, 0x1.68a288c3f12p-3,
> +      0x1.71da17bdf0d19p-3,   0x1.7af973608afd9p-3, 0x1.84011952a2579p-3,
> +      0x1.8cf1837a7ea6p-3,    0x1.95cb2891e43d6p-3, 0x1.9e8e7b0f869ep-3,
> +      0x1.a73beaa5db18dp-3,   0x1.afd3e394558d3p-3, 0x1.b856cf060d9f1p-3,
> +      0x1.c0c5134de1ffcp-3,   0x1.c91f1371bc99fp-3, 0x1.d1652ffcd3f53p-3,
> +      0x1.d997c6f635e75p-3,   0x1.e1b733ab90f3bp-3, 0x1.e9c3ceadac856p-3,
> +      0x1.f1bdeec43a305p-3,   0x1.f9a5e7a5fa3fep-3, 0x1.00be05ac02f2bp-2,
> +      0x1.04a054d81a2d4p-2,   0x1.087a0835957fbp-2, 0x1.0c4b457099517p-2,
> +      0x1.101431aa1fe51p-2,   0x1.13d4f08b98dd8p-2, 0x1.178da53edb892p-2,
> +      0x1.1b3e71e9f9d58p-2,   0x1.1ee777defdeedp-2, 0x1.2288d7b48e23bp-2,
> +      0x1.2622b0f52e49fp-2,   0x1.29b522a4c6314p-2, 0x1.2d404b0e30f8p-2,
> +      0x1.30c4478f3fbe5p-2,   0x1.34413509f7915p-2
> +    };
> +  static const union
> +  {
> +    float f;
> +    uint32_t u;
> +  } st[] =
> +  {
> +    { 0x1p+0 },        { 0x1.4p+3 },      { 0x1.9p+6 },       { 0x1.f4p+9 },
> +    { 0x1.388p+13 },   { 0x1.86ap+16 },   { 0x1.e848p+19 },   { 0x1.312dp+23 },
> +    { 0x1.7d784p+26 }, { 0x1.dcd65p+29 }, { 0x1.2a05f2p+33 }, { 0 },
> +    { 0 },             { 0 },             { 0 },              { 0 }
> +  };
> +  static const double b[] =
> +    { 
> +      0x1.bcb7b15c5a2f8p-2, -0x1.bcbb1dbb88ebap-3, 0x1.2871c39d521c6p-3
> +    };
> +  static const double c[] =
> +    {
> +      0x1.bcb7b1526e50ep-2,  -0x1.bcb7b1526e53dp-3, 0x1.287a7636f3fa2p-3,
> +      -0x1.bcb7b146a14b3p-4, 0x1.63c627d5219cbp-4,  -0x1.2880736c8762dp-4,
> +      0x1.fc1ecf913961ap-5
> +    };

Ok.

> +  uint32_t ux = asuint (x);
> +  if (__glibc_unlikely (ux < (1 << 23) || ux >= 0x7f800000u))
> +    {
> +      if (ux == 0 || ux >= 0x7f800000u)
> +	return as_special (x);
> +      /* subnormal */
> +      int n = __builtin_clz (ux) - 8;
> +      ux <<= n;
> +      ux -= n << 23;
> +    }

Ok.

> +  unsigned m = ux & ((1 << 23) - 1), j = (m + (1 << (23 - 7))) >> (23 - 6);
> +  double ix = tr[j], l = tl[j];
> +  int e = ((int) ux >> 23) - 127;
> +  unsigned je = e + 1;
> +  je = (je * 0x4d104d4) >> 28;
> +  if (__glibc_unlikely (ux == st[je].u))
> +    return je;

Ok.

> +  double tz = asdouble (((int64_t) m | ((int64_t) 1023 << 23)) << (52 - 23));
> +  double z = tz * ix - 1, z2 = z * z;
> +  double r
> +      = ((e * 0x1.34413509f79ffp-2 + l) + z * b[0]) + z2 * (b[1] + z * b[2]);
> +  float ub = r, lb = r + 0x1.b008p-34;
> +  if (__glibc_unlikely (ub != lb))
> +    {
> +      double f = z
> +		 * ((c[0] + z * c[1])
> +		    + z2
> +			  * ((c[2] + z * c[3])
> +			     + z2 * (c[4] + z * c[5] + z2 * c[6])));

Odd looking but OK.

> +      f -= 0x1.0cee0ed4ca7e9p-54 * e;
> +      f += l - tl[0];
> +      double el = e * 0x1.34413509f7ap-2;

Ok.

> +      r = el + f;
> +      ub = r;
> +      tz = r;
> +      if (__glibc_unlikely (!((asuint64 (tz) & ((1 << 28) - 1)))))
> +	{
> +	  double dr = (el - r) + f;
> +	  r += dr * 32;
> +	  ub = r;
>  	}
> -	if (__builtin_expect(hx >= 0x7f800000, 0)) return x+x;
> -	k += (hx>>23)-127;
> -	i  = ((uint32_t)k&0x80000000)>>31;
> -	hx = (hx&0x007fffff)|((0x7f-i)<<23);
> -	y  = (float)(k+i);
> -	if (FIX_INT_FP_CONVERT_ZERO && y == 0.0f)
> -	  y = 0.0f;
> -	SET_FLOAT_WORD(x,hx);
> -	z  = y*log10_2lo + ivln10*__ieee754_logf(x);
> -	return  z+y*log10_2hi;
> +    }
> +  return ub;
>  }
>  libm_alias_finite (__ieee754_log10f, __log10f)

Ok.

[PATCH 12/17] math: Use log10f from CORE-MATH

[Adhemerval Zanella]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain; charset=UTF-8, Size: 12012 bytes --]

The CORE-MATH implementation is correctly rounded (for any rounding mode)
and shows better performance compared to the generic log10f.

The code was adapted to glibc style and to use the definition of
math_config.h (to handle errno, overflow, and underflow).

Benchtest on x64_64 (Ryzen 9 5900X, gcc 14.2.1), aarch64 (Neoverse-N1,
gcc 13.3.1), and powerpc (POWER10, gcc 13.2.1):

Latency                      master        patched   improvement
x86_64                      49.9017        33.5143        32.84%
x86_64v2                    50.4878        33.5623        33.52%
x86_64v3                    50.0991        27.6078        44.89%
i686                        140.874        106.086        24.69%
aarch64                     19.2846        11.3573        41.11%
power10                     14.0994        7.7739        44.86%
powerpc                     14.2898        7.92497        44.54%

reciprocal-throughput        master        patched   improvement
x86_64                      17.8336        12.9074        27.62%
x86_64v2                    16.4418        11.3220        31.14%
x86_64v3                    15.6002        10.5158        32.59%
i686                        66.0678        80.2287        -21.43%
aarch64                      9.4906        6.8393        27.94%
power10                      7.5255        5.5084        26.80%
powerpc                      9.5204        6.98055        26.68%

The performance decrease for i686 is mostly due the use of x87 fpu,
when building with '-msse2 -mfpmath=sse':

                             master        patched   improvement
latency                     140.874        77.1137        45.26%
reciprocal-throughput        64.481        56.4397        12.47%

Signed-off-by: Alexei Sibidanov <sibid@uvic.ca>
Signed-off-by: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Signed-off-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
---
 sysdeps/i386/fpu/e_log10f.S       |  66 ----------
 sysdeps/ieee754/flt-32/e_log10f.c | 196 +++++++++++++++++++++++-------
 2 files changed, 152 insertions(+), 110 deletions(-)
 delete mode 100644 sysdeps/i386/fpu/e_log10f.S

diff --git a/sysdeps/i386/fpu/e_log10f.S b/sysdeps/i386/fpu/e_log10f.S
deleted file mode 100644
index 47f82e1993..0000000000
--- a/sysdeps/i386/fpu/e_log10f.S
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include <libm-alias-finite.h>
-
-	.section .rodata.cst8,"aM",@progbits,8
-
-	.p2align 3
-	.type one,@object
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-	/* It is not important that this constant is precise.  It is only
-	   a value which is known to be on the safe side for using the
-	   fyl2xp1 instruction.  */
-	.type limit,@object
-limit:	.double 0.29
-	ASM_SIZE_DIRECTIVE(limit)
-
-
-#ifdef PIC
-# define MO(op) op##@GOTOFF(%edx)
-#else
-# define MO(op) op
-#endif
-
-	.text
-ENTRY(__ieee754_log10f)
-	fldlg2			// log10(2)
-	flds	4(%esp)		// x : log10(2)
-#ifdef PIC
-	LOAD_PIC_REG (dx)
-#endif
-	fxam
-	fnstsw
-	fld	%st		// x : x : log10(2)
-	sahf
-	jc	3f		// in case x is NaN or ±Inf
-4:	fsubl	MO(one)		// x-1 : x : log10(2)
-	fld	%st		// x-1 : x-1 : x : log10(2)
-	fabs			// |x-1| : x-1 : x : log10(2)
-	fcompl	MO(limit)	// x-1 : x : log10(2)
-	fnstsw			// x-1 : x : log10(2)
-	andb	$0x45, %ah
-	jz	2f
-	fxam
-	fnstsw
-	andb	$0x45, %ah
-	cmpb	$0x40, %ah
-	jne	5f
-	fabs			// log10(1) is +0 in all rounding modes.
-5:	fstp	%st(1)		// x-1 : log10(2)
-	fyl2xp1			// log10(x)
-	ret
-
-2:	fstp	%st(0)		// x : log10(2)
-	fyl2x			// log10(x)
-	ret
-
-3:	jp	4b		// in case x is ±Inf
-	fstp	%st(1)
-	fstp	%st(1)
-	ret
-END (__ieee754_log10f)
-libm_alias_finite (__ieee754_log10f, __log10f)
diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c
index 791895e042..8620805634 100644
--- a/sysdeps/ieee754/flt-32/e_log10f.c
+++ b/sysdeps/ieee754/flt-32/e_log10f.c
@@ -1,54 +1,162 @@
-/* e_log10f.c -- float version of e_log10.c.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+/* Correctly-rounded radix-10 logarithm function for binary32 value.
+
+Copyright (c) 2022-2023 Alexei Sibidanov.
+
+This file is part of the CORE-MATH project
+project (file src/binary32/expm1/expm1f.c, revision bc385c2).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
 
 #include <math.h>
-#include <math_private.h>
-#include <fix-int-fp-convert-zero.h>
+#include <stdint.h>
 #include <libm-alias-finite.h>
+#include "math_config.h"
 
-static const float
-two25      =  3.3554432000e+07, /* 0x4c000000 */
-ivln10     =  4.3429449201e-01, /* 0x3ede5bd9 */
-log10_2hi  =  3.0102920532e-01, /* 0x3e9a2080 */
-log10_2lo  =  7.9034151668e-07; /* 0x355427db */
+static __attribute__ ((noinline)) float
+as_special (float x)
+{
+  uint32_t ux = asuint (x);
+  if (ux == 0x7f800000u)
+    return x; /* +inf */
+  uint32_t ax = ux << 1;
+  if (ax == 0u)
+    { /* -0.0 */
+      __math_divzerof (1);
+    }
+  if (ax > 0xff000000u)
+    return x + x; /* nan */
+  return __math_invalidf (x);
+}
 
 float
-__ieee754_log10f(float x)
+__ieee754_log10f (float x)
 {
-	float y,z;
-	int32_t i,k,hx;
-
-	GET_FLOAT_WORD(hx,x);
-
-	k=0;
-	if (hx < 0x00800000) {			/* x < 2**-126  */
-	    if (__builtin_expect((hx&0x7fffffff)==0, 0))
-	      return -two25/fabsf (x);	/* log(+-0)=-inf  */
-	    if (__builtin_expect(hx<0, 0))
-		return (x-x)/(x-x);	/* log(-#) = NaN */
-	    k -= 25; x *= two25; /* subnormal number, scale up x */
-	    GET_FLOAT_WORD(hx,x);
+  static const double tr[] =
+    {
+      0x1p+0,         0x1.f81f82p-1,  0x1.f07c1fp-1,  0x1.e9131acp-1,
+      0x1.e1e1e1ep-1, 0x1.dae6077p-1, 0x1.d41d41dp-1, 0x1.cd85689p-1,
+      0x1.c71c71cp-1, 0x1.c0e0704p-1, 0x1.bacf915p-1, 0x1.b4e81b5p-1,
+      0x1.af286bdp-1, 0x1.a98ef6p-1,  0x1.a41a41ap-1, 0x1.9ec8e95p-1,
+      0x1.999999ap-1, 0x1.948b0fdp-1, 0x1.8f9c19p-1,  0x1.8acb90fp-1,
+      0x1.8618618p-1, 0x1.8181818p-1, 0x1.7d05f41p-1, 0x1.78a4c81p-1,
+      0x1.745d174p-1, 0x1.702e05cp-1, 0x1.6c16c17p-1, 0x1.6816817p-1,
+      0x1.642c859p-1, 0x1.605816p-1,  0x1.5c9882cp-1, 0x1.58ed231p-1,
+      0x1.5555555p-1, 0x1.51d07ebp-1, 0x1.4e5e0a7p-1, 0x1.4afd6ap-1,
+      0x1.47ae148p-1, 0x1.446f865p-1, 0x1.4141414p-1, 0x1.3e22cbdp-1,
+      0x1.3b13b14p-1, 0x1.3813814p-1, 0x1.3521cfbp-1, 0x1.323e34ap-1,
+      0x1.2f684bep-1, 0x1.2c9fb4ep-1, 0x1.29e412ap-1, 0x1.27350b9p-1,
+      0x1.2492492p-1, 0x1.21fb781p-1, 0x1.1f7047ep-1, 0x1.1cf06aep-1,
+      0x1.1a7b961p-1, 0x1.1811812p-1, 0x1.15b1e5fp-1, 0x1.135c811p-1,
+      0x1.1111111p-1, 0x1.0ecf56cp-1, 0x1.0c9715p-1,  0x1.0a6810ap-1,
+      0x1.0842108p-1, 0x1.0624dd3p-1, 0x1.041041p-1,  0x1.0204081p-1,
+      0.5
+    };
+  static const double tl[] =
+    {
+      -0x1.d45fd6237ebe3p-47, 0x1.b947689311b6ep-8, 0x1.b5e909c96d7d5p-7,
+      0x1.45f4f59ed2165p-6,   0x1.af5f92cbd8f1ep-6, 0x1.0ba01a606de8cp-5,
+      0x1.3ed119b9a2b7bp-5,   0x1.714834298eec2p-5, 0x1.a30a9d98357fbp-5,
+      0x1.d41d512670813p-5,   0x1.02428c0f65519p-4, 0x1.1a23444eecc3ep-4,
+      0x1.31b30543f4cb4p-4,   0x1.48f3ed39bfd04p-4, 0x1.5fe8049a0e423p-4,
+      0x1.769140a6aa008p-4,   0x1.8cf1836c98cb3p-4, 0x1.a30a9d55541a1p-4,
+      0x1.b8de4d1ee823ep-4,   0x1.ce6e4202ca2e6p-4, 0x1.e3bc1accace07p-4,
+      0x1.f8c9683b5abd4p-4,   0x1.06cbd68ca9a6ep-3, 0x1.11142f19df73p-3,
+      0x1.1b3e71fa7a97fp-3,   0x1.254b4d37a46e3p-3, 0x1.2f3b6912cbf07p-3,
+      0x1.390f683115886p-3,   0x1.42c7e7fffc5a8p-3, 0x1.4c65808c78d3cp-3,
+      0x1.55e8c50751c55p-3,   0x1.5f52445dec3d8p-3, 0x1.68a288c3f12p-3,
+      0x1.71da17bdf0d19p-3,   0x1.7af973608afd9p-3, 0x1.84011952a2579p-3,
+      0x1.8cf1837a7ea6p-3,    0x1.95cb2891e43d6p-3, 0x1.9e8e7b0f869ep-3,
+      0x1.a73beaa5db18dp-3,   0x1.afd3e394558d3p-3, 0x1.b856cf060d9f1p-3,
+      0x1.c0c5134de1ffcp-3,   0x1.c91f1371bc99fp-3, 0x1.d1652ffcd3f53p-3,
+      0x1.d997c6f635e75p-3,   0x1.e1b733ab90f3bp-3, 0x1.e9c3ceadac856p-3,
+      0x1.f1bdeec43a305p-3,   0x1.f9a5e7a5fa3fep-3, 0x1.00be05ac02f2bp-2,
+      0x1.04a054d81a2d4p-2,   0x1.087a0835957fbp-2, 0x1.0c4b457099517p-2,
+      0x1.101431aa1fe51p-2,   0x1.13d4f08b98dd8p-2, 0x1.178da53edb892p-2,
+      0x1.1b3e71e9f9d58p-2,   0x1.1ee777defdeedp-2, 0x1.2288d7b48e23bp-2,
+      0x1.2622b0f52e49fp-2,   0x1.29b522a4c6314p-2, 0x1.2d404b0e30f8p-2,
+      0x1.30c4478f3fbe5p-2,   0x1.34413509f7915p-2
+    };
+  static const union
+  {
+    float f;
+    uint32_t u;
+  } st[] =
+  {
+    { 0x1p+0 },        { 0x1.4p+3 },      { 0x1.9p+6 },       { 0x1.f4p+9 },
+    { 0x1.388p+13 },   { 0x1.86ap+16 },   { 0x1.e848p+19 },   { 0x1.312dp+23 },
+    { 0x1.7d784p+26 }, { 0x1.dcd65p+29 }, { 0x1.2a05f2p+33 }, { 0 },
+    { 0 },             { 0 },             { 0 },              { 0 }
+  };
+  static const double b[] =
+    { 
+      0x1.bcb7b15c5a2f8p-2, -0x1.bcbb1dbb88ebap-3, 0x1.2871c39d521c6p-3
+    };
+  static const double c[] =
+    {
+      0x1.bcb7b1526e50ep-2,  -0x1.bcb7b1526e53dp-3, 0x1.287a7636f3fa2p-3,
+      -0x1.bcb7b146a14b3p-4, 0x1.63c627d5219cbp-4,  -0x1.2880736c8762dp-4,
+      0x1.fc1ecf913961ap-5
+    };
+  uint32_t ux = asuint (x);
+  if (__glibc_unlikely (ux < (1 << 23) || ux >= 0x7f800000u))
+    {
+      if (ux == 0 || ux >= 0x7f800000u)
+	return as_special (x);
+      /* subnormal */
+      int n = __builtin_clz (ux) - 8;
+      ux <<= n;
+      ux -= n << 23;
+    }
+  unsigned m = ux & ((1 << 23) - 1), j = (m + (1 << (23 - 7))) >> (23 - 6);
+  double ix = tr[j], l = tl[j];
+  int e = ((int) ux >> 23) - 127;
+  unsigned je = e + 1;
+  je = (je * 0x4d104d4) >> 28;
+  if (__glibc_unlikely (ux == st[je].u))
+    return je;
+
+  double tz = asdouble (((int64_t) m | ((int64_t) 1023 << 23)) << (52 - 23));
+  double z = tz * ix - 1, z2 = z * z;
+  double r
+      = ((e * 0x1.34413509f79ffp-2 + l) + z * b[0]) + z2 * (b[1] + z * b[2]);
+  float ub = r, lb = r + 0x1.b008p-34;
+  if (__glibc_unlikely (ub != lb))
+    {
+      double f = z
+		 * ((c[0] + z * c[1])
+		    + z2
+			  * ((c[2] + z * c[3])
+			     + z2 * (c[4] + z * c[5] + z2 * c[6])));
+      f -= 0x1.0cee0ed4ca7e9p-54 * e;
+      f += l - tl[0];
+      double el = e * 0x1.34413509f7ap-2;
+      r = el + f;
+      ub = r;
+      tz = r;
+      if (__glibc_unlikely (!((asuint64 (tz) & ((1 << 28) - 1)))))
+	{
+	  double dr = (el - r) + f;
+	  r += dr * 32;
+	  ub = r;
 	}
-	if (__builtin_expect(hx >= 0x7f800000, 0)) return x+x;
-	k += (hx>>23)-127;
-	i  = ((uint32_t)k&0x80000000)>>31;
-	hx = (hx&0x007fffff)|((0x7f-i)<<23);
-	y  = (float)(k+i);
-	if (FIX_INT_FP_CONVERT_ZERO && y == 0.0f)
-	  y = 0.0f;
-	SET_FLOAT_WORD(x,hx);
-	z  = y*log10_2lo + ivln10*__ieee754_logf(x);
-	return  z+y*log10_2hi;
+    }
+  return ub;
 }
 libm_alias_finite (__ieee754_log10f, __log10f)
-- 
2.43.0