caml-list - the Caml user's mailing list
 help / color / mirror / Atom feed
* [Caml-list] Complex Numbers
@ 2001-03-28  2:19 David McClain
  2001-03-28  3:38 ` [Caml-list] OCaml, where art thou? Anjayan Puvananathan
  0 siblings, 1 reply; 13+ messages in thread
From: David McClain @ 2001-03-28  2:19 UTC (permalink / raw)
  To: caml-list

"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


^ permalink raw reply	[flat|nested] 13+ messages in thread

end of thread, other threads:[~2001-04-03 18:06 UTC | newest]

Thread overview: 13+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2001-03-28  2:19 [Caml-list] Complex Numbers David McClain
2001-03-28  3:38 ` [Caml-list] OCaml, where art thou? Anjayan Puvananathan
2001-03-30 14:51   ` Xavier Leroy
2001-03-30 23:12     ` David Fox
2001-03-31  9:31       ` [Caml-list] Using HTML as a standard GUI for Ocaml Mattias Waldau
2001-04-01 12:52         ` Fergus Henderson
2001-04-01 20:11         ` Sven LUTHER
2001-04-01 20:35           ` Bruce Hoult
2001-04-02 10:09             ` Sven LUTHER
2001-04-02 15:53               ` CREGUT Pierre FTRD/DTL/LAN
2001-04-02 19:16         ` Gerd Stolpmann
2001-04-03 18:06           ` Mattias Waldau
2001-04-02 13:21       ` [Caml-list] Threads in OCaml Xavier Leroy

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).