1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
(*
* Copyright (c) 1997-1999, 2003 Massachusetts Institute of Technology
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*)
(* $Id: complex.ml,v 1.23 2003/03/16 23:43:46 stevenj Exp $ *)
(* abstraction layer for complex operations *)
(* type of complex expressions *)
open Exprdag
open Exprdag.LittleSimplifier
type expr = CE of node * node
let two = CE (makeNum Number.two, makeNum Number.zero)
let one = CE (makeNum Number.one, makeNum Number.zero)
let zero = CE (makeNum Number.zero, makeNum Number.zero)
let inverse_int n = CE (makeNum (Number.div Number.one
(Number.of_int n)),
makeNum Number.zero)
let times_4_2 (CE (a, b)) (CE (c, d)) =
CE (makePlus [makeTimes (a, c); makeUminus (makeTimes (b, d))],
makePlus [makeTimes (a, d); makeTimes (b, c)])
let simple = function
Num a -> Number.is_zero a or Number.is_one a or Number.is_mone a
| _ -> false
let rec times_3_3 (CE (a, b)) (CE (c, d)) =
(* refuse to do the 3-3 algorithm if a=1, i, -i, -1, etc. *)
if simple a or simple b or simple c or simple d then
times_4_2 (CE (c, d)) (CE (a, b))
else match a with
Num _ ->
let amb = makePlus [a; makeUminus b]
and cpd = makePlus [c; d]
and apb = makePlus [a; b]
in let apbc = makeTimes (apb, c)
and bcpd = makeTimes (b, cpd)
and ambd = makeTimes (amb, d)
in CE (makePlus [apbc; makeUminus bcpd],
makePlus [bcpd; ambd])
| _ -> match c with
Num _ -> times_3_3 (CE (c, d)) (CE (a, b))
| _ -> times_4_2 (CE (a, b)) (CE (c, d))
let times a b =
if !Magic.times_3_3 then
times_3_3 a b
else
times_4_2 a b
let uminus (CE (a, b)) = CE (makeUminus a, makeUminus b)
(* hack to swap real<->imaginary. Used by hc2hc codelets *)
let swap_re_im (CE (r, i)) = CE (i, r)
(* complex exponential (of root of unity); returns exp(2*pi*i/n * m) *)
let exp n i =
let (c, s) = Number.cexp n i
in CE (makeNum c, makeNum s)
(* complex sum *)
let plus a =
let rec unzip_complex = function
[] -> ([], [])
| ((CE (a, b)) :: s) ->
let (r,i) = unzip_complex s
in
(a::r), (b::i) in
let (c, d) = unzip_complex a in
CE (makePlus c, makePlus d)
(* extract real/imaginary *)
let real (CE (a, b)) = CE (a, makeNum Number.zero)
let imag (CE (a, b)) = CE (makeNum Number.zero, b)
let conj (CE (a, b)) = CE (a, makeUminus b)
let abs_sqr (CE (a, b)) = makePlus [makeTimes (a, a);
makeTimes (b, b)]
(*
* special cases for complex numbers w where |w| = 1
*)
(* (a + bi)^2 = (2a^2 - 1) + 2abi *)
let wsquare (CE (a, b)) =
let twoa = makeTimes (makeNum Number.two, a)
in let twoasq = makeTimes (twoa, a)
and twoab = makeTimes (twoa, b) in
CE (makePlus [twoasq; makeUminus (makeNum Number.one)], twoab)
(*
* compute w^n given w^{n-1}, w^{n-2}, and w, using the identity
*
* w^n + w^{n-2} = w^{n-1} (w + w^{-1}) = 2 w^{n-1} Re(w)
*)
let wthree (CE (an1, bn1)) wn2 (CE (a, b)) =
let twoa = makeTimes (makeNum Number.two, a)
in let twoa_wn1 = CE (makeTimes (twoa, an1),
makeTimes (twoa, bn1))
in plus [twoa_wn1; (uminus wn2)]
(* abstraction of sum_{i=0}^{n-1} *)
(* let sigma a b f = plus (Util.forall :: a b f) *)
let sigma a b f =
let rec loop a =
if (a >= b) then []
else (f a) :: (loop (a + 1))
in plus (loop a)
(* complex variables *)
type variable = CV of Variable.variable * Variable.variable
let load_var (CV (vr, vi)) =
CE (Load vr, Load vi)
let store_var (CV (vr, vi)) (CE (xr, xi)) =
[Store (vr, xr); Store (vi, xi)]
let store_real (CV (vr, vi)) (CE (xr, xi)) =
[Store (vr, xr)]
let store_imag (CV (vr, vi)) (CE (xr, xi)) =
[Store (vi, xi)]
let access what k =
let (r, i) = what k
in CV (r, i)
let access_input = access Variable.access_input
let access_output = access Variable.access_output
let access_twiddle = access Variable.access_twiddle
|
