1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/clojure/contrib/complex_numbers.clj	Sat Aug 21 06:25:44 2010 -0400
     1.3 @@ -0,0 +1,293 @@
     1.4 +;; Complex numbers
     1.5 +
     1.6 +;; by Konrad Hinsen
     1.7 +;; last updated May 4, 2009
     1.8 +
     1.9 +;; Copyright (c) Konrad Hinsen, 2009. All rights reserved.  The use
    1.10 +;; and distribution terms for this software are covered by the Eclipse
    1.11 +;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
    1.12 +;; which can be found in the file epl-v10.html at the root of this
    1.13 +;; distribution.  By using this software in any fashion, you are
    1.14 +;; agreeing to be bound by the terms of this license.  You must not
    1.15 +;; remove this notice, or any other, from this software.
    1.16 +
    1.17 +(ns
    1.18 +  ^{:author "Konrad Hinsen"
    1.19 +     :doc "Complex numbers
    1.20 +           NOTE: This library is in evolution. Most math functions are
    1.21 +                 not implemented yet."}
    1.22 +  clojure.contrib.complex-numbers
    1.23 +  (:refer-clojure :exclude (deftype))
    1.24 +  (:use [clojure.contrib.types :only (deftype)]
    1.25 +	[clojure.contrib.generic :only (root-type)])
    1.26 +  (:require [clojure.contrib.generic.arithmetic :as ga]
    1.27 +	    [clojure.contrib.generic.comparison :as gc]
    1.28 +	    [clojure.contrib.generic.math-functions :as gm]))
    1.29 +
    1.30 +;
    1.31 +; Complex numbers are represented as struct maps. The real and imaginary
    1.32 +; parts can be of any type for which arithmetic and maths functions
    1.33 +; are defined.
    1.34 +;
    1.35 +(defstruct complex-struct :real :imag)
    1.36 +
    1.37 +;
    1.38 +; The general complex number type
    1.39 +;
    1.40 +(deftype ::complex complex
    1.41 +  (fn [real imag] (struct complex-struct real imag))
    1.42 +  (fn [c] (vals c)))
    1.43 +
    1.44 +(derive ::complex root-type)
    1.45 +
    1.46 +;
    1.47 +; A specialized subtype for pure imaginary numbers. Introducing this type
    1.48 +; reduces the number of operations by eliminating additions with and
    1.49 +; multiplications by zero.
    1.50 +;
    1.51 +(deftype ::pure-imaginary imaginary
    1.52 +  (fn [imag] (struct complex-struct 0 imag))
    1.53 +  (fn [c] (list (:imag c))))
    1.54 +
    1.55 +(derive ::pure-imaginary ::complex)
    1.56 +
    1.57 +;
    1.58 +; Extraction of real and imaginary parts
    1.59 +;
    1.60 +(def real (accessor complex-struct :real))
    1.61 +(def imag (accessor complex-struct :imag))
    1.62 +
    1.63 +;
    1.64 +; Equality and zero test
    1.65 +;
    1.66 +(defmethod gc/zero? ::complex
    1.67 +  [x]
    1.68 +  (let [[rx ix] (vals x)]
    1.69 +    (and (zero? rx) (zero? ix))))
    1.70 +
    1.71 +(defmethod gc/= [::complex ::complex]
    1.72 +  [x y]
    1.73 +  (let [[rx ix] (vals x)
    1.74 +	[ry iy] (vals y)]
    1.75 +    (and (gc/= rx ry) (gc/= ix iy))))
    1.76 +
    1.77 +(defmethod gc/= [::pure-imaginary ::pure-imaginary]
    1.78 +  [x y]
    1.79 +  (gc/= (imag x) (imag y)))
    1.80 +
    1.81 +(defmethod gc/= [::complex ::pure-imaginary]
    1.82 +  [x y]
    1.83 +  (let [[rx ix] (vals x)]
    1.84 +    (and (gc/zero? rx) (gc/= ix (imag y)))))
    1.85 +
    1.86 +(defmethod gc/= [::pure-imaginary ::complex]
    1.87 +  [x y]
    1.88 +  (let [[ry iy] (vals y)]
    1.89 +    (and (gc/zero? ry) (gc/= (imag x) iy))))
    1.90 +
    1.91 +(defmethod gc/= [::complex root-type]
    1.92 +  [x y]
    1.93 +  (let [[rx ix] (vals x)]
    1.94 +    (and (gc/zero? ix) (gc/= rx y))))
    1.95 +
    1.96 +(defmethod gc/= [root-type ::complex]
    1.97 +  [x y]
    1.98 +  (let [[ry iy] (vals y)]
    1.99 +    (and (gc/zero? iy) (gc/= x ry))))
   1.100 +
   1.101 +(defmethod gc/= [::pure-imaginary root-type]
   1.102 +  [x y]
   1.103 +  (and (gc/zero? (imag x)) (gc/zero? y)))
   1.104 +
   1.105 +(defmethod gc/= [root-type ::pure-imaginary]
   1.106 +  [x y]
   1.107 +  (and (gc/zero? x) (gc/zero? (imag y))))
   1.108 +
   1.109 +;
   1.110 +; Addition
   1.111 +;
   1.112 +(defmethod ga/+ [::complex ::complex]
   1.113 +  [x y]
   1.114 +  (let [[rx ix] (vals x)
   1.115 +	[ry iy] (vals y)]
   1.116 +    (complex (ga/+ rx ry) (ga/+ ix iy))))
   1.117 +
   1.118 +(defmethod ga/+ [::pure-imaginary ::pure-imaginary]
   1.119 +  [x y]
   1.120 +  (imaginary (ga/+ (imag x) (imag y))))
   1.121 +
   1.122 +(defmethod ga/+ [::complex ::pure-imaginary]
   1.123 +  [x y]
   1.124 +  (let [[rx ix] (vals x)]
   1.125 +    (complex rx (ga/+ ix (imag y)))))
   1.126 +
   1.127 +(defmethod ga/+ [::pure-imaginary ::complex]
   1.128 +  [x y]
   1.129 +  (let [[ry iy] (vals y)]
   1.130 +    (complex ry (ga/+ (imag x) iy))))
   1.131 +
   1.132 +(defmethod ga/+ [::complex root-type]
   1.133 +  [x y]
   1.134 +  (let [[rx ix] (vals x)]
   1.135 +    (complex (ga/+ rx y) ix)))
   1.136 +
   1.137 +(defmethod ga/+ [root-type ::complex]
   1.138 +  [x y]
   1.139 +  (let [[ry iy] (vals y)]
   1.140 +    (complex (ga/+ x ry) iy)))
   1.141 +
   1.142 +(defmethod ga/+ [::pure-imaginary root-type]
   1.143 +  [x y]
   1.144 +  (complex y (imag x)))
   1.145 +
   1.146 +(defmethod ga/+ [root-type ::pure-imaginary]
   1.147 +  [x y]
   1.148 +  (complex x (imag y)))
   1.149 +
   1.150 +;
   1.151 +; Negation
   1.152 +;
   1.153 +(defmethod ga/- ::complex
   1.154 +  [x]
   1.155 +  (let [[rx ix] (vals x)]
   1.156 +    (complex (ga/- rx) (ga/- ix))))
   1.157 +
   1.158 +(defmethod ga/- ::pure-imaginary
   1.159 +  [x]
   1.160 +  (imaginary (ga/- (imag x))))
   1.161 +
   1.162 +;
   1.163 +; Subtraction is automatically supplied by ga/-, optimized implementations
   1.164 +; can be added later...
   1.165 +;
   1.166 +
   1.167 +;
   1.168 +; Multiplication
   1.169 +;
   1.170 +(defmethod ga/* [::complex ::complex]
   1.171 +  [x y]
   1.172 +  (let [[rx ix] (vals x)
   1.173 +	[ry iy] (vals y)]
   1.174 +    (complex (ga/- (ga/* rx ry) (ga/* ix iy))
   1.175 +	     (ga/+ (ga/* rx iy) (ga/* ix ry)))))
   1.176 +
   1.177 +(defmethod ga/* [::pure-imaginary ::pure-imaginary]
   1.178 +  [x y]
   1.179 +  (ga/- (ga/* (imag x) (imag y))))
   1.180 +
   1.181 +(defmethod ga/* [::complex ::pure-imaginary]
   1.182 +  [x y]
   1.183 +  (let [[rx ix] (vals x)
   1.184 +	iy (imag y)]
   1.185 +    (complex (ga/- (ga/* ix iy))
   1.186 +	     (ga/* rx iy))))
   1.187 +
   1.188 +(defmethod ga/* [::pure-imaginary ::complex]
   1.189 +  [x y]
   1.190 +  (let [ix (imag x)
   1.191 +	[ry iy] (vals y)]
   1.192 +    (complex (ga/- (ga/* ix iy))
   1.193 +	     (ga/* ix ry))))
   1.194 +
   1.195 +(defmethod ga/* [::complex root-type]
   1.196 +  [x y]
   1.197 +  (let [[rx ix] (vals x)]
   1.198 +    (complex (ga/* rx y) (ga/* ix y))))
   1.199 +
   1.200 +(defmethod ga/* [root-type ::complex]
   1.201 +  [x y]
   1.202 +  (let [[ry iy] (vals y)]
   1.203 +    (complex (ga/* x ry) (ga/* x iy))))
   1.204 +
   1.205 +(defmethod ga/* [::pure-imaginary root-type]
   1.206 +  [x y]
   1.207 +  (imaginary (ga/* (imag x) y)))
   1.208 +
   1.209 +(defmethod ga/* [root-type ::pure-imaginary]
   1.210 +  [x y]
   1.211 +  (imaginary (ga/* x (imag y))))
   1.212 +
   1.213 +;
   1.214 +; Inversion
   1.215 +;
   1.216 +(ga/defmethod* ga / ::complex
   1.217 +  [x]
   1.218 +  (let [[rx ix] (vals x)
   1.219 +	den ((ga/qsym ga /) (ga/+ (ga/* rx rx) (ga/* ix ix)))]
   1.220 +    (complex (ga/* rx den) (ga/- (ga/* ix den)))))
   1.221 +
   1.222 +(ga/defmethod* ga / ::pure-imaginary
   1.223 +  [x]
   1.224 +  (imaginary (ga/- ((ga/qsym ga /) (imag x)))))
   1.225 +
   1.226 +;
   1.227 +; Division is automatically supplied by ga//, optimized implementations
   1.228 +; can be added later...
   1.229 +;
   1.230 +
   1.231 +;
   1.232 +; Conjugation
   1.233 +;
   1.234 +(defmethod gm/conjugate ::complex
   1.235 +  [x]
   1.236 +  (let [[r i] (vals x)]
   1.237 +    (complex r (ga/- i))))
   1.238 +
   1.239 +(defmethod gm/conjugate ::pure-imaginary
   1.240 +  [x]
   1.241 +  (imaginary (ga/- (imag x))))
   1.242 +
   1.243 +;
   1.244 +; Absolute value
   1.245 +;
   1.246 +(defmethod gm/abs ::complex
   1.247 +  [x]
   1.248 +  (let [[r i] (vals x)]
   1.249 +    (gm/sqrt (ga/+ (ga/* r r) (ga/* i i)))))
   1.250 +
   1.251 +(defmethod gm/abs ::pure-imaginary
   1.252 +  [x]
   1.253 +  (gm/abs (imag x)))
   1.254 +
   1.255 +;
   1.256 +; Square root
   1.257 +;
   1.258 +(let [one-half   (/ 1 2)
   1.259 +      one-eighth (/ 1 8)]
   1.260 +  (defmethod gm/sqrt ::complex
   1.261 +    [x]
   1.262 +    (let [[r i] (vals x)]
   1.263 +      (if (and (gc/zero? r) (gc/zero? i))
   1.264 +        0
   1.265 +        (let [; The basic formula would say
   1.266 +              ;    abs (gm/sqrt (ga/+ (ga/* r r) (ga/* i i)))
   1.267 +	      ;    p   (gm/sqrt (ga/* one-half (ga/+ abs r)))
   1.268 +	      ; but the slightly more complicated one below
   1.269 +	      ; avoids overflow for large r or i.
   1.270 +	      ar  (gm/abs r)
   1.271 +	      ai  (gm/abs i)
   1.272 +	      r8  (ga/* one-eighth ar)
   1.273 +	      i8  (ga/* one-eighth ai)
   1.274 +	      abs (gm/sqrt (ga/+ (ga/* r8 r8) (ga/* i8 i8)))
   1.275 +	      p   (ga/* 2 (gm/sqrt (ga/+ abs r8)))
   1.276 +	      q   ((ga/qsym ga /) ai (ga/* 2 p))
   1.277 +	      s   (gm/sgn i)]
   1.278 +	  (if (gc/< r 0)
   1.279 +	    (complex q (ga/* s p))
   1.280 +	    (complex p (ga/* s q))))))))
   1.281 +
   1.282 +;
   1.283 +; Exponential function
   1.284 +;
   1.285 +(defmethod gm/exp ::complex
   1.286 +  [x]
   1.287 +  (let [[r i] (vals x)
   1.288 +	exp-r (gm/exp r)
   1.289 +	cos-i (gm/cos i)
   1.290 +	sin-i (gm/sin i)]
   1.291 +    (complex (ga/* exp-r cos-i) (ga/* exp-r sin-i))))
   1.292 +
   1.293 +(defmethod gm/exp ::pure-imaginary
   1.294 +  [x]
   1.295 +  (let [i (imag x)]
   1.296 +    (complex (gm/cos i) (gm/sin i))))