X-Recipient: archive-cygwin AT delorie DOT com
X-SWARE-Spam-Status: No, hits=-1.9 required=5.0 	tests=AWL,BAYES_00,SARE_MSGID_LONG40,SPF_PASS
X-Spam-Check-By: sourceware.org
MIME-Version: 1.0
In-Reply-To: <724619.16892.qm@web25501.mail.ukl.yahoo.com>
References: <200912242119 DOT nBOLJCr7019530 AT aludra DOT usc DOT edu> 	 <724619 DOT 16892 DOT qm AT web25501 DOT mail DOT ukl DOT yahoo DOT com>
Date: Fri, 25 Dec 2009 07:47:53 -0500
Message-ID: <b0ddfae80912250447u1d30ff42wbb5263727b424eb3@mail.gmail.com>
Subject: Re: floating-point math problem in Cygwin 1.7?
From: mike marchywka <marchywka AT gmail DOT com>
To: cygwin AT cygwin DOT com
Cc: Linh Phan <linhphan AT usc DOT edu>
Content-Type: text/plain; charset=windows-1252
Content-Transfer-Encoding: quoted-printable
X-IsSubscribed: yes
Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm
Precedence: bulk
List-Id: <cygwin.cygwin.com>
List-Unsubscribe: <mailto:cygwin-unsubscribe-archive-cygwin=delorie DOT com AT cygwin DOT com>
List-Subscribe: <mailto:cygwin-subscribe AT cygwin DOT com>
List-Archive: <http://sourceware.org/ml/cygwin/>
List-Post: <mailto:cygwin AT cygwin DOT com>
List-Help: <mailto:cygwin-help AT cygwin DOT com>, <http://sourceware.org/ml/#faqs>
Sender: cygwin-owner AT cygwin DOT com
Mail-Followup-To: cygwin AT cygwin DOT com
Delivered-To: mailing list cygwin AT cygwin DOT 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] =3D {0.41585680180318363,
>> 0.70823637211274604, 0.089200955545759006,
>> -0.56347808399291521};
>>   double qba[4] =3D {0.09648538897087118,
>> -0.37032037362514358,
>> 0.89403246842889339,   0.23293633462432001};
>>   double qca[4];
>>
>>   qca[0]=3D qcb[3]*qba[0] + qcb[2]*qba[1] -
>> qcb[1]*qba[2] + qcb[0]*qba[3];
>>   qca[1]=3D -qcb[2]*qba[0] + qcb[3]*qba[1] +
>> qcb[0]*qba[2] + qcb[1]*qba[3];
>>   qca[2]=3D qcb[1]*qba[0] - qcb[0]*qba[1] +
>> qcb[3]*qba[2] + qcb[2]*qba[3];
>>   qca[3]=3D -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 =3D [-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 =3D [%.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 =3D [-0.623718486146499718, 0.736824293298044886,
>> -0.260654850643024127, 0.011147182658310384] CORRECT
>> (LINUX/SOLARIS)
>> qca =3D [-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=92s 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=3D323

That should have never been accepted as a bug because it
relies on a floating point equality test, namely
if (y !=3D 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 margina=
lly
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
>
>


--=20
marchywka AT gmail DOT com
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
marchywka AT hotmail DOT com
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 marchywka AT yahoo DOT com

--
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