hypercomplex-complex division

Author: wedesoft

TODO.md

@@ -1,6 +1,6 @@
 ## Ready
 
-* hypercomplex division
+* real-hypercomplex division
 * complex, hypercomplex exp
 * operations on rgb values
 * operations on hypercomplex values

aiscm/hypercomplex.scm

@@ -25,7 +25,7 @@
             <hypercomplex<>>
             <hypercomplex<float>>  <meta<hypercomplex<float>>>  <hypercomplex<float<single>>> <meta<hypercomplex<float<single>>>>
             <hypercomplex<double>> <meta<hypercomplex<double>>> <hypercomplex<float<double>>> <meta<hypercomplex<float<double>>>>)
-  #:re-export (real-part imag-part equal? + - abs conj))
+  #:re-export (real-part imag-part equal? + - abs conj * /))
 
 
 (define-class <hypercomplex> ()
@@ -132,3 +132,33 @@
   (hypercomplex (* (real-part a) b) (* (imag-part a) b) (* (jmag-part a) b) (* (kmag-part a) b)))
 (define-method (* (a <scalar>) (b <hypercomplex<>>))
   (hypercomplex (* a (real-part b)) (* a (imag-part b)) (* a (jmag-part b)) (* a (kmag-part b))))
+
+(define (hypercomplex-inverse a)
+  (jit-let [(d1 (- (real-part a) (kmag-part a)))
+            (d2 (+ (imag-part a) (jmag-part a)))
+            (d3 (+ (real-part a) (kmag-part a)))
+            (d4 (- (imag-part a) (jmag-part a)))
+            (denom (* (+ (* d1 d1) (* d2 d2)) (+ (* d3 d3) (* d4 d4))))
+            (squares (+ (* (real-part a) (real-part a))
+                        (* (imag-part a) (imag-part a))
+                        (* (jmag-part a) (jmag-part a))
+                        (* (kmag-part a) (kmag-part a))))
+            (cross (- (* (real-part a) (kmag-part a)) (* (imag-part a) (jmag-part a))))
+            (c1 (* (kmag-part a) cross))
+            (c2 (* (jmag-part a) cross))
+            (c3 (* (imag-part a) cross))
+            (c4 (* (real-part a) cross))]
+    (hypercomplex (/ (- (* (real-part a) squares) (+ c1 c1)) denom)
+                  (/ (- (- (* (imag-part a) squares)) (+ c2 c2)) denom)
+                  (/ (- (- (* (jmag-part a) squares)) (+ c3 c3)) denom)
+                  (/ (- (* (kmag-part a) squares) (+ c4 c4)) denom))))
+(define (complex-inverse a)
+  (jit-let [(denom (+ (* (real-part a) (real-part a)) (* (imag-part a) (imag-part a))))]
+    (complex (/ (real-part a) denom) (/ (- (imag-part a)) denom))))
+
+(define-method (/ (a <hypercomplex<>>) (b <hypercomplex<>>))
+  (jit-let [(inverse (hypercomplex-inverse b))] (* a inverse)))
+(define-method (/ (a <hypercomplex<>>) (b <complex<>>))
+  (jit-let [(inverse (complex-inverse b))] (* a inverse)))
+(define-method (/ (a <complex<>>) (b <hypercomplex<>>))
+  (jit-let [(inverse (hypercomplex-inverse b))] (* a inverse)))

etc/playground.scm

@@ -1,5 +1,7 @@
 (use-modules (srfi srfi-1) (oop goops) (aiscm core) (aiscm util) (aiscm hypercomplex))
 
+; // http://mathworld.wolfram.com/HypercomplexNumber.html
+
 (define (hypercomplex-inverse a)
   (jit-let [(d1 (- (real-part a) (kmag-part a)))
             (d2 (+ (imag-part a) (jmag-part a)))

tests/test_hypercomplex.scm

@@ -166,6 +166,23 @@
   (test-equal "multiply hypercomplex and real number"
     (make-hypercomplex 4.0 6.0 10.0 14.0) ((jit (list (hypercomplex <float>) <float>) *) (make-hypercomplex 2 3 5 7) 2)))
 
+(test-group "hypercomplex division"
+  (let [(f (jit (list (hypercomplex <float>) (hypercomplex <float>)) /))]
+    (test-approximate "real part is one when dividing number by itself"
+      1.0 (real-part (f (make-hypercomplex 2 3 5 7) (make-hypercomplex 2 3 5 7))) 0.001)
+    (test-approximate "imaginary part is zero when dividing number by itself"
+      0.0 (imag-part (f (make-hypercomplex 2 3 5 7) (make-hypercomplex 2 3 5 7))) 0.001)
+    (test-approximate "jmaginary part is zero when dividing number by itself"
+      0.0 (jmag-part (f (make-hypercomplex 2 3 5 7) (make-hypercomplex 2 3 5 7))) 0.001)
+    (test-approximate "kmaginary part is zero when dividing number by itself"
+      0.0 (kmag-part (f (make-hypercomplex 2 3 5 7) (make-hypercomplex 2 3 5 7))) 0.001))
+  (let [(f (jit (list (hypercomplex <float>) (complex <float>)) /))]
+    (test-approximate "real part is one for trivial hypercomplex-complex division"
+      1.0 (real-part (f (make-hypercomplex 2 3 0 0) 2+3i)) 0.001))
+  (let [(f (jit (list (complex <float>) (hypercomplex <float>)) /))]
+    (test-approximate "real part is one for trivial complex-hypercomplex division"
+      1.0 (real-part (f 2+3i (make-hypercomplex 2 3 0 0))) 0.001)))
+
 (test-group "hypercomplex properties"
   (test-equal "hypercomplex absolute value"
     4.0 ((jit (list (hypercomplex <float>)) abs) (make-hypercomplex 2 2 -2 -2)))