← Back to branch summary

~diresu/blender/blender-command-port

« back to all changes in this revision

Viewing changes to extern/fftw/genfft/complex.ml

  • Committer: theeth
  • Date: 2008-10-14 16:52:04 UTC
  • Revision ID: vcs-imports@canonical.com-20081014165204-r32w2gm6s0osvdhn
copy back trunk

Show diffs side-by-side

added added

removed removed

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