[Camllist] Re: Complex Arithmetic

wester@i...
 David McClain
[
Home
]
[ Index:
by date

by threads
]
[ Message by date: previous  next ] [ Message in thread: previous  next ] [ Thread: previous  next ]
[ Message by date: previous  next ] [ Message in thread: previous  next ] [ Thread: previous  next ]
Date:   (:) 
From:  David McClain <dmcclain1@m...> 
Subject:  Re: [Camllist] Re: Complex Arithmetic 
> 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 simpleminded 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 camllistrequest@inria.fr. Archives: http://caml.inria.fr