[Caml-list] Complex Numbers

["David McClain" <dmcclain1@mindspring.com>]

"One of the articles by William Kahan (www.cs.berkeley.edu/~wkahan/)
"How JAVA's Floating-Point Hurts Everyone Everywhere" discusses
(among other things) some traps in complex number implementations,
including those in Fortran.
"

My NML, written in OCaml, paid particular attention to these issues after I
read everything I could from Kahan. To my knowledge, aside from Common Lisp
(perhaps), it is the only proper implementation of complex arithmetic
available.

In particular, Kahan's demo of "Borda's Pipe", shows that stream lines are
incorrectly computed in every language I have tried: Mathematica, RSI/IDL,
Fortran, C/C++, Basic, with the possible exception of Common Lisp. Only NML
has passed demonstrably passed this test.

The difficulty stems from the fact that, to paraphrase Kahan, "a tuple
representation of complex numbers is insufficient, given the Riemann
surfaces of common transcendental functions". There are two zeros defined in
the IEEE Floating Point Standard: 0+, and 0-.  Arithmetic in the floating
point domain do not obey classical algrebraic relations with respect to
conventional identity elements.

I sweated out two solid weeks in gaining a full understanding of his rather
cryptic statement. NML represents my final victory over this issue and I can
thank OCaml and the INRIA team for their wonderful implementation language
that made it all so readily possible.

- DM


Borda NML code follows:

--------------------------------------------------
(* borda.nml -- test of complex arithmetic in NML *)
(* DM/MCFA  09/99 *)

let g z = z^2 + z*sqrt(z^2+1)
let f z = 1 + g(z) + log(g z)

let ls =
  let ptor r theta =
    complex (r*cos theta) (r*sin theta)
  in
  let t = tan(pi/2*[0, 0.001, .. 1]) in
    map (fun ang ->
    ptor t ang)
      (lis (pi * [-0.5, -0.45, .. 0.5]))

let rec pshow (l, ... as args) =
  let clip x =
    x #< 400 #> -400
  in
  let z = f(l) in
    oplot(clip(re z), clip(im z) | args @ [clip: true])

let borda() =
  let s = 5 in
    axes(xrange: [-s,s], yrange: [-s,s],
 title: "The Correct Borda!");
    pshow(hd ls, color: red, thick: 2);
    app pshow (tl ls)

let _ = borda()



-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr