Subject: Re: floating-point math problem in Cygwin 1.7?
From: mike marchywka <[email-masked]>
To: [email-masked]
Cc: Linh Phan <[email-masked]>
Date: Fri, 25 Dec 2009 07:47:53 -0500
List-Id: <cygwin.cygwin.com>

On 12/25/09, Marco Atzeri > wrote:
> --- Gio 24/12/09, Linh Phan  ha scritto:
>
>> Hi,
>>
>>   I have wondering why when I do floating-point math
>> in Cygwin 1.7, Cygwin does not output the "correct" result,
>> eg:
>>
>> main () {
>>   double qcb[4] = {0.41585680180318363,
>> 0.70823637211274604, 0.089200955545759006,
>> -0.56347808399291521};
>>   double qba[4] = {0.09648538897087118,
>> -0.37032037362514358,
>> 0.89403246842889339,   0.23293633462432001};
>>   double qca[4];
>>
>>   qca[0]= qcb[3]*qba[0] + qcb[2]*qba[1] -
>> qcb[1]*qba[2] + qcb[0]*qba[3];
>>   qca[1]= -qcb[2]*qba[0] + qcb[3]*qba[1] +
>> qcb[0]*qba[2] + qcb[1]*qba[3];
>>   qca[2]= qcb[1]*qba[0] - qcb[0]*qba[1] +
>> qcb[3]*qba[2] + qcb[2]*qba[3];
>>   qca[3]= -qcb[0]*qba[0] - qcb[1]*qba[1] -
>> qcb[2]*qba[2] + qcb[3]*qba[3];
>>
>>   // I am expecting the result to be:
>>
>>   printf("qca = [-0.623718486146499718,
>> 0.736824293298044886, -0.260654850643024127,
>> 0.011147182658310384] CORRECT (LINUX/SOLARIS)\n");
>>
>>   // I get the above result on both Linux/Solaris AND
>> even when running under Cygwin's gdb and using gdb's printf
>>
>>   // but when I compile and run on cygwin, I get the
>> wrong result:
>
> s/wrong/different
>
>>   printf("qca = [%.18f, %.18f, %.18f, %.18f] NOT
>> CORRECT (CYGWIN)\n",qca[0],qca[1],qca[2],qca[3]);
>> }
>>
>> Here is the output when I compiled and run on cygwin:
>>
>> qca = [-0.623718486146499718, 0.736824293298044886,
>> -0.260654850643024127, 0.011147182658310384] CORRECT
>> (LINUX/SOLARIS)
>> qca = [-0.623718486146499607, 0.736824293298044886,
>> -0.260654850643024072, 0.011147182658310373] NOT CORRECT
>> (CYGWIN)
>>
>> Here's the output when running under cygwin's gdb (it
>> matches Linux/Solaris outputs but not cygwin’s output):
>>
>> (gdb) printf "%.18f\n", qcb[3]*qba[0] + qcb[2]*qba[1] -
>> qcb[1]*qba[2] + qcb[0]*qba[3]
>> -0.623718486146499718
>> (gdb) printf "%.18f\n", -qcb[2]*qba[0] + qcb[3]*qba[1] +
>> qcb[0]*qba[2] + qcb[1]*qba[3]
>> 0.736824293298044886
>> (gdb) printf "%.18f\n", qcb[1]*qba[0] - qcb[0]*qba[1] +
>> qcb[3]*qba[2] + qcb[2]*qba[3]
>> -0.260654850643024127
>> (gdb) printf "%.18f\n", -qcb[0]*qba[0] - qcb[1]*qba[1] -
>> qcb[2]*qba[2] + qcb[3]*qba[3]
>> 0.011147182658310384
>> (gdb)
>>
>> These errors accumulate and eventually builds up so much
>> that the my simulation result is totally wrong when running
>> under cygwin.

What exactly is the magnitude of the error? Sorry, it is hard to pick
out from here. Generally you design numerical code such that it
knows or discovers various limitations.

>>
>> Please help.
>
> I will bet
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323

That should have never been accepted as a bug because it
relies on a floating point equality test, namely
if (y != y2) printf("error\n");
Whlie this test may make FP look stupid, consider the case where
one of the constants is expressed as 1/3 or something known
to be hard to get exact.

Close is good enough in horseshoes, handgrenades, and floating point.
As the bug report suggests, look at the books on IA-32 FPU at intel. MSFT
IIRC even has options for "improving floating point consistency" but you
don't want that if you care about performance and robust floating point.
I've seen filter design code that makes various comments about compilers
that believe that fp associates and commutes- anything that is only marginally
stable needs to be aware of numerical issues. I've also run into this
problem with "arithmetic encoders" but if you really want repro
performance you are stuck with something like java.

>
> see
> http://hal.archives-ouvertes.fr/hal-00128124
> for further detailed info on rounding
>
> if so it is not a bug is a feature
> :-)
>
> by the way on (cygwin) octave that round as you expect I have
>  printf("%.18g\n",qca)
> -0.623718486146499718
> 0.736824293298044886
> -0.260654850643024127
> 0.0111471826583103839
>
>
>
>> Thanks,
>>
>> Linh
>>
>
> Marco
>
>
>
>
>
> --
> Problem reports:       http://cygwin.com/problems.html
> FAQ:                   http://cygwin.com/faq/
> Documentation:         http://cygwin.com/docs.html
> Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple
>
>

-- 
[email-masked]
note new address 2009-12-16:
Mike Marchywka
1975 Village Round
Marietta GA 30064
415-264-8477 (w)<- use this
404-788-1216 (C)<- leave message
989-348-4796 (P)<- emergency only
[email-masked]
Note: If I am asking for free stuff, I normally use for hobby/non-profit
information but may use in investment forums, public and private.
Please indicate any concerns if applicable.
Note: hotmail is censoring incoming mail using random criteria beyond
my control and often hangs my browser
but all my subscriptions are here..., try also [email-masked]

--
Problem reports:       http://cygwin.com/problems.html
FAQ:                   http://cygwin.com/faq/
Documentation:         http://cygwin.com/docs.html
Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple