Browse thread
[Caml-list] Complex numbers in OCaml
-
wester@i...
- Markus Mottl
- Michael Hohn
- Xavier Leroy
[
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: | Michael Hohn <hohn@m...> |
| Subject: | Re: [Caml-list] Complex numbers in OCaml |
>> ...
>>
>> 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