anil.recoil.org
My working unpac repository
1(**************************************************************************)
2(* *)
3(* OCaml *)
4(* *)
5(* Xavier Leroy, projet Cristal, INRIA Rocquencourt *)
6(* *)
7(* Copyright 2002 Institut National de Recherche en Informatique et *)
8(* en Automatique. *)
9(* *)
10(* All rights reserved. This file is distributed under the terms of *)
11(* the GNU Lesser General Public License version 2.1, with the *)
12(* special exception on linking described in the file LICENSE. *)
13(* *)
14(**************************************************************************)
15
16(* Complex numbers *)
17
18type t = { re: float; im: float }
19
20let zero = { re = 0.0; im = 0.0 }
21let one = { re = 1.0; im = 0.0 }
22let i = { re = 0.0; im = 1.0 }
23
24let add x y = { re = x.re +. y.re; im = x.im +. y.im }
25
26let sub x y = { re = x.re -. y.re; im = x.im -. y.im }
27
28let neg x = { re = -. x.re; im = -. x.im }
29
30let conj x = { re = x.re; im = -. x.im }
31
32let mul x y = { re = x.re *. y.re -. x.im *. y.im;
33 im = x.re *. y.im +. x.im *. y.re }
34
35let div x y =
36 if abs_float y.re >= abs_float y.im then
37 let r = y.im /. y.re in
38 let d = y.re +. r *. y.im in
39 { re = (x.re +. r *. x.im) /. d;
40 im = (x.im -. r *. x.re) /. d }
41 else
42 let r = y.re /. y.im in
43 let d = y.im +. r *. y.re in
44 { re = (r *. x.re +. x.im) /. d;
45 im = (r *. x.im -. x.re) /. d }
46
47let inv x = div one x
48
49let norm2 x = x.re *. x.re +. x.im *. x.im
50
51let norm x = Float.hypot x.re x.im
52
53let arg x = atan2 x.im x.re
54
55let polar n a = { re = cos a *. n; im = sin a *. n }
56
57let sqrt x =
58 if x.re = 0.0 && x.im = 0.0 then { re = 0.0; im = 0.0 }
59 else begin
60 let r = abs_float x.re and i = abs_float x.im in
61 let w =
62 if r >= i then begin
63 let q = i /. r in
64 sqrt(r) *. sqrt(0.5 *. (1.0 +. sqrt(1.0 +. q *. q)))
65 end else begin
66 let q = r /. i in
67 sqrt(i) *. sqrt(0.5 *. (q +. sqrt(1.0 +. q *. q)))
68 end in
69 if x.re >= 0.0
70 then { re = w; im = 0.5 *. x.im /. w }
71 else { re = 0.5 *. i /. w; im = if x.im >= 0.0 then w else -. w }
72 end
73
74let exp x =
75 let e = exp x.re in { re = e *. cos x.im; im = e *. sin x.im }
76
77let log x = { re = log (norm x); im = atan2 x.im x.re }
78
79let pow x y = exp (mul y (log x))