cannam@95: (*
cannam@95:  * Copyright (c) 1997-1999 Massachusetts Institute of Technology
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  *)
cannam@95: 
cannam@95: (* abstraction layer for complex operations *)
cannam@95: open Littlesimp
cannam@95: open Expr
cannam@95: 
cannam@95: (* type of complex expressions *)
cannam@95: type expr = CE of Expr.expr * Expr.expr
cannam@95: 
cannam@95: let two = CE (makeNum Number.two, makeNum Number.zero)
cannam@95: let one = CE (makeNum Number.one, makeNum Number.zero)
cannam@95: let i = CE (makeNum Number.zero, makeNum Number.one)
cannam@95: let zero = CE (makeNum Number.zero, makeNum Number.zero)
cannam@95: let make (r, i) = CE (r, i)
cannam@95: 
cannam@95: let uminus (CE (a, b)) =  CE (makeUminus a, makeUminus b)
cannam@95: 
cannam@95: let inverse_int n = CE (makeNum (Number.div Number.one (Number.of_int n)),
cannam@95: 			makeNum Number.zero)
cannam@95: 
cannam@95: let inverse_int_sqrt n = 
cannam@95:   CE (makeNum (Number.div Number.one (Number.sqrt (Number.of_int n))),
cannam@95:       makeNum Number.zero)
cannam@95: let int_sqrt n = 
cannam@95:   CE (makeNum (Number.sqrt (Number.of_int n)),
cannam@95:       makeNum Number.zero)
cannam@95: 
cannam@95: let nan x = CE (NaN x, makeNum Number.zero)
cannam@95: 
cannam@95: let half = inverse_int 2
cannam@95: 
cannam@95: let times3x3 (CE (a, b)) (CE (c, d)) = 
cannam@95:   CE (makePlus [makeTimes (c, makePlus [a; makeUminus (b)]);
cannam@95: 	        makeTimes (b, makePlus [c; makeUminus (d)])],
cannam@95:       makePlus [makeTimes (a, makePlus [c; d]);
cannam@95: 	        makeUminus(makeTimes (c, makePlus [a; makeUminus (b)]))])
cannam@95: 
cannam@95: let times (CE (a, b)) (CE (c, d)) = 
cannam@95:   if not !Magic.threemult then
cannam@95:     CE (makePlus [makeTimes (a, c); makeUminus (makeTimes (b, d))],
cannam@95:         makePlus [makeTimes (a, d); makeTimes (b, c)])
cannam@95:   else if is_constant c && is_constant d then
cannam@95:     times3x3 (CE (a, b)) (CE (c, d))
cannam@95:   else (* hope a and b are constant expressions *)
cannam@95:     times3x3 (CE (c, d)) (CE (a, b))
cannam@95: 
cannam@95: let ctimes (CE (a, _)) (CE (c, _)) = 
cannam@95:   CE (CTimes (a, c), makeNum Number.zero)
cannam@95: 
cannam@95: let ctimesj (CE (a, _)) (CE (c, _)) = 
cannam@95:   CE (CTimesJ (a, c), makeNum Number.zero)
cannam@95:       
cannam@95: (* complex exponential (of root of unity); returns exp(2*pi*i/n * m) *)
cannam@95: let exp n i =
cannam@95:   let (c, s) = Number.cexp n i
cannam@95:   in CE (makeNum c, makeNum s)
cannam@95: 
cannam@95: (* various trig functions evaluated at (2*pi*i/n * m) *)
cannam@95: let sec n m =
cannam@95:   let (c, s) = Number.cexp n m
cannam@95:   in CE (makeNum (Number.div Number.one c), makeNum Number.zero)
cannam@95: let csc n m =
cannam@95:   let (c, s) = Number.cexp n m
cannam@95:   in CE (makeNum (Number.div Number.one s), makeNum Number.zero)
cannam@95: let tan n m =
cannam@95:   let (c, s) = Number.cexp n m
cannam@95:   in CE (makeNum (Number.div s c), makeNum Number.zero)
cannam@95: let cot n m =
cannam@95:   let (c, s) = Number.cexp n m
cannam@95:   in CE (makeNum (Number.div c s), makeNum Number.zero)
cannam@95:     
cannam@95: (* complex sum *)
cannam@95: let plus a =
cannam@95:   let rec unzip_complex = function
cannam@95:       [] -> ([], [])
cannam@95:     | ((CE (a, b)) :: s) ->
cannam@95:         let (r,i) = unzip_complex s
cannam@95: 	in
cannam@95: 	(a::r), (b::i) in
cannam@95:   let (c, d) = unzip_complex a in
cannam@95:   CE (makePlus c, makePlus d)
cannam@95: 
cannam@95: (* extract real/imaginary *)
cannam@95: let real (CE (a, b)) = CE (a, makeNum Number.zero)
cannam@95: let imag (CE (a, b)) = CE (b, makeNum Number.zero)
cannam@95: let iimag (CE (a, b)) = CE (makeNum Number.zero, b)
cannam@95: let conj (CE (a, b)) = CE (a, makeUminus b)
cannam@95: 
cannam@95:     
cannam@95: (* abstraction of sum_{i=0}^{n-1} *)
cannam@95: let sigma a b f = plus (List.map f (Util.interval a b))
cannam@95: 
cannam@95: (* store and assignment operations *)
cannam@95: let store_real v (CE (a, b)) = Expr.Store (v, a)
cannam@95: let store_imag v (CE (a, b)) = Expr.Store (v, b)
cannam@95: let store (vr, vi) x = (store_real vr x, store_imag vi x)
cannam@95: 
cannam@95: let assign_real v (CE (a, b)) = Expr.Assign (v, a)
cannam@95: let assign_imag v (CE (a, b)) = Expr.Assign (v, b)
cannam@95: let assign (vr, vi) x = (assign_real vr x, assign_imag vi x)
cannam@95: 
cannam@95: 
cannam@95: (************************
cannam@95:    shortcuts 
cannam@95:  ************************)
cannam@95: let (@*) = times
cannam@95: let (@+) a b = plus [a; b]
cannam@95: let (@-) a b = plus [a; uminus b]
cannam@95: 
cannam@95: (* type of complex signals *)
cannam@95: type signal = int -> expr
cannam@95: 
cannam@95: (* make a finite signal infinite *)
cannam@95: let infinite n signal i = if ((0 <= i) && (i < n)) then signal i else zero
cannam@95: 
cannam@95: let hermitian n a =
cannam@95:   Util.array n (fun i ->
cannam@95:     if (i = 0) then real (a 0)
cannam@95:     else if (i < n - i)  then (a i)
cannam@95:     else if (i > n - i)  then conj (a (n - i))
cannam@95:     else real (a i))
cannam@95: 
cannam@95: let antihermitian n a =
cannam@95:   Util.array n (fun i ->
cannam@95:     if (i = 0) then iimag (a 0)
cannam@95:     else if (i < n - i)  then (a i)
cannam@95:     else if (i > n - i)  then uminus (conj (a (n - i)))
cannam@95:     else iimag (a i))