Mercurial > lasercutter

     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))))
author	Robert McIntyre <rlm@mit.edu>
date	Sat, 21 Aug 2010 06:25:44 -0400
parents
children