diff options
| -rw-r--r-- | src/clojure/contrib/complex_numbers.clj | 48 | ||||
| -rw-r--r-- | src/clojure/contrib/test_contrib/complex_numbers.clj | 30 |
2 files changed, 66 insertions, 12 deletions
diff --git a/src/clojure/contrib/complex_numbers.clj b/src/clojure/contrib/complex_numbers.clj index 1a05934a..942051da 100644 --- a/src/clojure/contrib/complex_numbers.clj +++ b/src/clojure/contrib/complex_numbers.clj @@ -1,7 +1,7 @@ ;; Complex numbers ;; by Konrad Hinsen -;; last updated March 23, 2009 +;; last updated April 2, 2009 ;; Copyright (c) Konrad Hinsen, 2009. All rights reserved. The use ;; and distribution terms for this software are covered by the Eclipse @@ -249,10 +249,42 @@ ; ; Square root ; -(defmethod gm/sqrt ::complex - [{r :real i :imag}] - (let [abs (gm/sqrt (ga/+ (ga/* r r) (ga/* i i))) - p (gm/sqrt ((ga/qsym ga /) (ga/+ abs r) 2)) - q (ga/* (gm/sgn i) - (gm/sqrt ((ga/qsym ga /) (ga/- abs r) 2)))] - (complex p q))) +(let [one-half (/ 1 2) + one-eighth (/ 1 8)] + (defmethod gm/sqrt ::complex + [x] + (let [[r i] (vals x)] + (if (and (gc/zero? r) (gc/zero? i)) + 0 + (let [; The basic formula would say + ; abs (gm/sqrt (ga/+ (ga/* r r) (ga/* i i))) + ; p (gm/sqrt (ga/* one-half (ga/+ abs r))) + ; but the slightly more complicated one below + ; avoids overflow for large r or i. + ar (gm/abs r) + ai (gm/abs i) + r8 (ga/* one-eighth ar) + i8 (ga/* one-eighth ai) + abs (gm/sqrt (ga/+ (ga/* r8 r8) (ga/* i8 i8))) + p (ga/* 2 (gm/sqrt (ga/+ abs r8))) + q ((ga/qsym ga /) ai (ga/* 2 p)) + s (gm/sgn i)] + (if (gc/< r 0) + (complex q (ga/* s p)) + (complex p (ga/* s q)))))))) + +; +; Exponential function +; +(defmethod gm/exp ::complex + [x] + (let [[r i] (vals x) + exp-r (gm/exp r) + cos-i (gm/cos i) + sin-i (gm/sin i)] + (complex (ga/* exp-r cos-i) (ga/* exp-r sin-i)))) + +(defmethod gm/exp ::pure-imaginary + [x] + (let [i (imag x)] + (complex (gm/cos i) (gm/sin i)))) diff --git a/src/clojure/contrib/test_contrib/complex_numbers.clj b/src/clojure/contrib/test_contrib/complex_numbers.clj index b8b03f49..2ac2ff86 100644 --- a/src/clojure/contrib/test_contrib/complex_numbers.clj +++ b/src/clojure/contrib/test_contrib/complex_numbers.clj @@ -1,7 +1,7 @@ ;; Test routines for complex-numbers.clj ;; by Konrad Hinsen -;; last updated March 23, 2009 +;; last updated April 2, 2009 ;; Copyright (c) Konrad Hinsen, 2008. All rights reserved. The use ;; and distribution terms for this software are covered by the Eclipse @@ -20,7 +20,7 @@ [clojure.contrib.generic.comparison :only (= < > <= >=)] [clojure.contrib.generic.math-functions - :only (abs approx= conjugate sqr sqrt)] + :only (abs approx= conjugate exp sqr sqrt)] [clojure.contrib.complex-numbers :only (complex imaginary real imag)])) @@ -279,13 +279,35 @@ (is (= (conjugate (imaginary 5)) (imaginary -5)))) (deftest complex-abs - (doseq [c [(complex 1 2) (complex -3 -7) (imaginary -2) (imaginary 5)]] + (doseq [c [(complex 1 2) (complex -2 3) (complex 4 -2) + (complex -3 -7) (imaginary -2) (imaginary 5)]] (is (approx= (* c (conjugate c)) (sqr (abs c)) 1e-14)))) (deftest complex-sqrt - (doseq [c [(complex 1 2) (complex -3 -7) (imaginary -2) (imaginary 5)]] + (doseq [c [(complex 1 2) (complex -2 3) (complex 4 -2) + (complex -3 -7) (imaginary -2) (imaginary 5)]] (let [r (sqrt c)] (is (approx= c (sqr r) 1e-14)) (is (>= (real r) 0))))) + +(deftest complex-exp + (is (approx= (exp (complex 1 2)) + (complex -1.1312043837568135 2.4717266720048188) + 1e-14)) + (is (approx= (exp (complex 2 3)) + (complex -7.3151100949011028 1.0427436562359045) + 1e-14)) + (is (approx= (exp (complex 4 -2)) + (complex -22.720847417619233 -49.645957334580565) + 1e-14)) + (is (approx= (exp (complex 3 -7)) + (complex 15.142531566086868 -13.195928586605717) + 1e-14)) + (is (approx= (exp (imaginary -2)) + (complex -0.41614683654714241 -0.90929742682568171) + 1e-14)) + (is (approx= (exp (imaginary 5)) + (complex 0.2836621854632263 -0.95892427466313845) + 1e-14))) |