|
rlm@10
|
1 ;; Complex numbers
|
|
rlm@10
|
2
|
|
rlm@10
|
3 ;; by Konrad Hinsen
|
|
rlm@10
|
4 ;; last updated May 4, 2009
|
|
rlm@10
|
5
|
|
rlm@10
|
6 ;; Copyright (c) Konrad Hinsen, 2009. All rights reserved. The use
|
|
rlm@10
|
7 ;; and distribution terms for this software are covered by the Eclipse
|
|
rlm@10
|
8 ;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
|
|
rlm@10
|
9 ;; which can be found in the file epl-v10.html at the root of this
|
|
rlm@10
|
10 ;; distribution. By using this software in any fashion, you are
|
|
rlm@10
|
11 ;; agreeing to be bound by the terms of this license. You must not
|
|
rlm@10
|
12 ;; remove this notice, or any other, from this software.
|
|
rlm@10
|
13
|
|
rlm@10
|
14 (ns
|
|
rlm@10
|
15 ^{:author "Konrad Hinsen"
|
|
rlm@10
|
16 :doc "Complex numbers
|
|
rlm@10
|
17 NOTE: This library is in evolution. Most math functions are
|
|
rlm@10
|
18 not implemented yet."}
|
|
rlm@10
|
19 clojure.contrib.complex-numbers
|
|
rlm@10
|
20 (:refer-clojure :exclude (deftype))
|
|
rlm@10
|
21 (:use [clojure.contrib.types :only (deftype)]
|
|
rlm@10
|
22 [clojure.contrib.generic :only (root-type)])
|
|
rlm@10
|
23 (:require [clojure.contrib.generic.arithmetic :as ga]
|
|
rlm@10
|
24 [clojure.contrib.generic.comparison :as gc]
|
|
rlm@10
|
25 [clojure.contrib.generic.math-functions :as gm]))
|
|
rlm@10
|
26
|
|
rlm@10
|
27 ;
|
|
rlm@10
|
28 ; Complex numbers are represented as struct maps. The real and imaginary
|
|
rlm@10
|
29 ; parts can be of any type for which arithmetic and maths functions
|
|
rlm@10
|
30 ; are defined.
|
|
rlm@10
|
31 ;
|
|
rlm@10
|
32 (defstruct complex-struct :real :imag)
|
|
rlm@10
|
33
|
|
rlm@10
|
34 ;
|
|
rlm@10
|
35 ; The general complex number type
|
|
rlm@10
|
36 ;
|
|
rlm@10
|
37 (deftype ::complex complex
|
|
rlm@10
|
38 (fn [real imag] (struct complex-struct real imag))
|
|
rlm@10
|
39 (fn [c] (vals c)))
|
|
rlm@10
|
40
|
|
rlm@10
|
41 (derive ::complex root-type)
|
|
rlm@10
|
42
|
|
rlm@10
|
43 ;
|
|
rlm@10
|
44 ; A specialized subtype for pure imaginary numbers. Introducing this type
|
|
rlm@10
|
45 ; reduces the number of operations by eliminating additions with and
|
|
rlm@10
|
46 ; multiplications by zero.
|
|
rlm@10
|
47 ;
|
|
rlm@10
|
48 (deftype ::pure-imaginary imaginary
|
|
rlm@10
|
49 (fn [imag] (struct complex-struct 0 imag))
|
|
rlm@10
|
50 (fn [c] (list (:imag c))))
|
|
rlm@10
|
51
|
|
rlm@10
|
52 (derive ::pure-imaginary ::complex)
|
|
rlm@10
|
53
|
|
rlm@10
|
54 ;
|
|
rlm@10
|
55 ; Extraction of real and imaginary parts
|
|
rlm@10
|
56 ;
|
|
rlm@10
|
57 (def real (accessor complex-struct :real))
|
|
rlm@10
|
58 (def imag (accessor complex-struct :imag))
|
|
rlm@10
|
59
|
|
rlm@10
|
60 ;
|
|
rlm@10
|
61 ; Equality and zero test
|
|
rlm@10
|
62 ;
|
|
rlm@10
|
63 (defmethod gc/zero? ::complex
|
|
rlm@10
|
64 [x]
|
|
rlm@10
|
65 (let [[rx ix] (vals x)]
|
|
rlm@10
|
66 (and (zero? rx) (zero? ix))))
|
|
rlm@10
|
67
|
|
rlm@10
|
68 (defmethod gc/= [::complex ::complex]
|
|
rlm@10
|
69 [x y]
|
|
rlm@10
|
70 (let [[rx ix] (vals x)
|
|
rlm@10
|
71 [ry iy] (vals y)]
|
|
rlm@10
|
72 (and (gc/= rx ry) (gc/= ix iy))))
|
|
rlm@10
|
73
|
|
rlm@10
|
74 (defmethod gc/= [::pure-imaginary ::pure-imaginary]
|
|
rlm@10
|
75 [x y]
|
|
rlm@10
|
76 (gc/= (imag x) (imag y)))
|
|
rlm@10
|
77
|
|
rlm@10
|
78 (defmethod gc/= [::complex ::pure-imaginary]
|
|
rlm@10
|
79 [x y]
|
|
rlm@10
|
80 (let [[rx ix] (vals x)]
|
|
rlm@10
|
81 (and (gc/zero? rx) (gc/= ix (imag y)))))
|
|
rlm@10
|
82
|
|
rlm@10
|
83 (defmethod gc/= [::pure-imaginary ::complex]
|
|
rlm@10
|
84 [x y]
|
|
rlm@10
|
85 (let [[ry iy] (vals y)]
|
|
rlm@10
|
86 (and (gc/zero? ry) (gc/= (imag x) iy))))
|
|
rlm@10
|
87
|
|
rlm@10
|
88 (defmethod gc/= [::complex root-type]
|
|
rlm@10
|
89 [x y]
|
|
rlm@10
|
90 (let [[rx ix] (vals x)]
|
|
rlm@10
|
91 (and (gc/zero? ix) (gc/= rx y))))
|
|
rlm@10
|
92
|
|
rlm@10
|
93 (defmethod gc/= [root-type ::complex]
|
|
rlm@10
|
94 [x y]
|
|
rlm@10
|
95 (let [[ry iy] (vals y)]
|
|
rlm@10
|
96 (and (gc/zero? iy) (gc/= x ry))))
|
|
rlm@10
|
97
|
|
rlm@10
|
98 (defmethod gc/= [::pure-imaginary root-type]
|
|
rlm@10
|
99 [x y]
|
|
rlm@10
|
100 (and (gc/zero? (imag x)) (gc/zero? y)))
|
|
rlm@10
|
101
|
|
rlm@10
|
102 (defmethod gc/= [root-type ::pure-imaginary]
|
|
rlm@10
|
103 [x y]
|
|
rlm@10
|
104 (and (gc/zero? x) (gc/zero? (imag y))))
|
|
rlm@10
|
105
|
|
rlm@10
|
106 ;
|
|
rlm@10
|
107 ; Addition
|
|
rlm@10
|
108 ;
|
|
rlm@10
|
109 (defmethod ga/+ [::complex ::complex]
|
|
rlm@10
|
110 [x y]
|
|
rlm@10
|
111 (let [[rx ix] (vals x)
|
|
rlm@10
|
112 [ry iy] (vals y)]
|
|
rlm@10
|
113 (complex (ga/+ rx ry) (ga/+ ix iy))))
|
|
rlm@10
|
114
|
|
rlm@10
|
115 (defmethod ga/+ [::pure-imaginary ::pure-imaginary]
|
|
rlm@10
|
116 [x y]
|
|
rlm@10
|
117 (imaginary (ga/+ (imag x) (imag y))))
|
|
rlm@10
|
118
|
|
rlm@10
|
119 (defmethod ga/+ [::complex ::pure-imaginary]
|
|
rlm@10
|
120 [x y]
|
|
rlm@10
|
121 (let [[rx ix] (vals x)]
|
|
rlm@10
|
122 (complex rx (ga/+ ix (imag y)))))
|
|
rlm@10
|
123
|
|
rlm@10
|
124 (defmethod ga/+ [::pure-imaginary ::complex]
|
|
rlm@10
|
125 [x y]
|
|
rlm@10
|
126 (let [[ry iy] (vals y)]
|
|
rlm@10
|
127 (complex ry (ga/+ (imag x) iy))))
|
|
rlm@10
|
128
|
|
rlm@10
|
129 (defmethod ga/+ [::complex root-type]
|
|
rlm@10
|
130 [x y]
|
|
rlm@10
|
131 (let [[rx ix] (vals x)]
|
|
rlm@10
|
132 (complex (ga/+ rx y) ix)))
|
|
rlm@10
|
133
|
|
rlm@10
|
134 (defmethod ga/+ [root-type ::complex]
|
|
rlm@10
|
135 [x y]
|
|
rlm@10
|
136 (let [[ry iy] (vals y)]
|
|
rlm@10
|
137 (complex (ga/+ x ry) iy)))
|
|
rlm@10
|
138
|
|
rlm@10
|
139 (defmethod ga/+ [::pure-imaginary root-type]
|
|
rlm@10
|
140 [x y]
|
|
rlm@10
|
141 (complex y (imag x)))
|
|
rlm@10
|
142
|
|
rlm@10
|
143 (defmethod ga/+ [root-type ::pure-imaginary]
|
|
rlm@10
|
144 [x y]
|
|
rlm@10
|
145 (complex x (imag y)))
|
|
rlm@10
|
146
|
|
rlm@10
|
147 ;
|
|
rlm@10
|
148 ; Negation
|
|
rlm@10
|
149 ;
|
|
rlm@10
|
150 (defmethod ga/- ::complex
|
|
rlm@10
|
151 [x]
|
|
rlm@10
|
152 (let [[rx ix] (vals x)]
|
|
rlm@10
|
153 (complex (ga/- rx) (ga/- ix))))
|
|
rlm@10
|
154
|
|
rlm@10
|
155 (defmethod ga/- ::pure-imaginary
|
|
rlm@10
|
156 [x]
|
|
rlm@10
|
157 (imaginary (ga/- (imag x))))
|
|
rlm@10
|
158
|
|
rlm@10
|
159 ;
|
|
rlm@10
|
160 ; Subtraction is automatically supplied by ga/-, optimized implementations
|
|
rlm@10
|
161 ; can be added later...
|
|
rlm@10
|
162 ;
|
|
rlm@10
|
163
|
|
rlm@10
|
164 ;
|
|
rlm@10
|
165 ; Multiplication
|
|
rlm@10
|
166 ;
|
|
rlm@10
|
167 (defmethod ga/* [::complex ::complex]
|
|
rlm@10
|
168 [x y]
|
|
rlm@10
|
169 (let [[rx ix] (vals x)
|
|
rlm@10
|
170 [ry iy] (vals y)]
|
|
rlm@10
|
171 (complex (ga/- (ga/* rx ry) (ga/* ix iy))
|
|
rlm@10
|
172 (ga/+ (ga/* rx iy) (ga/* ix ry)))))
|
|
rlm@10
|
173
|
|
rlm@10
|
174 (defmethod ga/* [::pure-imaginary ::pure-imaginary]
|
|
rlm@10
|
175 [x y]
|
|
rlm@10
|
176 (ga/- (ga/* (imag x) (imag y))))
|
|
rlm@10
|
177
|
|
rlm@10
|
178 (defmethod ga/* [::complex ::pure-imaginary]
|
|
rlm@10
|
179 [x y]
|
|
rlm@10
|
180 (let [[rx ix] (vals x)
|
|
rlm@10
|
181 iy (imag y)]
|
|
rlm@10
|
182 (complex (ga/- (ga/* ix iy))
|
|
rlm@10
|
183 (ga/* rx iy))))
|
|
rlm@10
|
184
|
|
rlm@10
|
185 (defmethod ga/* [::pure-imaginary ::complex]
|
|
rlm@10
|
186 [x y]
|
|
rlm@10
|
187 (let [ix (imag x)
|
|
rlm@10
|
188 [ry iy] (vals y)]
|
|
rlm@10
|
189 (complex (ga/- (ga/* ix iy))
|
|
rlm@10
|
190 (ga/* ix ry))))
|
|
rlm@10
|
191
|
|
rlm@10
|
192 (defmethod ga/* [::complex root-type]
|
|
rlm@10
|
193 [x y]
|
|
rlm@10
|
194 (let [[rx ix] (vals x)]
|
|
rlm@10
|
195 (complex (ga/* rx y) (ga/* ix y))))
|
|
rlm@10
|
196
|
|
rlm@10
|
197 (defmethod ga/* [root-type ::complex]
|
|
rlm@10
|
198 [x y]
|
|
rlm@10
|
199 (let [[ry iy] (vals y)]
|
|
rlm@10
|
200 (complex (ga/* x ry) (ga/* x iy))))
|
|
rlm@10
|
201
|
|
rlm@10
|
202 (defmethod ga/* [::pure-imaginary root-type]
|
|
rlm@10
|
203 [x y]
|
|
rlm@10
|
204 (imaginary (ga/* (imag x) y)))
|
|
rlm@10
|
205
|
|
rlm@10
|
206 (defmethod ga/* [root-type ::pure-imaginary]
|
|
rlm@10
|
207 [x y]
|
|
rlm@10
|
208 (imaginary (ga/* x (imag y))))
|
|
rlm@10
|
209
|
|
rlm@10
|
210 ;
|
|
rlm@10
|
211 ; Inversion
|
|
rlm@10
|
212 ;
|
|
rlm@10
|
213 (ga/defmethod* ga / ::complex
|
|
rlm@10
|
214 [x]
|
|
rlm@10
|
215 (let [[rx ix] (vals x)
|
|
rlm@10
|
216 den ((ga/qsym ga /) (ga/+ (ga/* rx rx) (ga/* ix ix)))]
|
|
rlm@10
|
217 (complex (ga/* rx den) (ga/- (ga/* ix den)))))
|
|
rlm@10
|
218
|
|
rlm@10
|
219 (ga/defmethod* ga / ::pure-imaginary
|
|
rlm@10
|
220 [x]
|
|
rlm@10
|
221 (imaginary (ga/- ((ga/qsym ga /) (imag x)))))
|
|
rlm@10
|
222
|
|
rlm@10
|
223 ;
|
|
rlm@10
|
224 ; Division is automatically supplied by ga//, optimized implementations
|
|
rlm@10
|
225 ; can be added later...
|
|
rlm@10
|
226 ;
|
|
rlm@10
|
227
|
|
rlm@10
|
228 ;
|
|
rlm@10
|
229 ; Conjugation
|
|
rlm@10
|
230 ;
|
|
rlm@10
|
231 (defmethod gm/conjugate ::complex
|
|
rlm@10
|
232 [x]
|
|
rlm@10
|
233 (let [[r i] (vals x)]
|
|
rlm@10
|
234 (complex r (ga/- i))))
|
|
rlm@10
|
235
|
|
rlm@10
|
236 (defmethod gm/conjugate ::pure-imaginary
|
|
rlm@10
|
237 [x]
|
|
rlm@10
|
238 (imaginary (ga/- (imag x))))
|
|
rlm@10
|
239
|
|
rlm@10
|
240 ;
|
|
rlm@10
|
241 ; Absolute value
|
|
rlm@10
|
242 ;
|
|
rlm@10
|
243 (defmethod gm/abs ::complex
|
|
rlm@10
|
244 [x]
|
|
rlm@10
|
245 (let [[r i] (vals x)]
|
|
rlm@10
|
246 (gm/sqrt (ga/+ (ga/* r r) (ga/* i i)))))
|
|
rlm@10
|
247
|
|
rlm@10
|
248 (defmethod gm/abs ::pure-imaginary
|
|
rlm@10
|
249 [x]
|
|
rlm@10
|
250 (gm/abs (imag x)))
|
|
rlm@10
|
251
|
|
rlm@10
|
252 ;
|
|
rlm@10
|
253 ; Square root
|
|
rlm@10
|
254 ;
|
|
rlm@10
|
255 (let [one-half (/ 1 2)
|
|
rlm@10
|
256 one-eighth (/ 1 8)]
|
|
rlm@10
|
257 (defmethod gm/sqrt ::complex
|
|
rlm@10
|
258 [x]
|
|
rlm@10
|
259 (let [[r i] (vals x)]
|
|
rlm@10
|
260 (if (and (gc/zero? r) (gc/zero? i))
|
|
rlm@10
|
261 0
|
|
rlm@10
|
262 (let [; The basic formula would say
|
|
rlm@10
|
263 ; abs (gm/sqrt (ga/+ (ga/* r r) (ga/* i i)))
|
|
rlm@10
|
264 ; p (gm/sqrt (ga/* one-half (ga/+ abs r)))
|
|
rlm@10
|
265 ; but the slightly more complicated one below
|
|
rlm@10
|
266 ; avoids overflow for large r or i.
|
|
rlm@10
|
267 ar (gm/abs r)
|
|
rlm@10
|
268 ai (gm/abs i)
|
|
rlm@10
|
269 r8 (ga/* one-eighth ar)
|
|
rlm@10
|
270 i8 (ga/* one-eighth ai)
|
|
rlm@10
|
271 abs (gm/sqrt (ga/+ (ga/* r8 r8) (ga/* i8 i8)))
|
|
rlm@10
|
272 p (ga/* 2 (gm/sqrt (ga/+ abs r8)))
|
|
rlm@10
|
273 q ((ga/qsym ga /) ai (ga/* 2 p))
|
|
rlm@10
|
274 s (gm/sgn i)]
|
|
rlm@10
|
275 (if (gc/< r 0)
|
|
rlm@10
|
276 (complex q (ga/* s p))
|
|
rlm@10
|
277 (complex p (ga/* s q))))))))
|
|
rlm@10
|
278
|
|
rlm@10
|
279 ;
|
|
rlm@10
|
280 ; Exponential function
|
|
rlm@10
|
281 ;
|
|
rlm@10
|
282 (defmethod gm/exp ::complex
|
|
rlm@10
|
283 [x]
|
|
rlm@10
|
284 (let [[r i] (vals x)
|
|
rlm@10
|
285 exp-r (gm/exp r)
|
|
rlm@10
|
286 cos-i (gm/cos i)
|
|
rlm@10
|
287 sin-i (gm/sin i)]
|
|
rlm@10
|
288 (complex (ga/* exp-r cos-i) (ga/* exp-r sin-i))))
|
|
rlm@10
|
289
|
|
rlm@10
|
290 (defmethod gm/exp ::pure-imaginary
|
|
rlm@10
|
291 [x]
|
|
rlm@10
|
292 (let [i (imag x)]
|
|
rlm@10
|
293 (complex (gm/cos i) (gm/sin i))))
|
