;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; BASIC STORY ON RECURSIVE DEFINITION OF DIVISION ALGEBRAS

;; Building algebras recursively on top of the reals.  
;; Cayley shows how to get reals (R), complex numbers (C), 
;; quaternions (H), octonions (O), etc. using a simple 
;; recursive multiplication rule:

;; 	(a,b) (c,d) = (ac - db*, a*d + cb)

;; for suitable recursive definitions of negation and conjugation.
;; Here the notation c = (a,b) represents a Cayley number c, of
;; order n, composed of two Cayley numbers (a and b), of order (n-1).
;; The hierarchy bottoms out with the reals (Cayley number 0).

;; This works for any n, although the algebras constructed
;; get weirder and weirder. To quote:

;;	There are exactly four normed division algebras: 
;;	the real numbers (R), complex numbers (C), 
;;	quaternions (H), and octonions (O). 
;;	The real numbers are the dependable breadwinner of the family, 
;;	the complete ordered field we all rely on. 
;;	The complex numbers are a slightly flashier, but still respectable
;;	younger brother: not ordered, but algebraically complete. 
;;	The quaternions, being noncommutative, are the eccentric cousin 
;;	who is shunned at important family gatherings. 
;;	But the	octonions are the crazy old uncle nobody lets out of 
;;	the attic: they are nonassociative.

;; from The Octonions by John C. Baez.  

;; It does go beyond those four, next are the sedenions.
;; But these algebras tend to become somewhat less interesting since 
;; they are then no longer division algebras ---
;; that is, ab=0 no longer implies that either a=0 or b=0.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Constructor (the term "number" here is taken to mean a Cayley number)

(define make-number list)

;; Note: can't use make-number for reals, i.e. Cayley number of order 0.
;; Those are just represented as Scheme numbers for convenience.

;; Selectors (the terms "real" and "imaginary" here are generalized)

(define real car)

(define imaginary cadr)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define z (make-number 3 4))
(real z)
(imaginary z)

(define q (make-number (make-number 1 2) 
		       (make-number 3 4)))
(real q)
(real (real q))

(define o (make-number (make-number (make-number 1 2)
				    (make-number 3 4))
		       (make-number (make-number 5 6)
				    (make-number 7 8))))
(imaginary o)
(real (imaginary o))
(imaginary (real (imaginary o)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Some handy short-hand constructors

(define make-complex make-number)

(define z (make-complex 3 4))

(define (make-quaternion a b c d) 
  (make-number (make-complex a b)
	       (make-complex c d)))

(define q (make-quaternion 1 2 3 4))

(define (make-octonion a b c d e f g h) 
  (make-number (make-quaternion a b c d) 
	       (make-quaternion e f g h)))

(define o (make-octonion 1 2 3 4 5 6 7 8))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Simple operations:

;; Negation -(a, b) = (-a, -b)

(define (negate x)
  (if (number? x) 
      (- x)
      (make-number (negate (real x))
		   (negate (imaginary x)))))

;; Be careful not to violate the abstraction barrier!
;; For example, do not use - instead of negate

(negate z)
(negate q)
(negate o)

;; Conjugation (a, b)* = (a*, -b)
;; Note generalization of conjugation for complex numbers

(define (conjugate x)
  (if (number? x)
      x
      (make-number (conjugate (real x))
		   (negate (imaginary x)))))

;; Be careful not to violate the abstraction barrier!
;; For example, do not use - instead of negate

(conjugate z)
(conjugate q)
(conjugate o)

;; Addition

(define (add x y)
  (if (and (number? x) (number? y))
      (+ x y)
      (make-number (add (real x) (real y))
		   (add (imaginary x) (imaginary y)))))

;; Be careful not to violate the abstraction barrier!
;; For example, do not use + instead of add

(add q q)

;; Subtraction

(define (subtract x y)
  (if (and (number? x) (number? y))
      (- x y)
      (make-number (subtract (real x) (real y))
		   (subtract (imaginary x) (imaginary y)))))

(subtract o o)

;; Scale a number (second argument is a real scale factor)

(define (scale x k)
  (if (number? x)
      (* x k)
      (make-number (scale (real x) k)
		   (scale (imaginary x) k))))

(scale o 0.5)

;; Multiplication --- this is the heart of the algebra!
;; Be careful about the order of multiplicant and multiplier
;; since multiplication is not guaranteed to be commutative!

;; 	(a,b) (c,d) = (ac - db*, a*d + cb)

(define (multiply x y)
  (if (and (number? x) (number? y))
      (* x y)
      (make-number
       (subtract (multiply (real x) (real y))
		 (multiply (imaginary y) (conjugate (imaginary x))))
       (add (multiply (conjugate (real x)) (imaginary y))
	    (multiply (real y) (imaginary x))))))

(multiply q (conjugate q))

(multiply o (conjugate o))

;; Real part of real part of real part ...

(define (real-part x)
  (if (number? x)
      x
      (real-part (real x))))

(real-part (multiply o (conjugate o)))

;; Magnitude squared ||a||^2 = a a* = a* a

(define (magnitude-squared x)
  (if (number? x)
      (* x x)
      (real-part (multiply x (conjugate x)))))

(magnitude-squared q)

(define (magnitude x) (sqrt (magnitude-squared x)))

(magnitude z)

;; In division algebras numbers (other than zero) have inverses 
;; Multiplicative inverse a^{-1} = a* / ||a||^2, 

(define (inverse x)
  (if (number? x) 
      (if (zero? x)
	  (error "ERROR: division by zero")
	  (/ 1 x))
      (let ((magsq (magnitude-squared x)))
	(if (zero? magsq)
	    (error "ERROR: division by zero")
	    (scale (conjugate x) (/ 1 magsq))))))

(inverse z)

(multiply z (inverse z))

(inverse q)

(multiply q (inverse q))

(inverse o)

(multiply o (inverse o))

;; END OF BASIC STORY ON DIVISION ALGEBRAS

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Some extra stuff that may be useful ...

;; Make Cayley number zero of order n
;; Try and be efficient by sharing

(define (make-zero n)
  (if (zero? n)
      0
      (let ((lower (make-zero (- n 1))))
	(make-number lower lower))))

(make-zero 4)

;; Make Cayley number one of order n

(define (make-one n)
  (if (zero? n)
      1
      (make-number (make-one (- n 1))
		   (make-zero (- n 1)))))

(make-one 4)

;; Check whether equal to zero

(define (is-zero? x)
  (if (number? x)
      (zero? x)
      (and (is-zero? (real x))
	   (is-zero? (imaginary x)))))

(is-zero? (make-zero 4))

;; Check whether equal to one

(define (is-one? x)
  (if (number? x)
      (= x 1)
      (and (is-one? (real x))
	   (is-zero? (imaginary x)))))

(is-one? (make-one 4))

;; Normalize Cayley number to unit magnitude

(define (normalize x)
  (if (number? x)
      (if (zero? x)
	  (error "ERROR: division by zero")
	  1 )
      (let ((mag (magnitude x)))
	(if (zero? mag)
	    (error "ERROR: division by zero")
	    (scale x (/ 1 mag))))))

(magnitude (normalize o))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Quaternions and higher do not commute

;; Useful test procedure
;; Commutator [u; v] = u * v - v * u
;; Should be zero for order n < 2 (real and complex)

(define (commutator u v)
  (subtract (multiply u v)
	    (multiply v u)))

(define z1 (make-complex 1 2))
(define z2 (make-complex 5 -1))

(commutator z1 z2)

(define q1 (make-quaternion 1 2 3 4))
(define q2 (make-quaternion 8 -3 7 5))

(commutator q1 q2)
(commutator q1 q1)

;; Octonions and higher are not associative

;; Useful test procedure
;; Associator [u; v; w] = (uv)w - u(vw)
;; Should be zero for n < 3 (real, complex, and quaternion)

(define (associator u v w)
  (subtract (multiply (multiply u v) w)
	    (multiply u (multiply v w))))

(define z3 (make-complex 3 1))

(associator z1 z2 z3)

(define q3 (make-quaternion -1 2 9 13))

(associator q1 q2 q3)

(define o1 (make-octonion 1 2 3 4 5 6 7 8))
(define o2 (make-octonion 6 5 -8 7 2 -1 4 3))
(define o3 (make-octonion 6 -15 8 17 -2 11 4 -13))

(associator o1 o2 o3)
(associator o1 o2 o2)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; DANGER: Leaving the world of division algebras...

(define (make-sedenion a b c d e f g h i j k l m n o p)
  (make-number (make-octonion a b c d e f g h)
	       (make-octonion i j k l m n o p)))

(define s1 (make-sedenion 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6))

(define s2 (make-sedenion 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0))

(define s3 (make-sedenion 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0))

(multiply s1 s2)

(commutator s1 s2)

(associator s1 s2 s3)

;; In abstract algebra, an algebra is called alternative if the 
;; subalgebra generated by any two of its elements is associative.

;; That is if x(xy) = (xx)y or equivalently (xy)y = x(yy)

;; Useful test function

(define (alternative x y)
  (subtract (multiply x (multiply x y))
	    (multiply (multiply x x) y)))

(alternative z1 z2)
(alternative q1 q2)
(alternative o1 o2)
(alternative s1 s2)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; EXTENSIONS: presently this assumes binary operators
;; are given Cayley numbers of the same order as arguments
;; Extend to allow mixed arguments by "coercing" or "promoting"
;; the lower order Cayley to the order of the higher.

;; Get the order of a Cayley number --- assumed "balanced"
;; (that is, the real and imaginary parts have the same order)

(define (order x)
  (if (number? x) 
      0
      (+ 1 (order (real x)))))

(order 4)
(order z1)
(order q1)
(order o1)
(order s1)

;; "Coerce" a to be same order Cayley number as b
;; assume order of b is equal to or greater than that of a

(define (coerce a b)
  (if (= (order a) (order b)) 
      a
      (coerce (make-number a (coerce 0 a)) b)))

(coerce 4 z1)
(coerce 4 q1)
(coerce 4 o1)
(coerce 4 s1)

(coerce z1 4)

;; Now could modify all arithmetic procedure taking arguments a and b
;; to instead us (coerce a b) and (coerce b a)

;; Further, allow for unbalanced Cayley numbers...

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; For application of unit quaternions to rotation in 3D, see
;; Quaternions_for_Rotation_in_Scheme.scm

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;