The state of glibc libm

The state of glibc libm

[Joseph S. Myers]

I've reviewed many (not yet all) of glibc's open "math" component bugs.  I 
hope some actual summary information on what the current state of libm 
looks like may be of interest to the people involved in the various past 
discussions of better libm for GCC or glibc - and those interested in 
fixing things, whether through patches to existing code, new 
implementations or both.

I would say the actual main issues with the present glibc libm are the 
following (if there are others, they are not well-reflected in the open 
bugs).  When I refer to a function by name I'm generally referred to all 
of the float, double and various long double implementations (x86 
extended, IEEE binary128, IBM double+double) - and fixes will often need 
to fix multiple versions of a function in similar ways.

(a) Most libm functions are not correctly rounded - and do not make an 
attempt at being correctly rounded.

A full fix would likely require new (automatically generated and tuned) 
implementations such as proposed at 
<http://gcc.gnu.org/ml/gcc/2012-02/msg00298.html>.  As I understand it, 
correct rounding (proved correct) is generally feasible for functions of 
one binary32, binary64 or x86 extended argument.  For functions of two 
arguments, or one binary128 argument, it may not be feasible to search for 
worst cases for correct rounding, although it may be possible to produce 
implementations that are "probably" correctly rounding.  For functions of 
complex arguments or IBM long double, correct rounding may be less 
feasible (even complex multiplication and division, and all of +-*/ on IBM 
long double, are not correctly rounding, and correct rounding isn't so 
well-defined for IBM long double with its possibility of discontiguous 
mantissa bits).

Known inaccuracies in functions are indicated by the libm-test-ulps files 
in glibc.  Given my patch 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00770.html> to reduce some 
old figures that look like having been left behind after tests or the 
library were fixed, all listed errors are 24ulp or below (9ulp or below 
for x86 and x86_64; I haven't tested the larger errors on other 
architectures to see if they are actually still applicable).

(b) Where functions do make attempts at being correctly rounded 
(especially the IBM Accurate Mathematical Library functions), they tend to 
be sufficiently slow that the slowness attracts bug reports.  Again, this 
would likely be addressed by new implementations that use careful error 
bounds and information about worst cases to reduce the cost of being 
correctly rounding.

(c) Various functions do not set errno correctly (many cases) or raise the 
proper floating-point exceptions (a smaller number of cases - both 
spurious exceptions where not permitted by ISO C, and failing to raise 
required overflow/underflow exceptions).  In general this is a separate 
bug for each function (filed as many separate bugs in glibc Bugzilla) and 
can be fixed by a separate local patch for each function (adding a 
testcase, of course - note that glibc's main libm-test.inc presently only 
tests invalid and divide-by-zero exceptions, so if working on these error 
handling issues it might be useful to extend it to cover other exceptions 
as well as errno values).

(d) There are some specific bugs filed for functions such as nexttoward 
whose results are precisely specified by ISO C; in general these should be 
fixable by local patches.

(e) Various functions, mainly IBM Accurate Mathematical Library ones, 
produce wildly wrong results or crash in non-default rounding modes.  I 
have a patch for exp pending review at 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00748.html> and I expect 
others can be fixed similarly.

(f) Various trigonometrical functions are inaccurate on x86 and x86_64 
(see glibc bug 13658 and recent discussions).

(g) Bessel function implementations handle large inputs in different ways 
to other functions, as I discuss at 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00512.html>.

(h) Various complex functions have problems such as inaccuracy or wrong 
branch cuts.

(i) Some real functions have particular issues (which should be fixable by 
local changes short of new correctly rounded implementations):

  - erfc (my patch 
    <http://sourceware.org/ml/libc-alpha/2012-02/msg00640.html> is pending 
    review).

  - pow (bugs 369, 706, 2678, 3866).  The assembly implementations may 
    complicate fixing these issues, though it's probably possible to fix 
    only some bugs (in all relevant implementations, with plenty of 
    testcases) rather than a patch needing to fix both all issues and all 
    implementations at once.

  - lgamma, in cases where the result is close to 0 and there is a lot of 
    cancellation in the present calculations.

  - tgamma, in cases of results of large magnitude (where the approach of 
    using exp (lgamma) leads to large errors).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[David Miller]

From: "Joseph S. Myers" <joseph@codesourcery.com>
Date: Wed, 29 Feb 2012 17:17:17 +0000 (UTC)

Thanks for looking into all of these issues.

> (c) Various functions do not set errno correctly (many cases) or raise the 
> proper floating-point exceptions (a smaller number of cases - both 
> spurious exceptions where not permitted by ISO C, and failing to raise 
> required overflow/underflow exceptions).  In general this is a separate 
> bug for each function (filed as many separate bugs in glibc Bugzilla) and 
> can be fixed by a separate local patch for each function (adding a 
> testcase, of course - note that glibc's main libm-test.inc presently only 
> tests invalid and divide-by-zero exceptions, so if working on these error 
> handling issues it might be useful to extend it to cover other exceptions 
> as well as errno values).

This reminds me that there are math tests local to the powerpc port
that therefore only run on powerpc.  For example, I've looked at
sysdeps/powerpc/test-arith{,f}.c and they don't seem so non-portable
that we couldn't run them everywhere with just some small tweaks.

Re: The state of glibc libm

[Vincent Lefevre]

On 2012-02-29 17:17:17 +0000, Joseph S. Myers wrote:
> (a) Most libm functions are not correctly rounded - and do not make an 
> attempt at being correctly rounded.
> 
> A full fix would likely require new (automatically generated and tuned) 
> implementations such as proposed at 
> <http://gcc.gnu.org/ml/gcc/2012-02/msg00298.html>.  As I understand it, 
> correct rounding (proved correct) is generally feasible for functions of 
> one binary32, binary64 or x86 extended argument.  For functions of two 
> arguments, or one binary128 argument, it may not be feasible to search for 
> worst cases for correct rounding, although it may be possible to produce 
> implementations that are "probably" correctly rounding.  For functions of 
> complex arguments or IBM long double, correct rounding may be less 
> feasible (even complex multiplication and division, and all of +-*/ on IBM 
> long double, are not correctly rounding, and correct rounding isn't so 
> well-defined for IBM long double with its possibility of discontiguous 
> mantissa bits).

For double-double (IBM long double), I don't think the notion of
correct rounding makes much sense anyway. Actually the double-double
arithmetic is mainly useful for the basic operations in order to be
able to implement elementary functions accurately (first step in
Ziv's strategy, possibly a second step as well). IMHO, on such a
platform, if expl() (for instance) just calls exp(), this is OK.

> (b) Where functions do make attempts at being correctly rounded 
> (especially the IBM Accurate Mathematical Library functions), they tend to 
> be sufficiently slow that the slowness attracts bug reports.  Again, this 
> would likely be addressed by new implementations that use careful error 
> bounds and information about worst cases to reduce the cost of being 
> correctly rounding.

I'm not sure that the complaints are about worst cases. More probably
software implementation vs hardware implementation in the average
case. But a new software implementation (better in average) could
help.

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: The state of glibc libm

[Joseph S. Myers]

On Wed, 14 Mar 2012, Vincent Lefevre wrote:

> For double-double (IBM long double), I don't think the notion of
> correct rounding makes much sense anyway. Actually the double-double
> arithmetic is mainly useful for the basic operations in order to be
> able to implement elementary functions accurately (first step in
> Ziv's strategy, possibly a second step as well). IMHO, on such a
> platform, if expl() (for instance) just calls exp(), this is OK.

expl just calling exp - losing 53 bits of precision - seems rather 
extreme.  But I'd think it would be fine to say: when asked to compute 
f(x), take x' within 10ulp of x, and return a number within 10ulp of 
f(x'), where ulp is interpreted as if the mantissa were a fixed 106 bits 
(fewer bits for subnormals, of course).  (And as a consequence, accurate 
range reduction for large arguments would be considered not to matter for 
IBM long double; sin and cos could return any value in the range [-1, 1] 
for sufficiently large arguments.)

> > (b) Where functions do make attempts at being correctly rounded 
> > (especially the IBM Accurate Mathematical Library functions), they tend to 
> > be sufficiently slow that the slowness attracts bug reports.  Again, this 
> > would likely be addressed by new implementations that use careful error 
> > bounds and information about worst cases to reduce the cost of being 
> > correctly rounding.
> 
> I'm not sure that the complaints are about worst cases. More probably
> software implementation vs hardware implementation in the average
> case. But a new software implementation (better in average) could
> help.

Various bugs do complain about particular cases being slow (as well as 
about such things as sinf being slower than sin - there, if you 
automatically generate functions based not just on the type for the 
function being generated but also on what wider types are available and 
efficient in hardware, you could generate a version of sinf that uses 
double or long double computations internally to speed things up).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Jeff Law]

On 03/14/2012 08:30 AM, Vincent Lefevre wrote:
>
>> (b) Where functions do make attempts at being correctly rounded
>> (especially the IBM Accurate Mathematical Library functions), they tend to
>> be sufficiently slow that the slowness attracts bug reports.  Again, this
>> would likely be addressed by new implementations that use careful error
>> bounds and information about worst cases to reduce the cost of being
>> correctly rounding.
>
> I'm not sure that the complaints are about worst cases. More probably
> software implementation vs hardware implementation in the average
> case. But a new software implementation (better in average) could
> help.
The complaints I typically see are about worst case performance.  I 
occasionally see requests for better performance with the potential loss 
of accuracy.

jeff

Re: The state of glibc libm

[Joseph S. Myers]

On Wed, 14 Mar 2012, Jeff Law wrote:

> On 03/14/2012 08:30 AM, Vincent Lefevre wrote:
> > 
> > > (b) Where functions do make attempts at being correctly rounded
> > > (especially the IBM Accurate Mathematical Library functions), they tend to
> > > be sufficiently slow that the slowness attracts bug reports.  Again, this
> > > would likely be addressed by new implementations that use careful error
> > > bounds and information about worst cases to reduce the cost of being
> > > correctly rounding.
> > 
> > I'm not sure that the complaints are about worst cases. More probably
> > software implementation vs hardware implementation in the average
> > case. But a new software implementation (better in average) could
> > help.
> The complaints I typically see are about worst case performance.  I
> occasionally see requests for better performance with the potential loss of
> accuracy.

I'd say that "better performance with the potential loss of accuracy" 
should be covered by -ffast-math - that GCC should generate direct use of 
fsin/fcos instructions for sin/cos for -O2 -funsafe-math-optimizations on 
x86_64, as it does on x86, unless there is some reason to think they would 
perform worse than the out-of-line implementation.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Jeff Law]

On 03/14/2012 10:30 AM, Joseph S. Myers wrote:
>
> I'd say that "better performance with the potential loss of accuracy"
> should be covered by -ffast-math - that GCC should generate direct use of
> fsin/fcos instructions for sin/cos for -O2 -funsafe-math-optimizations on
> x86_64, as it does on x86, unless there is some reason to think they would
> perform worse than the out-of-line implementation.
The better performance with potential loss of accuracy is an across the 
board request, it's not just a sin/cos issue.  All the trig, exp, pow, 
log, etc, which I don't think are necessarily covered by using the old 
x87 fp unit.

Jeff

Re: The state of glibc libm

[Andi Kleen]

Jeff Law <law@redhat.com> writes:

> On 03/14/2012 10:30 AM, Joseph S. Myers wrote:
>>
>> I'd say that "better performance with the potential loss of accuracy"
>> should be covered by -ffast-math - that GCC should generate direct use of
>> fsin/fcos instructions for sin/cos for -O2 -funsafe-math-optimizations on
>> x86_64, as it does on x86, unless there is some reason to think they would
>> perform worse than the out-of-line implementation.
> The better performance with potential loss of accuracy is an across
> the board request, it's not just a sin/cos issue.  All the trig, exp,
> pow, log, etc, which I don't think are necessarily covered by using
> the old x87 fp unit.

Note that various interesting CPUs (like Atom) have a very slow x87
units and it's generally deprecated and discouraged.

So if you go for a new libm it's better to be SSE only (ideally 
AVX too) and avoid any usage of x87.

One big win alone on 32bit x86 would be to use a SSE ABI for libm
by default.

You may still need some fallback for old CPUs, but this doesn't
need to be terrible optimized.

-Andi
-- 
ak@linux.intel.com -- Speaking for myself only

Re: The state of glibc libm

[Marc Glisse]

On Wed, 14 Mar 2012, Joseph S. Myers wrote:

> I'd say that "better performance with the potential loss of accuracy"
> should be covered by -ffast-math - that GCC should generate direct use of
> fsin/fcos instructions for sin/cos for -O2 -funsafe-math-optimizations on
> x86_64, as it does on x86, unless there is some reason to think they would
> perform worse than the out-of-line implementation.

Last time I did some timings (maybe 4 years ago), for double, fsin was 
slower than the libm software implementation compiled for x87, which was 
itself slower than the same implementation compiled for sse. And the 
software implementation was more precise than fsin. My conclusion was to 
ignore fsin from then on.

-- 
Marc Glisse

Re: The state of glibc libm

[Joseph S. Myers]

On Wed, 14 Mar 2012, Andi Kleen wrote:

> One big win alone on 32bit x86 would be to use a SSE ABI for libm
> by default.

I haven't checked, but I'd hope x32 does that as a better 32-bit ABI for 
newer x86 processors.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Joseph S. Myers]

On Wed, 14 Mar 2012, Marc Glisse wrote:

> On Wed, 14 Mar 2012, Joseph S. Myers wrote:
> 
> > I'd say that "better performance with the potential loss of accuracy"
> > should be covered by -ffast-math - that GCC should generate direct use of
> > fsin/fcos instructions for sin/cos for -O2 -funsafe-math-optimizations on
> > x86_64, as it does on x86, unless there is some reason to think they would
> > perform worse than the out-of-line implementation.
> 
> Last time I did some timings (maybe 4 years ago), for double, fsin was slower
> than the libm software implementation compiled for x87, which was itself
> slower than the same implementation compiled for sse. And the software
> implementation was more precise than fsin. My conclusion was to ignore fsin
> from then on.

Interesting - hopefully that means the glibc changes Andreas Jaeger and I 
have been working on to stop using fsin etc. (in the interests of accuracy 
and fixing the sincos issues discussed lately) won't actually make 
performance worse.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Andi Kleen]

On Wed, Mar 14, 2012 at 09:04:53PM +0000, Joseph S. Myers wrote:
> On Wed, 14 Mar 2012, Andi Kleen wrote:
> 
> > One big win alone on 32bit x86 would be to use a SSE ABI for libm
> > by default.
> 
> I haven't checked, but I'd hope x32 does that as a better 32-bit ABI for 
> newer x86 processors.

x32 is not a mainstream option.

-Andi
-- 
ak@linux.intel.com -- Speaking for myself only.

Re: The state of glibc libm

[David Miller]

From: Jeff Law <law@redhat.com>
Date: Wed, 14 Mar 2012 10:41:27 -0600

> The better performance with potential loss of accuracy is an across
> the board request, it's not just a sin/cos issue.  All the trig, exp,
> pow, log, etc, which I don't think are necessarily covered by using
> the old x87 fp unit.

This is no surprise to me, vendors have commonly provided some kind of
"optimized libm" library that provided better performance with less
accuracy.

For example, Sun has always had a variant of libm optimized in this
way, called "libmopt", with it's SunPRO compilers.

Re: The state of glibc libm

[Vincent Lefevre]

On 2012-03-14 14:40:06 +0000, Joseph S. Myers wrote:
> On Wed, 14 Mar 2012, Vincent Lefevre wrote:
> 
> > For double-double (IBM long double), I don't think the notion of
> > correct rounding makes much sense anyway. Actually the double-double
> > arithmetic is mainly useful for the basic operations in order to be
> > able to implement elementary functions accurately (first step in
> > Ziv's strategy, possibly a second step as well). IMHO, on such a
> > platform, if expl() (for instance) just calls exp(), this is OK.
> 
> expl just calling exp - losing 53 bits of precision - seems rather 
> extreme.  But I'd think it would be fine to say: when asked to compute 
> f(x), take x' within 10ulp of x, and return a number within 10ulp of 
> f(x'), where ulp is interpreted as if the mantissa were a fixed 106 bits 
> (fewer bits for subnormals, of course).  (And as a consequence, accurate 
> range reduction for large arguments would be considered not to matter for 
> IBM long double; sin and cos could return any value in the range [-1, 1] 
> for sufficiently large arguments.)

After thinking about this, you could assume that you have a 106-bit
floating-point system (BTW, LDBL_MANT_DIG = 106) and use the same
method to generate code that provides an accurate implementation
(if the code generator doesn't assume an IEEE 754 compatible FP
system). Concerning sin and cos, I think there should be a minimum
of specification and some consistency (such as sin(x)² + cos(x)²
being close to 1).

> Various bugs do complain about particular cases being slow (as well as 
> about such things as sinf being slower than sin - there, if you 
> automatically generate functions based not just on the type for the 
> function being generated but also on what wider types are available and 
> efficient in hardware, you could generate a version of sinf that uses 
> double or long double computations internally to speed things up).

sinf being slower than sin is surprising (but I know that sinl
could be faster than sin on x86_64 because sin uses the accurate
IBM implementation, while sinl uses the hardware instruction).

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: The state of glibc libm

[Richard Guenther]

On Wed, Mar 14, 2012 at 10:47 PM, Andi Kleen <andi@firstfloor.org> wrote:
> On Wed, Mar 14, 2012 at 09:04:53PM +0000, Joseph S. Myers wrote:
>> On Wed, 14 Mar 2012, Andi Kleen wrote:
>>
>> > One big win alone on 32bit x86 would be to use a SSE ABI for libm
>> > by default.
>>
>> I haven't checked, but I'd hope x32 does that as a better 32-bit ABI for
>> newer x86 processors.
>
> x32 is not a mainstream option.

SSE ABI entries for i?86 in glibc were rejected.  I proposed them like
4-5 years ago to make -mfpmath=sse not suck.

Richard.

> -Andi
> --
> ak@linux.intel.com -- Speaking for myself only.

Re: The state of glibc libm

[Andi Kleen]

> SSE ABI entries for i?86 in glibc were rejected.  I proposed them like
> 4-5 years ago to make -mfpmath=sse not suck.

With the new libm hopefully this can be revisited.

But there's the ABI and there's the internal implementation.

My point was just that relying on x87 fully again does not really make
sense anymore in 2012.

-andi

-- 
ak@linux.intel.com -- Speaking for myself only.

Re: The state of glibc libm

[Richard Guenther]

On Thu, Mar 15, 2012 at 3:17 PM, Andi Kleen <andi@firstfloor.org> wrote:
>> SSE ABI entries for i?86 in glibc were rejected.  I proposed them like
>> 4-5 years ago to make -mfpmath=sse not suck.
>
> With the new libm hopefully this can be revisited.
>
> But there's the ABI and there's the internal implementation.
>
> My point was just that relying on x87 fully again does not really make
> sense anymore in 2012.

That's true, the return value shuffle to/from the x87 stack might not be
too bad for performance, nor argument passing on the stack.

Richard.

Re: The state of glibc libm

[James Cloos]

>>>>> "VL" == Vincent Lefevre <vincent+gcc@vinc17.org> writes:

VL> sinf being slower than sin is surprising

Some weeks back I wrote a simple app to output a couple of waveforms as
float samples (I encased it in a tiny script which piped the output to
sox to create a wav).  I found that converting it to work on and output
doubles improved performance from around 2.2 sec per minute of output to
around 0.3 sec per minute for a simple tan(sin(t)) waveform.

That is using glibc-2.14.1 on amd64 on a 2800MHz K10, compiled with
gcc-4.6.2 with -O2 and the graphite -floop-... options.

Some improvements have landed for the float versions of libm since.

-JimC
-- 
James Cloos <cloos@jhcloos.com>         OpenPGP: 1024D/ED7DAEA6

Re: The state of glibc libm

[Steven Munroe]

On Thu, 2012-03-15 at 03:07 +0100, Vincent Lefevre wrote:
> On 2012-03-14 14:40:06 +0000, Joseph S. Myers wrote:
> > On Wed, 14 Mar 2012, Vincent Lefevre wrote:
> > 
> > > For double-double (IBM long double), I don't think the notion of
> > > correct rounding makes much sense anyway. Actually the double-double
> > > arithmetic is mainly useful for the basic operations in order to be
> > > able to implement elementary functions accurately (first step in
> > > Ziv's strategy, possibly a second step as well). IMHO, on such a
> > > platform, if expl() (for instance) just calls exp(), this is OK.
> > 
Why would that be OK? If we have higher precision long double then the
libm should deliver that higher precision.


> > expl just calling exp - losing 53 bits of precision - seems rather 
> > extreme.  But I'd think it would be fine to say: when asked to compute 
> > f(x), take x' within 10ulp of x, and return a number within 10ulp of 
> > f(x'), where ulp is interpreted as if the mantissa were a fixed 106 bits 
> > (fewer bits for subnormals, of course).  (And as a consequence, accurate 
> > range reduction for large arguments would be considered not to matter for 
> > IBM long double; sin and cos could return any value in the range [-1, 1] 
> > for sufficiently large arguments.)
> 
> After thinking about this, you could assume that you have a 106-bit
> floating-point system (BTW, LDBL_MANT_DIG = 106) and use the same
> method to generate code that provides an accurate implementation
> (if the code generator doesn't assume an IEEE 754 compatible FP
> system). Concerning sin and cos, I think there should be a minimum
> of specification and some consistency (such as sin(x)² + cos(x)²
> being close to 1).
> 
actually back in 2007 I overrode slowexp and slowpow for powerpc/power4
(and later) to use expl, powl on the slow path of exp/pow, instead of
the mpa.h implementation. This provide a nice performance improvement
but does imply some rounding mode issues. 

Re: The state of glibc libm

[Vincent Lefevre]

On 2012-03-15 10:52:05 -0500, Steven Munroe wrote:
> On Thu, 2012-03-15 at 03:07 +0100, Vincent Lefevre wrote:
> > On 2012-03-14 14:40:06 +0000, Joseph S. Myers wrote:
> > > On Wed, 14 Mar 2012, Vincent Lefevre wrote:
> > > 
> > > > For double-double (IBM long double), I don't think the notion of
> > > > correct rounding makes much sense anyway. Actually the double-double
> > > > arithmetic is mainly useful for the basic operations in order to be
> > > > able to implement elementary functions accurately (first step in
> > > > Ziv's strategy, possibly a second step as well). IMHO, on such a
> > > > platform, if expl() (for instance) just calls exp(), this is OK.
> > > 
> Why would that be OK? If we have higher precision long double then the
> libm should deliver that higher precision.

I initially thought that the only goal of a double-double format
(instead of the standard quadruple precision) was to get an
accurate implementation of the elementary functions in double
precision (BTW, that's probably why expl() and so on didn't
exist before C99).

Now, since expl() now exists, if the user calls it, perhaps his
goal is to get more precision, so finally I agree that expl()
should really have an accuracy close to LDBL_MANT_DIG. However
this is quite useless in portable programs, where long double
can have the same precision as double (as this is the case on
ARM).

For the same reason, if the user chose long double instead of
double, this may be because he wanted more precision than double.
So, in the long term, the ABI should probably be changed to have
long double = quadruple precision (binary128).

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: The state of glibc libm

[Joseph S. Myers]

On Thu, 22 Mar 2012, Vincent Lefevre wrote:

> I initially thought that the only goal of a double-double format
> (instead of the standard quadruple precision) was to get an
> accurate implementation of the elementary functions in double
> precision (BTW, that's probably why expl() and so on didn't
> exist before C99).

I'd imagine the point is to be faster than pure software quad-precision 
floating point.

> Now, since expl() now exists, if the user calls it, perhaps his
> goal is to get more precision, so finally I agree that expl()
> should really have an accuracy close to LDBL_MANT_DIG. However
> this is quite useless in portable programs, where long double
> can have the same precision as double (as this is the case on
> ARM).

(FWIW, the ABI for AArch64 - ARMv8 in 64-bit mode - uses quad precision 
for long double.)

> For the same reason, if the user chose long double instead of
> double, this may be because he wanted more precision than double.

You mean range?  IBM long double provides more precision, but not more 
range.

> So, in the long term, the ABI should probably be changed to have
> long double = quadruple precision (binary128).

The ABI for Power Architecture changed away from quad precision to using 
IBM long double (the original SysV ABI for PowerPC used quad precision, 
the current ABI uses IBM long double)....

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Vincent Lefevre]

On 2012-03-22 16:29:00 +0000, Joseph S. Myers wrote:
> On Thu, 22 Mar 2012, Vincent Lefevre wrote:
> > For the same reason, if the user chose long double instead of
> > double, this may be because he wanted more precision than double.
> 
> You mean range?  IBM long double provides more precision, but not more 
> range.

Well, precision and/or range. If double precision format is sufficient
for his application, the user can just choose the "double" type. So,
I don't think that it is useful to have long double = double.

Then concerning double-double vs quad (binary128) for the "long double"
type, I think that quad would be more useful, in particular because
it has been standardized and it is a true FP format. If need be (for
efficiency reasons), double-double could still be implemented using
the "double" type, via a library or ad-hoc code (that does something
more clever, taking the context into account). And the same code (with
just a change of the base type) could be reused to get a double-quad
(i.e. quad + quad) arithmetic, that can be useful to implement the
"long double" versions of the math functions (expl, and so on).

> > So, in the long term, the ABI should probably be changed to have
> > long double = quadruple precision (binary128).
> 
> The ABI for Power Architecture changed away from quad precision to using 
> IBM long double (the original SysV ABI for PowerPC used quad precision, 
> the current ABI uses IBM long double)....

Perhaps they could change back to quad precision.

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: The state of glibc libm

[Steven Munroe]

On Mon, 2012-03-26 at 12:26 +0200, Vincent Lefevre wrote:
> On 2012-03-22 16:29:00 +0000, Joseph S. Myers wrote:
> > On Thu, 22 Mar 2012, Vincent Lefevre wrote:
> > > For the same reason, if the user chose long double instead of
> > > double, this may be because he wanted more precision than double.
> > 
> > You mean range?  IBM long double provides more precision, but not more 
> > range.
> 
> Well, precision and/or range. If double precision format is sufficient
> for his application, the user can just choose the "double" type. So,
> I don't think that it is useful to have long double = double.
> 
> Then concerning double-double vs quad (binary128) for the "long double"
> type, I think that quad would be more useful, in particular because
> it has been standardized and it is a true FP format. If need be (for
> efficiency reasons), double-double could still be implemented using
> the "double" type, via a library or ad-hoc code (that does something
> more clever, taking the context into account). And the same code (with
> just a change of the base type) could be reused to get a double-quad
> (i.e. quad + quad) arithmetic, that can be useful to implement the
> "long double" versions of the math functions (expl, and so on).
> 
This is much easier said then done. In practice it is a major ABI change
and would have to be staged over multiple (7-10) years.

> > > So, in the long term, the ABI should probably be changed to have
> > > long double = quadruple precision (binary128).
> > 
> > The ABI for Power Architecture changed away from quad precision to using 
> > IBM long double (the original SysV ABI for PowerPC used quad precision, 
> > the current ABI uses IBM long double)....
> 
> Perhaps they could change back to quad precision.
> 
That is not the feedback we get from our customers. No one will use
software IEEE binary128 and we don't have hardware binary128. So far
there is abstract interest but no strong demand for this. So there is no
incentive to change.

Re: The state of glibc libm

[Vincent Lefevre]

On 2012-03-26 11:15:37 -0500, Steven Munroe wrote:
> On Mon, 2012-03-26 at 12:26 +0200, Vincent Lefevre wrote:
> > Then concerning double-double vs quad (binary128) for the "long double"
> > type, I think that quad would be more useful, in particular because
> > it has been standardized and it is a true FP format. If need be (for
> > efficiency reasons), double-double could still be implemented using
> > the "double" type, via a library or ad-hoc code (that does something
> > more clever, taking the context into account). And the same code (with
> > just a change of the base type) could be reused to get a double-quad
> > (i.e. quad + quad) arithmetic, that can be useful to implement the
> > "long double" versions of the math functions (expl, and so on).
> > 
> This is much easier said then done. In practice it is a major ABI change
> and would have to be staged over multiple (7-10) years.

Yes, I meant in the long term.

> > > > So, in the long term, the ABI should probably be changed to have
> > > > long double = quadruple precision (binary128).
> > > 
> > > The ABI for Power Architecture changed away from quad precision to using 
> > > IBM long double (the original SysV ABI for PowerPC used quad precision, 
> > > the current ABI uses IBM long double)....
> > 
> > Perhaps they could change back to quad precision.
> > 
> That is not the feedback we get from our customers. No one will use
> software IEEE binary128 and we don't have hardware binary128. So far
> there is abstract interest but no strong demand for this. So there is no
> incentive to change.

But perhaps one reason why there is no hardware binary128 is that
there is no software that uses it.

BTW, what is the feedback from your customers concerning the
long double math functions (expl, etc.)?

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: The state of glibc libm

[Joseph S. Myers]

On Wed, 29 Feb 2012, Joseph S. Myers wrote:

> I've reviewed many (not yet all) of glibc's open "math" component bugs.  I 
> hope some actual summary information on what the current state of libm 
> looks like may be of interest to the people involved in the various past 
> discussions of better libm for GCC or glibc - and those interested in 
> fixing things, whether through patches to existing code, new 
> implementations or both.

Here is an update on this analysis 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00772.html>, having now 
reviewed all the open "math" component bugs and fixed quite a few of them 
(and filed more bugs for issues found in the fixing process).  Any glibc 
libm bugs you know of that don't currently have an open bug filed in 
Bugzilla should be filed there.

> (a) Most libm functions are not correctly rounded - and do not make an 
> attempt at being correctly rounded.

Still present.

> (b) Where functions do make attempts at being correctly rounded 
> (especially the IBM Accurate Mathematical Library functions), they tend to 
> be sufficiently slow that the slowness attracts bug reports.  Again, this 
> would likely be addressed by new implementations that use careful error 
> bounds and information about worst cases to reduce the cost of being 
> correctly rounding.

Still present.

> (c) Various functions do not set errno correctly (many cases) or raise the 
> proper floating-point exceptions (a smaller number of cases - both 
> spurious exceptions where not permitted by ISO C, and failing to raise 
> required overflow/underflow exceptions).  In general this is a separate 
> bug for each function (filed as many separate bugs in glibc Bugzilla) and 
> can be fixed by a separate local patch for each function (adding a 
> testcase, of course - note that glibc's main libm-test.inc presently only 
> tests invalid and divide-by-zero exceptions, so if working on these error 
> handling issues it might be useful to extend it to cover other exceptions 
> as well as errno values).

Still present.  libm-test.inc now tests for overflow exceptions as well as 
invalid and divide-by-zero.

> (d) There are some specific bugs filed for functions such as nexttoward 
> whose results are precisely specified by ISO C; in general these should be 
> fixable by local patches.

Still present (some fixed, some other such bugs found).

> (e) Various functions, mainly IBM Accurate Mathematical Library ones, 
> produce wildly wrong results or crash in non-default rounding modes.  I 
> have a patch for exp pending review at 
> <http://sourceware.org/ml/libc-alpha/2012-02/msg00748.html> and I expect 
> others can be fixed similarly.

Fixed.

> (f) Various trigonometrical functions are inaccurate on x86 and x86_64 
> (see glibc bug 13658 and recent discussions).

Fixed.

> (g) Bessel function implementations handle large inputs in different ways 
> to other functions, as I discuss at 
> <http://sourceware.org/ml/libc-alpha/2012-02/msg00512.html>.

Fixed.

> (h) Various complex functions have problems such as inaccuracy or wrong 
> branch cuts.

Some fixed, others still present (though I don't think there are any wrong 
branch cut issues after all).

> (i) Some real functions have particular issues (which should be fixable by 
> local changes short of new correctly rounded implementations):
> 
>   - erfc (my patch 
>     <http://sourceware.org/ml/libc-alpha/2012-02/msg00640.html> is pending 
>     review).

Fixed.

>   - pow (bugs 369, 706, 2678, 3866).  The assembly implementations may 
>     complicate fixing these issues, though it's probably possible to fix 
>     only some bugs (in all relevant implementations, with plenty of 
>     testcases) rather than a patch needing to fix both all issues and all 
>     implementations at once.

Parts fixed, parts pending (bugs 706 and 13881).

>   - lgamma, in cases where the result is close to 0 and there is a lot of 
>     cancellation in the present calculations.
> 
>   - tgamma, in cases of results of large magnitude (where the approach of 
>     using exp (lgamma) leads to large errors).

Still present.  Some other issues with particular functions have also been 
found since then.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: The state of glibc libm

[Uros Bizjak]

Hello!

>>> SSE ABI entries for i?86 in glibc were rejected. ?I proposed them like
>>> 4-5 years ago to make -mfpmath=sse not suck.
>>
>> With the new libm hopefully this can be revisited.
>>
>> But there's the ABI and there's the internal implementation.
>>
>> My point was just that relying on x87 fully again does not really make
>> sense anymore in 2012.
>
> That's true, the return value shuffle to/from the x87 stack might not be
> too bad for performance, nor argument passing on the stack.

Can this issue be solved with alternative function entry points, in
the same way ICC does? This way, we can pass values in SSE, not on the
stack.

Uros.