Mercurial > lasercutter

comparison src/clojure/contrib/complex_numbers.clj @ 10:ef7dbbd6452c

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