libm accuracy, eglibc compared to glibc

libm accuracy, eglibc compared to glibc

[Mats Kärrman]

Hi,

I face a problem with libm. With my Dora build the following assertion doesn't hold:

sqrt( pow( sqrt( 2.0 ), 4.0 ) ) == 2.0

1) With my old OE-classic / glibc-2.9 / PowerPC-hf it holds.
2) With my Debian desktop PC / glibc-2.17 / amd64 it holds.
3) With my new OE-core Dora / eglibc-2.18 / PowerPC-hf it doesn't hold.

Close enough one can think but the standardized test case I run requires the result to be 2.0.
A simple test program:
--------------------------------------------------------------------
#include <stdio.h>
#include <math.h>

int main(){

	double two = 2.0;
	double four = 4.0;

	double s1 = sqrt(two);
	printf("sqrt(%.20f) = %.20f\n", two, s1);
	double p   = pow(s1, four);
	printf("pow(%.20f, %.20f) = %.20f\n", s1, four, p);
	double s2 = sqrt(p);
	printf("sqrt(%.20f) = %.20f\n", p, s2);

}
--------------------------------------------------------------------

results from 1 and 2 (identical):
sqrt(2.00000000000000000000) = 1.41421356237309514547
pow(1.41421356237309514547, 4.00000000000000000000) = 4.00000000000000088818
sqrt(4.00000000000000088818) = 2.00000000000000000000

results from 3:
sqrt(2.00000000000000000000) = 1.41421356237309492343
pow(1.41421356237309492343, 4.00000000000000000000) = 3.99999999999999866773
sqrt(3.99999999999999866773) = 1.99999999999999955591

Does anyone know if this is a known "feature" of eglibc or know any other reason
for this difference?

Best Regards,
Mats

Re: libm accuracy, eglibc compared to glibc

[Stanacar, StefanX]

Hello,

I just want to point out that this looks very similar to a qemuppc
machine bug we've seen before with libm and floor() call:
https://bugzilla.yoctoproject.org/show_bug.cgi?id=4854

It was fixed with:
http://git.yoctoproject.org/cgit/cgit.cgi/poky/commit/?id=2a48e9007cc92fedfdb6919411859fb943eacc94

HTH,
Stefan

On Wed, 2014-03-12 at 15:30 +0000, Mats Kärrman wrote:
> Hi,
> 
> I face a problem with libm. With my Dora build the following assertion doesn't hold:
> 
> sqrt( pow( sqrt( 2.0 ), 4.0 ) ) == 2.0
> 
> 1) With my old OE-classic / glibc-2.9 / PowerPC-hf it holds.
> 2) With my Debian desktop PC / glibc-2.17 / amd64 it holds.
> 3) With my new OE-core Dora / eglibc-2.18 / PowerPC-hf it doesn't hold.
> 
> Close enough one can think but the standardized test case I run requires the result to be 2.0.
> A simple test program:
> --------------------------------------------------------------------
> #include <stdio.h>
> #include <math.h>
> 
> int main(){
> 
> 	double two = 2.0;
> 	double four = 4.0;
> 
> 	double s1 = sqrt(two);
> 	printf("sqrt(%.20f) = %.20f\n", two, s1);
> 	double p   = pow(s1, four);
> 	printf("pow(%.20f, %.20f) = %.20f\n", s1, four, p);
> 	double s2 = sqrt(p);
> 	printf("sqrt(%.20f) = %.20f\n", p, s2);
> 
> }
> --------------------------------------------------------------------
> 
> results from 1 and 2 (identical):
> sqrt(2.00000000000000000000) = 1.41421356237309514547
> pow(1.41421356237309514547, 4.00000000000000000000) = 4.00000000000000088818
> sqrt(4.00000000000000088818) = 2.00000000000000000000
> 
> results from 3:
> sqrt(2.00000000000000000000) = 1.41421356237309492343
> pow(1.41421356237309492343, 4.00000000000000000000) = 3.99999999999999866773
> sqrt(3.99999999999999866773) = 1.99999999999999955591
> 
> Does anyone know if this is a known "feature" of eglibc or know any other reason
> for this difference?
> 
> Best Regards,
> Mats

Re: libm accuracy, eglibc compared to glibc

[Burton, Ross]

On 12 March 2014 15:30, Mats Kärrman <Mats.Karrman@tritech.se> wrote:
> Does anyone know if this is a known "feature" of eglibc or know any other reason
> for this difference?

I expect the standard response from upstream here would be that
floating point mathematics isn't exact and if you want exact answers
then you should use a different representation. The differences start
at the first sqrt():

Test 1 and 2:
sqrt(2.00000000000000000000) = 1.41421356237309514547
Test 3:
sqrt(2.00000000000000000000) = 1.41421356237309492343
Test using bc and scale=30 (bc is fixed-point not floating-point):
sqrt(2.00000000000000000000) = 1.414213562373095048801688724209

The output from bc is a truncation of the "real" value of sqrt(2)
whereas both values from C's sqrt() function using floats are
approximations, both equally close.  Why does your test case expect
the answer to be 2.0?  That's a massive assumption to make, floats can
trivially lose accuracy rapidly as
http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems
explains.

Ross

Re: libm accuracy, eglibc compared to glibc

[Phil Blundell]

On Wed, 2014-03-12 at 15:30 +0000, Mats Kärrman wrote:
> Does anyone know if this is a known "feature" of eglibc or know any other reason
> for this difference?

Does eglibc's own libm testsuite pass on your platform?  I wonder if
oe-core is accidentally building eglibc without "big" libm for you, or
something like that.  Might be worth checking your option-groups.config
to make sure.

p.

Re: libm accuracy, eglibc compared to glibc

[Khem Raj]

Hi Mats

On Wed, Mar 12, 2014 at 8:30 AM, Mats Kärrman <Mats.Karrman@tritech.se> wrote:
> Hi,
>
> I face a problem with libm. With my Dora build the following assertion doesn't hold:
>
> sqrt( pow( sqrt( 2.0 ), 4.0 ) ) == 2.0
>
> 1) With my old OE-classic / glibc-2.9 / PowerPC-hf it holds.
> 2) With my Debian desktop PC / glibc-2.17 / amd64 it holds.
> 3) With my new OE-core Dora / eglibc-2.18 / PowerPC-hf it doesn't hold.
>
> Close enough one can think but the standardized test case I run requires the result to be 2.0.
> A simple test program:
> --------------------------------------------------------------------
> #include <stdio.h>
> #include <math.h>
>
> int main(){
>
>         double two = 2.0;
>         double four = 4.0;
>
>         double s1 = sqrt(two);
>         printf("sqrt(%.20f) = %.20f\n", two, s1);
>         double p   = pow(s1, four);
>         printf("pow(%.20f, %.20f) = %.20f\n", s1, four, p);
>         double s2 = sqrt(p);
>         printf("sqrt(%.20f) = %.20f\n", p, s2);
>
> }
> --------------------------------------------------------------------
>
> results from 1 and 2 (identical):
> sqrt(2.00000000000000000000) = 1.41421356237309514547
> pow(1.41421356237309514547, 4.00000000000000000000) = 4.00000000000000088818
> sqrt(4.00000000000000088818) = 2.00000000000000000000
>
> results from 3:
> sqrt(2.00000000000000000000) = 1.41421356237309492343
> pow(1.41421356237309492343, 4.00000000000000000000) = 3.99999999999999866773
> sqrt(3.99999999999999866773) = 1.99999999999999955591
>
> Does anyone know if this is a known "feature" of eglibc or know any other reason
> for this difference?
>
We have done changes since OE-Classic days in the FPU area for ppc.
what ppc machine is it ? I have access to p2020 based machine (e500v2)
and it works as expected with
dora/eglibc-2.18

> Best Regards,
> Mats
> --
> _______________________________________________
> Openembedded-core mailing list
> Openembedded-core@lists.openembedded.org
> http://lists.openembedded.org/mailman/listinfo/openembedded-core

Re: libm accuracy, eglibc compared to glibc

[Mats Kärrman]

On Wednesday, March 12, 2014 4:58 PM, Burton, Ross wrote:
> On 12 March 2014 15:30, Mats Kärrman <Mats.Karrman@tritech.se> wrote:
>> Does anyone know if this is a known "feature" of eglibc or know any other reason
>> for this difference?
>
> I expect the standard response from upstream here would be that
> floating point mathematics isn't exact and if you want exact answers
> then you should use a different representation. The differences start
> at the first sqrt():
> 
> Test 1 and 2:
> sqrt(2.00000000000000000000) = 1.41421356237309514547
> Test 3:
> sqrt(2.00000000000000000000) = 1.41421356237309492343
> Test using bc and scale=30 (bc is fixed-point not floating-point):
> sqrt(2.00000000000000000000) = 1.414213562373095048801688724209
> 
> The output from bc is a truncation of the "real" value of sqrt(2)
> whereas both values from C's sqrt() function using floats are
> approximations, both equally close.  Why does your test case expect
> the answer to be 2.0?  That's a massive assumption to make, floats can
> trivially lose accuracy rapidly as
> http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems
> explains.

I am aware of the problems with floating point math. The test cases are from IBM
and used to validate their J9 JVM on new platforms. AfaIk they normally work (as
they do for my old OE-classic system and for my Debian PC).

Any reason that EGLIBC should be different from GLIBC?

// Mats

Re: libm accuracy, eglibc compared to glibc

[Mats Kärrman]

On Wednesday, March 12, 2014 6:37 PM, Phil Blundell wrote:
> On Wed, 2014-03-12 at 15:30 +0000, Mats Kärrman wrote:
>> Does anyone know if this is a known "feature" of eglibc or know any other reason
>> for this difference?
>
> Does eglibc's own libm testsuite pass on your platform?  I wonder if
> oe-core is accidentally building eglibc without "big" libm for you, or
> something like that.  Might be worth checking your option-groups.config
> to make sure.

Thanks for the tip, I'll have a look.
// Mats

Re: libm accuracy, eglibc compared to glibc

[Mats Kärrman]

On Wednesday, March 12, 2014 7:09 PM, Khem Raj wrote:
> On Wed, Mar 12, 2014 at 8:30 AM, Mats Kärrman <Mats.Karrman@tritech.se> wrote:
>> I face a problem with libm. With my Dora build the following assertion doesn't hold:
>>
>> sqrt( pow( sqrt( 2.0 ), 4.0 ) ) == 2.0
>>
<snip>
>> Does anyone know if this is a known "feature" of eglibc or know any other reason
>> for this difference?
>>
> We have done changes since OE-Classic days in the FPU area for ppc.
> what ppc machine is it ? I have access to p2020 based machine (e500v2)
> and it works as expected with
> dora/eglibc-2.18

My machine is an e300c4.
I will try running the libm internal tests on my machine as Phil Blundell suggested.

BR // Mats

Re: libm accuracy, eglibc compared to glibc

[Mats Kärrman]

On Wednesday, March 12, 2014 6:37 PM, Phil Blundell wrote:
> On Wed, 2014-03-12 at 15:30 +0000, Mats Kärrman wrote:
>> Does anyone know if this is a known "feature" of eglibc or know any other reason
>> for this difference?
>
> Does eglibc's own libm testsuite pass on your platform?  I wonder if
> oe-core is accidentally building eglibc without "big" libm for you, or
> something like that.  Might be worth checking your option-groups.config
> to make sure.

I have no NFS set up so I have not managed to run the eglibc target tests yet.

I ran the test program I mentioned before on an i.MX6 (arm cortexa9-hf) system
built on Poky/Dora (the PowerPC system is bare-bone OE-core/Dora) and got the
correct results there (i.e identical to my PC et.c.).

Comparing:
.../eglibc/2.18-r0/build-powerpc-oe-linux/option-groups.config (my failing system)
.../eglibc/2.18-r0/build-arm-poky-linux-gnueabi/option-groups.config (arm i.MX6)

They are both identical and both have all options set "= Y".

My "home made" hard float tune for PowerPC looks like this:

------------------------------------------------------------------
# Tune for the e300c3 core
require conf/machine/include/tune-ppce300c3.inc

# Use hardware floating point
AVAILTUNES += "ppce300c3hf"
TUNE_FEATURES_tune-ppce300c3hf = "m32 fpu-hard ppce300c3"
TUNE_PKGARCH_tune-ppce300c3hf = "ppce300c3hf"
PACKAGE_EXTRA_ARCHS_tune-ppce300c3hf = "${PACKAGE_EXTRA_ARCHS_tune-powerpc} ppce300c3hf"
DEFAULTTUNE = "ppce300c3hf"
------------------------------------------------------------------

Your comments are welcome.

BR // Mats