--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2009, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
)abbrev category QUATCAT QuaternionCategory
++ Author: Robert S. Sutor
++ Date Created: 23 May 1990
++ Change History:
++   10 September 1990
++ Basic Operations: (Algebra)
++   abs, conjugate, imagI, imagJ, imagK, norm, quatern, rational,
++   rational?, real
++ Related Constructors: Quaternion, QuaternionCategoryFunctions2
++ Also See: DivisionRing
++ AMS Classifications: 11R52
++ Keywords: quaternions, division ring, algebra
++ Description:
++   \spadtype{QuaternionCategory} describes the category of quaternions
++   and implements functions that are not representation specific.
 
QuaternionCategory(R: CommutativeRing): Category ==
  Join(Algebra R, FullyRetractableTo R, DifferentialExtension R,
   FullyEvalableOver R, FullyLinearlyExplicitRingOver R) with
 
     conjugate: $ -> $
       ++ conjugate(q) negates the imaginary parts of quaternion \spad{q}.
     imagI:   $ -> R
       ++ imagI(q) extracts the imaginary i part of quaternion \spad{q}.
     imagJ:   $ -> R
       ++ imagJ(q) extracts the imaginary j part of quaternion \spad{q}.
     imagK:   $ -> R
       ++ imagK(q) extracts the imaginary k part of quaternion \spad{q}.
     norm:    $ -> R
       ++ norm(q) computes the norm of \spad{q} (the sum of the
       ++ squares of the components).
     quatern: (R,R,R,R) -> $
       ++ quatern(r,i,j,k) constructs a quaternion from scalars.
     real:    $ -> R
       ++ real(q) extracts the real part of quaternion \spad{q}.
 
     if R has EntireRing then EntireRing
     if R has OrderedSet then OrderedSet
     if R has Field then DivisionRing
     if R has ConvertibleTo InputForm then ConvertibleTo InputForm
     if R has CharacteristicZero then CharacteristicZero
     if R has CharacteristicNonZero then CharacteristicNonZero
     if R has RealNumberSystem then
       abs    : $ -> R
         ++ abs(q) computes the absolute value of quaternion \spad{q}
         ++ (sqrt of norm).
     if R has IntegerNumberSystem then
       rational?    : $ -> Boolean
         ++ rational?(q) returns {\it true} if all the imaginary
         ++ parts of \spad{q} are zero and the real part can be
         ++ converted into a rational number, and {\it false}
         ++ otherwise.
       rational     : $ -> Fraction Integer
         ++ rational(q) tries to convert \spad{q} into a
         ++ rational number. Error: if this is not
         ++ possible. If \spad{rational?(q)} is true, the
         ++ conversion will be done and the rational number returned.
       rationalIfCan: $ -> Union(Fraction Integer, "failed")
         ++ rationalIfCan(q) returns \spad{q} as a rational number,
         ++ or "failed" if this is not possible.
         ++ Note: if \spad{rational?(q)} is true, the conversion
         ++ can be done and the rational number will be returned.
 
 add
 
       characteristic ==
         characteristic$R
       conjugate x      ==
         quatern(real x, - imagI x, - imagJ x, - imagK x)
       map(fn, x)       ==
         quatern(fn real x, fn imagI x, fn imagJ x, fn imagK x)
       norm x ==
         real x * real x + imagI x * imagI x +
           imagJ x * imagJ x + imagK x * imagK x
       x = y            ==
         (real x = real y) and (imagI x = imagI y) and
           (imagJ x = imagJ y) and (imagK x = imagK y)
       x + y            ==
         quatern(real x + real y, imagI x + imagI y,
           imagJ x + imagJ y, imagK x + imagK y)
       x - y            ==
         quatern(real x - real y, imagI x - imagI y,
           imagJ x - imagJ y, imagK x - imagK y)
       - x              ==
         quatern(- real x, - imagI x, - imagJ x, - imagK x)
       r:R * x:$        ==
         quatern(r * real x, r * imagI x, r * imagJ x, r * imagK x)
       n:Integer * x:$  ==
         quatern(n * real x, n * imagI x, n * imagJ x, n * imagK x)
       differentiate(x:$, d:R -> R) ==
         quatern(d real x, d imagI x, d imagJ x, d imagK x)
       coerce(r:R)      ==
         quatern(r,0$R,0$R,0$R)
       coerce(n:Integer)      ==
         quatern(n :: R,0$R,0$R,0$R)
       one? x ==
         one? real x and zero? imagI x and
           zero? imagJ x and zero? imagK x
       zero? x ==
         zero? real x and zero? imagI x and
           zero? imagJ x and zero? imagK x
       retract(x):R ==
         not (zero? imagI x and zero? imagJ x and zero? imagK x) =>
           error "Cannot retract quaternion."
         real x
       retractIfCan(x):Union(R,"failed") ==
         not (zero? imagI x and zero? imagJ x and zero? imagK x) =>
           "failed"
         real x
 
       coerce(x:$):OutputForm ==
         part,z : OutputForm
         y : $
         zero? x => (0$R) :: OutputForm
         not zero?(real x) =>
           y := quatern(0$R,imagI(x),imagJ(x),imagK(x))
           zero? y => real(x) :: OutputForm
           (real(x) :: OutputForm) + (y :: OutputForm)
         -- we know that the real part is 0
         not zero?(imagI(x)) =>
           y := quatern(0$R,0$R,imagJ(x),imagK(x))
           z :=
             part := 'i::OutputForm
             one? imagI(x) => part
             (imagI(x) :: OutputForm) * part
           zero? y => z
           z + (y :: OutputForm)
         -- we know that the real part and i part are 0
         not zero?(imagJ(x)) =>
           y := quatern(0$R,0$R,0$R,imagK(x))
           z :=
             part := 'j::OutputForm
             one? imagJ(x) => part
             (imagJ(x) :: OutputForm) * part
           zero? y => z
           z + (y :: OutputForm)
         -- we know that the real part and i and j parts are 0
         part := 'k::OutputForm
         one? imagK(x) => part
         (imagK(x) :: OutputForm) * part
 
       if R has Field then
         inv x ==
           norm x = 0 => error "This quaternion is not invertible."
           (inv norm x) * conjugate x
 
       if R has ConvertibleTo InputForm then
         convert(x:$):InputForm ==
           l : List InputForm := [convert("quatern" :: Symbol),
             convert(real x)$R, convert(imagI x)$R, convert(imagJ x)$R,
               convert(imagK x)$R]
           convert(l)$InputForm
 
       if R has OrderedSet then
         x < y ==
           real x = real y =>
             imagI x = imagI y =>
               imagJ x = imagJ y =>
                 imagK x < imagK y
               imagJ x < imagJ y
             imagI x < imagI y
           real x < real y
 
       if R has RealNumberSystem then
         abs x == sqrt norm x
 
       if R has IntegerNumberSystem then
         rational? x ==
           (zero? imagI x) and (zero? imagJ x) and (zero? imagK x)
         rational  x ==
           rational? x => rational real x
           error "Not a rational number"
         rationalIfCan x ==
           rational? x => rational real x
           "failed"

)abbrev domain QUAT Quaternion
++ Author: Robert S. Sutor
++ Date Created: 23 May 1990
++ Change History:
++   10 September 1990
++ Basic Operations: (Algebra)
++   abs, conjugate, imagI, imagJ, imagK, norm, quatern, rational,
++   rational?, real
++ Related Constructors: QuaternionCategoryFunctions2
++ Also See: QuaternionCategory, DivisionRing
++ AMS Classifications: 11R52
++ Keywords: quaternions, division ring, algebra
++ Description: \spadtype{Quaternion} implements quaternions over a
++   commutative ring. The main constructor function is \spadfun{quatern}
++   which takes 4 arguments: the real part, the i imaginary part, the j
++   imaginary part and the k imaginary part.
 
Quaternion(R:CommutativeRing): QuaternionCategory(R) == add
  Rep := Record(r:R,i:R,j:R,k:R)
 
  0 == [0,0,0,0]
  1 == [1,0,0,0]
 
  a,b,c,d : R
  x,y : $
 
  real  x == x.r
  imagI x == x.i
  imagJ x == x.j
  imagK x == x.k
 
  quatern(a,b,c,d) == [a,b,c,d]
 
  x * y == [x.r*y.r-x.i*y.i-x.j*y.j-x.k*y.k,
               x.r*y.i+x.i*y.r+x.j*y.k-x.k*y.j,
                 x.r*y.j+x.j*y.r+x.k*y.i-x.i*y.k,
                   x.r*y.k+x.k*y.r+x.i*y.j-x.j*y.i]

)abbrev package QUATCT2 QuaternionCategoryFunctions2
++ Author: Robert S. Sutor
++ Date Created: 23 May 1990
++ Change History:
++   23 May 1990
++ Basic Operations: map
++ Related Constructors: QuaternionCategory, Quaternion
++ Also See:
++ AMS Classifications: 11R52
++ Keywords: quaternions, division ring, map
++ Description:
++    \spadtype{QuaternionCategoryFunctions2} implements functions between
++    two quaternion domains.  The function \spadfun{map} is used by
++    the system interpreter to coerce between quaternion types.
 
QuaternionCategoryFunctions2(QR,R,QS,S) : Exports ==
  Implementation where
    R  : CommutativeRing
    S  : CommutativeRing
    QR : QuaternionCategory R
    QS : QuaternionCategory S
    Exports == with
      map:     (R -> S, QR) -> QS
        ++ map(f,u) maps f onto the component parts of the quaternion
        ++ u.
    Implementation == add
      map(fn : R -> S, u : QR): QS ==
        quatern(fn real u, fn imagI u, fn imagJ u, fn imagK u)$QS

