Re: [Caml-list] Complex numbers in OCaml

[Michael Hohn <hohn@math.utah.edu>]

>> ...
>> 
>> I need complex numbers. I tried it this way:
>> 
>> ...

I also needed complex numbers a while ago.
There are several options:
o   Make a FORTRAN binding.
    This is tedious, and for simple low-level operations (add, etc.)
    will probably not help performance.  If 

o   Implement them in ocaml.
    This is very simple when using user-defined operators.
    Also, it is fast enough for all simple applications


>> ...
>> This works but is quite cumbersome to use. Is there any other way to do it
>> ...

Below is the module I used, based on a struct-like product type instead
of a sum type.

>> ...
>> and how can arrays of complex numbers be implemented efficiently in OCaml? 
>> ...

Arrays of  (C) doubles are optimized by ocamlopt, but for complex
types, which are just structs, I do not know.

For the ultimate in efficiency of complex arrays, there is always
FORTRASH :(

Cheers,
Michael


(* $Id: cmplx.ml,v 1.6.2.1 2000/04/05 16:29:57 hohn Exp $
   Elementary complex number support functions.
   Taken from libf2c/libF77.
*)

type cmplx = { re : float; im : float}

let cx re im = { re = re ; im = im}

let zarg z = atan2 z.im z.re 

let conj z = { re = z.re ; im = -. z.im}

let (+:) a b = {re =  (a.re +. b.re); im =  (a.im +. b.im) }

let (-:) a b = {re = (a.re -. b.re); im = (a.im -. b.im) }

let ( *: ) a b =
    {
      re = (  a.re *. b.re -. a.im *. b.im ) ; 
      im = (a.re *. b.im +. a.im *. b.re )
    } 

let f__cabs (real : float) (imag : float) = 
    let real = if (real < 0.0) then (-. real) else real in
    let imag = if (imag < 0.0) then (-. imag) else imag in
    let (real, imag) = if(imag > real) then (imag, real) else (real, imag) in
    if ((real +. imag) = real) then real
    else (
    let temp = imag /.real in
    real *. (sqrt (1.0 +. temp*.temp))
       )

let zabs z = {re = f__cabs z.re z.im ; im = 0.0}

let dabs x = abs_float x 

let ( /: ) a b = 
    let abr = if( b.re < 0.0) then (-. b.re) else (b.re)  in 
    let abi = if( b.im < 0.0) then (-. b.im) else b.im in 
    if ( abr <= abi ) then
    (
         if (abi = 0.0) then
             raise Division_by_zero
	      else
             let ratio = b.re /. b.im in
             let den = b.im *. (1.0 +. ratio*.ratio) in 
             {
	            re = (a.re *. ratio +. a.im) /. den ; 
               im = (a.im *. ratio -. a.re) /. den 
	            } 
        )
    else
        (
         let ratio = b.im /. b.re in 
         let den = b.re *. (1.0 +. ratio*.ratio) in
         {
	    re = (a.re +. a.im *. ratio) /. den ;
           im = (a.im -. a.re *. ratio) /. den
	    } 
        )

	    
let zexp z =
    let zi = z.im in 
    let expx = exp z.re in 
    {
      re = expx *. cos(zi);
      im = expx *. sin(zi)
    } 

let zlog z = 
    let zi = z.im and zr = z.re in 
    {
      im = atan2 zi zr;
      re = (log ( f__cabs zr zi ) )
    } 

(* Integer power z**i, for small i *)
let zipow z i =
    let rv z = 
    if (i < 0) then ((cx 1.0 0.0) /: z)
    else z in
    let i = if i < 0 then -i else i in
    match i with
    | 0 -> cx 1.0 0.0
    | 1 -> rv z
    | 2 -> rv (z *: z)
    | 3 -> rv (z *: z *: z)
    | 4 -> rv (z *: z *: z *: z)
    | _ ->
    let p = ref z in 
    for k=2 to i do p := !p *: z done ;
    rv !p
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr