Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/genfft/complex.ml @ 84:08ae793730bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add null config files
author Chris Cannam
date Mon, 02 Mar 2020 14:03:47 +0000
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 (*
Chris@42 2 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
Chris@42 3 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 4 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 5 *
Chris@42 6 * This program is free software; you can redistribute it and/or modify
Chris@42 7 * it under the terms of the GNU General Public License as published by
Chris@42 8 * the Free Software Foundation; either version 2 of the License, or
Chris@42 9 * (at your option) any later version.
Chris@42 10 *
Chris@42 11 * This program is distributed in the hope that it will be useful,
Chris@42 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 14 * GNU General Public License for more details.
Chris@42 15 *
Chris@42 16 * You should have received a copy of the GNU General Public License
Chris@42 17 * along with this program; if not, write to the Free Software
Chris@42 18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 19 *
Chris@42 20 *)
Chris@42 21
Chris@42 22 (* abstraction layer for complex operations *)
Chris@42 23 open Littlesimp
Chris@42 24 open Expr
Chris@42 25
Chris@42 26 (* type of complex expressions *)
Chris@42 27 type expr = CE of Expr.expr * Expr.expr
Chris@42 28
Chris@42 29 let two = CE (makeNum Number.two, makeNum Number.zero)
Chris@42 30 let one = CE (makeNum Number.one, makeNum Number.zero)
Chris@42 31 let i = CE (makeNum Number.zero, makeNum Number.one)
Chris@42 32 let zero = CE (makeNum Number.zero, makeNum Number.zero)
Chris@42 33 let make (r, i) = CE (r, i)
Chris@42 34
Chris@42 35 let uminus (CE (a, b)) = CE (makeUminus a, makeUminus b)
Chris@42 36
Chris@42 37 let inverse_int n = CE (makeNum (Number.div Number.one (Number.of_int n)),
Chris@42 38 makeNum Number.zero)
Chris@42 39
Chris@42 40 let inverse_int_sqrt n =
Chris@42 41 CE (makeNum (Number.div Number.one (Number.sqrt (Number.of_int n))),
Chris@42 42 makeNum Number.zero)
Chris@42 43 let int_sqrt n =
Chris@42 44 CE (makeNum (Number.sqrt (Number.of_int n)),
Chris@42 45 makeNum Number.zero)
Chris@42 46
Chris@42 47 let nan x = CE (NaN x, makeNum Number.zero)
Chris@42 48
Chris@42 49 let half = inverse_int 2
Chris@42 50
Chris@42 51 let times3x3 (CE (a, b)) (CE (c, d)) =
Chris@42 52 CE (makePlus [makeTimes (c, makePlus [a; makeUminus (b)]);
Chris@42 53 makeTimes (b, makePlus [c; makeUminus (d)])],
Chris@42 54 makePlus [makeTimes (a, makePlus [c; d]);
Chris@42 55 makeUminus(makeTimes (c, makePlus [a; makeUminus (b)]))])
Chris@42 56
Chris@42 57 let times (CE (a, b)) (CE (c, d)) =
Chris@42 58 if not !Magic.threemult then
Chris@42 59 CE (makePlus [makeTimes (a, c); makeUminus (makeTimes (b, d))],
Chris@42 60 makePlus [makeTimes (a, d); makeTimes (b, c)])
Chris@42 61 else if is_constant c && is_constant d then
Chris@42 62 times3x3 (CE (a, b)) (CE (c, d))
Chris@42 63 else (* hope a and b are constant expressions *)
Chris@42 64 times3x3 (CE (c, d)) (CE (a, b))
Chris@42 65
Chris@42 66 let ctimes (CE (a, _)) (CE (c, _)) =
Chris@42 67 CE (CTimes (a, c), makeNum Number.zero)
Chris@42 68
Chris@42 69 let ctimesj (CE (a, _)) (CE (c, _)) =
Chris@42 70 CE (CTimesJ (a, c), makeNum Number.zero)
Chris@42 71
Chris@42 72 (* complex exponential (of root of unity); returns exp(2*pi*i/n * m) *)
Chris@42 73 let exp n i =
Chris@42 74 let (c, s) = Number.cexp n i
Chris@42 75 in CE (makeNum c, makeNum s)
Chris@42 76
Chris@42 77 (* various trig functions evaluated at (2*pi*i/n * m) *)
Chris@42 78 let sec n m =
Chris@42 79 let (c, s) = Number.cexp n m
Chris@42 80 in CE (makeNum (Number.div Number.one c), makeNum Number.zero)
Chris@42 81 let csc n m =
Chris@42 82 let (c, s) = Number.cexp n m
Chris@42 83 in CE (makeNum (Number.div Number.one s), makeNum Number.zero)
Chris@42 84 let tan n m =
Chris@42 85 let (c, s) = Number.cexp n m
Chris@42 86 in CE (makeNum (Number.div s c), makeNum Number.zero)
Chris@42 87 let cot n m =
Chris@42 88 let (c, s) = Number.cexp n m
Chris@42 89 in CE (makeNum (Number.div c s), makeNum Number.zero)
Chris@42 90
Chris@42 91 (* complex sum *)
Chris@42 92 let plus a =
Chris@42 93 let rec unzip_complex = function
Chris@42 94 [] -> ([], [])
Chris@42 95 | ((CE (a, b)) :: s) ->
Chris@42 96 let (r,i) = unzip_complex s
Chris@42 97 in
Chris@42 98 (a::r), (b::i) in
Chris@42 99 let (c, d) = unzip_complex a in
Chris@42 100 CE (makePlus c, makePlus d)
Chris@42 101
Chris@42 102 (* extract real/imaginary *)
Chris@42 103 let real (CE (a, b)) = CE (a, makeNum Number.zero)
Chris@42 104 let imag (CE (a, b)) = CE (b, makeNum Number.zero)
Chris@42 105 let iimag (CE (a, b)) = CE (makeNum Number.zero, b)
Chris@42 106 let conj (CE (a, b)) = CE (a, makeUminus b)
Chris@42 107
Chris@42 108
Chris@42 109 (* abstraction of sum_{i=0}^{n-1} *)
Chris@42 110 let sigma a b f = plus (List.map f (Util.interval a b))
Chris@42 111
Chris@42 112 (* store and assignment operations *)
Chris@42 113 let store_real v (CE (a, b)) = Expr.Store (v, a)
Chris@42 114 let store_imag v (CE (a, b)) = Expr.Store (v, b)
Chris@42 115 let store (vr, vi) x = (store_real vr x, store_imag vi x)
Chris@42 116
Chris@42 117 let assign_real v (CE (a, b)) = Expr.Assign (v, a)
Chris@42 118 let assign_imag v (CE (a, b)) = Expr.Assign (v, b)
Chris@42 119 let assign (vr, vi) x = (assign_real vr x, assign_imag vi x)
Chris@42 120
Chris@42 121
Chris@42 122 (************************
Chris@42 123 shortcuts
Chris@42 124 ************************)
Chris@42 125 let (@*) = times
Chris@42 126 let (@+) a b = plus [a; b]
Chris@42 127 let (@-) a b = plus [a; uminus b]
Chris@42 128
Chris@42 129 (* type of complex signals *)
Chris@42 130 type signal = int -> expr
Chris@42 131
Chris@42 132 (* make a finite signal infinite *)
Chris@42 133 let infinite n signal i = if ((0 <= i) && (i < n)) then signal i else zero
Chris@42 134
Chris@42 135 let hermitian n a =
Chris@42 136 Util.array n (fun i ->
Chris@42 137 if (i = 0) then real (a 0)
Chris@42 138 else if (i < n - i) then (a i)
Chris@42 139 else if (i > n - i) then conj (a (n - i))
Chris@42 140 else real (a i))
Chris@42 141
Chris@42 142 let antihermitian n a =
Chris@42 143 Util.array n (fun i ->
Chris@42 144 if (i = 0) then iimag (a 0)
Chris@42 145 else if (i < n - i) then (a i)
Chris@42 146 else if (i > n - i) then uminus (conj (a (n - i)))
Chris@42 147 else iimag (a i))