From mboxrd@z Thu Jan 1 00:00:00 1970 Received: (from majordomo@localhost) by pauillac.inria.fr (8.7.6/8.7.3) id LAA27397; Wed, 28 Mar 2001 11:20:41 +0200 (MET DST) X-Authentication-Warning: pauillac.inria.fr: majordomo set sender to owner-caml-list@pauillac.inria.fr using -f Received: from concorde.inria.fr (concorde.inria.fr [192.93.2.39]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id LAA27393 for ; Wed, 28 Mar 2001 11:20:40 +0200 (MET DST) Received: from smtp6.mindspring.com (smtp6.mindspring.com [207.69.200.110]) by concorde.inria.fr (8.11.1/8.10.0) with ESMTP id f2S9KcT04220 for ; Wed, 28 Mar 2001 11:20:38 +0200 (MET DST) Received: from dylan (1Cust105.tnt4.tucson.az.da.uu.net [63.11.145.105]) by smtp6.mindspring.com (8.9.3/8.8.5) with SMTP id EAA26193; Wed, 28 Mar 2001 04:20:33 -0500 (EST) Message-ID: <000901c0b768$64b9bdb0$210148bf@dylan> From: "David McClain" To: , References: <200103280844.KAA29428@ilt.fhg.de> Subject: Re: [Caml-list] Re: Complex Arithmetic Date: Wed, 28 Mar 2001 02:20:53 -0700 MIME-Version: 1.0 Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: 7bit X-Priority: 3 X-MSMail-Priority: Normal X-Mailer: Microsoft Outlook Express 5.00.2919.6600 X-MimeOLE: Produced By Microsoft MimeOLE V5.00.2919.6600 Sender: owner-caml-list@pauillac.inria.fr Precedence: bulk > as far as I understood Kahan says that the complex arithmetic in > FORTRAN and many C++ libraries is wrong. I tested: > > sqrt( (0.0 , -1.0)**2 ) vs. sqrt( (0.0 , 1.0)**2 ) > > on a DEC Alpha with FORTRAN and C++ (source and output below). > They both give the correct result. A naive approach abviously gives > the same result (0.0, 1.0) for both expressions. Can you tell us how you > implemented the correct behaviour in OCaml (or Nml)? Well, well... it appears that DEC Alpha has taken Kahan to heart! I have never used an Alpha, and perhaps I should. I understood that their implementation of FP traps is a bit strange, using trap gates interspersed with the code so you don't actually get the trap at the location of the error unless you pepper your code to death with trap gates... But I have heard many great things about the speed of these boxes! > BTW what about sqrt( (-1.0 , -1.0)**2 ) and sqrt( (-1.0 , 1.0)**2 ), which > result in (1.0, 1.0) and (1.0, -1.0) on the DEC Alpha. What should the > results be in this case? These cases again illustrate Kahan's principal about the limitations of the language of [rectangular] pairs... (-1,-1) when squared exceeds the principal Riemann sheet and so all attempts at representing complex arithmetic with rectilinear coordinates of the principal sheet are doomed to failure. Same for (-1,1). So in these cases, even NML will fail like all the others. Only polar representation which can display arbitrary levels of the Riemann surface could possibly succeed in such cases. My NML implementation [correct for the principal sheet arithmetic only] is available to anyone, and I simply (ha!) made sure that all the algebraic relations between 0+ and 0- were properly observed. For example, IIRC, x + 0+ != x (this fails in the case where x = 0-. It succeeds everywhere else). x + 0- == x x * 0- = 0- unless x = 0- in which case 0- * 0- = 0+ (sheesh!!) (of course I may be mistaken here in my mental recall. I have to reexamine my code to be certain) When implementing NML, I actually wrote a simple-minded algebraic theorem prover in OCaml to assist me in this process. That's another terrific thing about OCaml - you can perform symbolic calculations very easily, perhaps even more easily than in Lisp (although I also love Lisp). The whole prover took less than a hundred lines of OCaml or so... and it ran correctly the first time out! - DM ------------------- To unsubscribe, mail caml-list-request@inria.fr. Archives: http://caml.inria.fr