[Caml-list] Complex numbers in ocaml

[Caml-list] Complex numbers in ocaml

[Post Office!~/sentmail]

I'm trying to write a system for complex numbers in ocaml for some
numerical work I'm doing, and I thought I'd ask if anyone else has such
code already existant.

Otherwise, given an abstract for numbers
	type number = Int of int | Float of float | Complex of float *
float | Error;;

is there an easy way to coerce numbers back out of type number in order to
operate on them?

Fred Ross

-------------------
Bug reports: http://caml.inria.fr/bin/caml-bugs  FAQ: http://caml.inria.fr/FAQ/
To unsubscribe, mail caml-list-request@inria.fr  Archives: http://caml.inria.fr

Re: [Caml-list] Complex numbers in ocaml

[John Max Skaller]

"Post Office!~/sentmail" wrote:
> 
> I'm trying to write a system for complex numbers in ocaml for some
> numerical work I'm doing, and I thought I'd ask if anyone else has such
> code already existant.
> 
> Otherwise, given an abstract for numbers
>         type number = Int of int | Float of float | Complex of float *
> float | Error;;
> 
> is there an easy way to coerce numbers back out of type number in order to
> operate on them?

	match n with 
	| Int i -> ..
	| Float x -> ...
	| Complex (x,y) -> ...
	| Error -> ...

-- 
John (Max) Skaller, mailto:skaller@maxtal.com.au 
10/1 Toxteth Rd Glebe NSW 2037 Australia voice: 61-2-9660-0850
New generation programming language Felix  http://felix.sourceforge.net
Literate Programming tool Interscript     
http://Interscript.sourceforge.net
-------------------
Bug reports: http://caml.inria.fr/bin/caml-bugs  FAQ: http://caml.inria.fr/FAQ/
To unsubscribe, mail caml-list-request@inria.fr  Archives: http://caml.inria.fr

Re: [Caml-list] Complex numbers in OCaml

[Arturo Borquez]

On Tue, 27 March 2001, wester@ilt.fhg.de wrote:

> 
> Hi,
> 
> I need complex numbers. I tried it this way:
> 
> type number  = Real of float | Complex of (float*float);;
> 
> let add_num a b = 
>   match (a,b) with
>       (Real a, Real b)          -> Real (a +. b)
>     | (Real a, Complex (bx,by)) -> Complex (a +. bx, by) .................

Hi: Rolf:
Here is a little example of an alternative to operate with complex numbers. The key is to define a complex a a duple and define new opertators to do computations, ie +! for addition -! for substraction *! /!...

let (+!) (xr, xi) (yr, yi) = (xr +. yr), (xi +. yi);;
let a = (1.0, 2.0);;
let b = (0.0, -4.75);;
let z = a +! b;;
z;;
val z : float * float = 1, -2.75

The same applies for sub, mpy and div ops.
hope this helps.

Best regards.
Arturo Borquez


Find the best deals on the web at AltaVista Shopping!
http://www.shopping.altavista.com
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr

[Caml-list] Complex numbers in OCaml

[wester]

Hi,

I need complex numbers. I tried it this way:

type number  = Real of float | Complex of (float*float);;

let add_num a b = 
  match (a,b) with
	  (Real a, Real b)          -> Real (a +. b)
	| (Real a, Complex (bx,by)) -> Complex (a +. bx, by)
	| (Complex (ax,ay), Real b) -> Complex (ax +. b, ay)
	| (Complex (ax,ay), Complex (bx,by)) -> Complex (ax +. bx, ay +. by);;
let sub_num a b = 
  match (a,b) with
	  (Real a, Real b)          -> Real (a -. b)
	| (Real a, Complex (bx,by)) -> Complex (a -. bx, by)
	| (Complex (ax,ay), Real b) -> Complex (ax -. b, ay)
	| (Complex (ax,ay), Complex (bx,by)) -> Complex (ax -. bx, ay -. by);;
let mul_num a b = 
  match (a,b) with
	  (Real a, Real b)          -> Real (a *. b)
	| (Real a, Complex (bx,by)) -> Complex (a *. bx, a *. by)
	| (Complex (ax,ay), Real b) -> Complex (ax *. b, ay)
	| (Complex (ax,ay), Complex (bx,by)) -> Complex (ax *. bx -. ay *. by, ax *. by +. ay *. bx);;
let div_num a b = 
  match (a,b) with
	  (Real a, Real b)          -> Real (a /. b)
	| (Real a, Complex (bx,by)) -> Complex (a *. bx /. (bx**2.0 +. by**2.0), 
                                            -. a *. by /. (bx**2.0 +. by**2.0))
	| (Complex (ax,ay), Real b) -> Complex (ax /. b, ay /. b)
	| (Complex (ax,ay), Complex (bx,by)) -> let n = bx**2.0 +. by**2.0 in
                                            Complex ((ax *. bx +. ay *. by)/.n, 
                                                     (ay *. bx -. ax *. by)/.n);;
let conj a = 
  match a with
	  Real a -> Real a;
	| Complex (ax,ay) -> Complex (ax, -. ay);;
let real a = 
  match a with
	  Real a -> a
	|Complex (ax, ay) -> ax;;
let imag a = 
  match a with
	  Real a -> 0.0
	|Complex (ax, ay) -> ay;;

div_num (Complex (2.0,3.0)) (Complex (3.0,2.0));;

This works but is quite cumbersome to use. Is there any other way to do it and how can
arrays of complex numbers be implemented efficiently in OCaml? 

Thanks in advance.

Rolf Wester

 

-------------------------------------
Rolf Wester
wester@ilt.fhg.de
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr

Re: [Caml-list] Complex numbers in OCaml

[Markus Mottl]

wester@ilt.fhg.de schrieb am Tuesday, den 27. March 2001:
> This works but is quite cumbersome to use. Is there any other way to
> do it and how can arrays of complex numbers be implemented efficiently
> in OCaml?

It would be indeed nice to see Bigarrays support complex numbers so that
we can interface to more functions in Fortran-libraries...

Regards,
Markus Mottl

-- 
Markus Mottl, mottl@miss.wu-wien.ac.at, http://miss.wu-wien.ac.at/~mottl
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr

Re: [Caml-list] Complex numbers in OCaml

[Jan Skibinski]

On Tue, 27 Mar 2001, Markus Mottl wrote:

> wester@ilt.fhg.de schrieb am Tuesday, den 27. March 2001:
> > This works but is quite cumbersome to use. Is there any other way to
> > do it and how can arrays of complex numbers be implemented efficiently
> > in OCaml?
> 
> It would be indeed nice to see Bigarrays support complex numbers so that
> we can interface to more functions in Fortran-libraries...
> 
	
Or not. :-)

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.

Jan
 
	

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

Re: [Caml-list] Complex numbers in OCaml

[Markus Mottl]

On Tue, 27 Mar 2001, Jan Skibinski wrote:
> Or not. :-)
> 
> 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.

Really interesting article! I always believed that doing floating-point
arithmetic correctly was a tricky business - now I am convinced... ;)

Good to know for us that OCaml is at least extremely close to the IEEE
754 standard, though there is still criticism in this article that is
relevant for OCaml, too (e.g. things like non-availability of extended
double-precision; lack of operator overloading; etc.). Luckily, my needs
for numerical computation are not this high...

Regards,
Markus Mottl

-- 
Markus Mottl, mottl@miss.wu-wien.ac.at, http://miss.wu-wien.ac.at/~mottl
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr

Re: [Caml-list] Complex numbers in OCaml

[Michael Hohn]

>> ...
>> 
>> 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

Re: [Caml-list] Complex numbers in OCaml

[Xavier Leroy]

> I need complex numbers. I tried it this way:
> type number  = Real of float | Complex of (float*float);;
> [...]
> This works but is quite cumbersome to use.
> Is there any other way to do it

Unless you have strong reasons for doing so, I'd advise against the
mixed-mode arithmetic (i.e. having a special case for Real): it makes
operations slower and more complex.

In terms of performance, the best representation is a record type:

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

since the compiler will then unbox the two floats.  This would not be
the case for a pair of floats:

        type complex = float * float

The concrete syntax is a bit heavy, e.g. {re = 1.0; im = 0.0}, but at
least highly non-ambiguous :-)  If you're desperate for a more compact
notation, Camlp4 could help.

One last performance hint: write b *. b rather than b ** 2.0,
it's much, much faster.

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

Just regular arrays of complex numbers, e.g. the type "complex array",
should be reasonably efficient, although you pay one extra indirection
compared with the equivalent C or Fortran array.

As Markus mentioned, I'm considering extending the Bigarray module to
handle complex numbers as well.  This would help interfacing with
Fortran, and provide essentially the same data representation than in
C and Fortran.

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