mjo/octonions.py

   1 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
   2 from sage.combinat.free_module import CombinatorialFreeModule
   3 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
   4 from sage.categories.magmatic_algebras import MagmaticAlgebras
   5 from sage.rings.all import AA, ZZ
   6 from sage.matrix.matrix_space import MatrixSpace
   7 from sage.misc.table import table
   8
   9 class Octonion(IndexedFreeModuleElement):
  10     def conjugate(self):
  11         r"""
  12         SETUP::
  13
  14             sage: from mjo.octonions import Octonions
  15
  16         EXAMPLES::
  17
  18             sage: O = Octonions()
  19             sage: x = sum(O.gens())
  20             sage: x.conjugate()
  21             e0 - e1 - e2 - e3 - e4 - e5 - e6 - e7
  22
  23         TESTS::
  24
  25         Conjugating twice gets you the original element::
  26
  27             sage: set_random_seed()
  28             sage: O = Octonions()
  29             sage: x = O.random_element()
  30             sage: x.conjugate().conjugate() == x
  31             True
  32
  33         """
  34         C = MatrixSpace(ZZ,8).diagonal_matrix((1,-1,-1,-1,-1,-1,-1,-1))
  35         return self.parent().from_vector(C*self.to_vector())
  36
  37     def real(self):
  38         r"""
  39         Return the real part of this octonion.
  40
  41         The real part of an octonion is its projection onto the span
  42         of the first generator. In other words, the "first dimension"
  43         is real and the others are imaginary.
  44
  45         SETUP::
  46
  47             sage: from mjo.octonions import Octonions
  48
  49         EXAMPLES::
  50
  51             sage: O = Octonions()
  52             sage: x = sum(O.gens())
  53             sage: x.real()
  54             e0
  55
  56         TESTS:
  57
  58         This method is idempotent::
  59
  60             sage: set_random_seed()
  61             sage: O = Octonions()
  62             sage: x = O.random_element()
  63             sage: x.real().real() == x.real()
  64             True
  65
  66         """
  67         return (self + self.conjugate())/2
  68
  69     def imag(self):
  70         r"""
  71         Return the imaginary part of this octonion.
  72
  73         The imaginary part of an octonion is its projection onto the
  74         orthogonal complement of the span of the first generator. In
  75         other words, the "first dimension" is real and the others are
  76         imaginary.
  77
  78         SETUP::
  79
  80             sage: from mjo.octonions import Octonions
  81
  82         EXAMPLES::
  83
  84             sage: O = Octonions()
  85             sage: x = sum(O.gens())
  86             sage: x.imag()
  87             e1 + e2 + e3 + e4 + e5 + e6 + e7
  88
  89         TESTS:
  90
  91         This method is idempotent::
  92
  93             sage: set_random_seed()
  94             sage: O = Octonions()
  95             sage: x = O.random_element()
  96             sage: x.imag().imag() == x.imag()
  97             True
  98
  99         """
 100         return (self - self.conjugate())/2
 101
 102     def _norm_squared(self):
 103         return (self*self.conjugate()).coefficient(0)
 104
 105     def norm(self):
 106         r"""
 107         Return the norm of this octonion.
 108
 109         SETUP::
 110
 111             sage: from mjo.octonions import Octonions
 112
 113         EXAMPLES::
 114
 115             sage: O = Octonions()
 116             sage: O.one().norm()
 117             1
 118
 119         TESTS:
 120
 121         The norm is nonnegative and belongs to the base field::
 122
 123             sage: set_random_seed()
 124             sage: O = Octonions()
 125             sage: n = O.random_element().norm()
 126             sage: n >= 0 and n in O.base_ring()
 127             True
 128
 129         The norm is homogeneous::
 130
 131             sage: set_random_seed()
 132             sage: O = Octonions()
 133             sage: x = O.random_element()
 134             sage: alpha = O.base_ring().random_element()
 135             sage: (alpha*x).norm() == alpha.abs()*x.norm()
 136             True
 137
 138         """
 139         return self._norm_squared().sqrt()
 140
 141     def inverse(self):
 142         r"""
 143         Return the inverse of this element if it exists.
 144
 145         SETUP::
 146
 147             sage: from mjo.octonions import Octonions
 148
 149         EXAMPLES::
 150
 151             sage: O = Octonions()
 152             sage: x = sum(O.gens())
 153             sage: x*x.inverse() == O.one()
 154             True
 155
 156         ::
 157
 158             sage: O = Octonions()
 159             sage: O.one().inverse() == O.one()
 160             True
 161
 162         TESTS::
 163
 164             sage: set_random_seed()
 165             sage: O = Octonions()
 166             sage: x = O.random_element()
 167             sage: x.is_zero() or ( x*x.inverse() == O.one() )
 168             True
 169
 170         """
 171         if self.is_zero():
 172             raise ValueError("zero is not invertible")
 173         return self.conjugate()/self._norm_squared()
 174
 175
 176     def cayley_dickson(self, Q=None):
 177         r"""
 178         Return the Cayley-Dickson representation of this element in terms
 179         of the quaternion algebra ``Q``.
 180
 181         The Cayley-Dickson representation is an identification of
 182         octionions `x` and `y` with pairs of quaternions `(a,b)` and
 183         `(c,d)` respectively such that:
 184
 185         * `x + y = (a+b, c+d)`
 186         * `xy` = (ac - \bar{d}*b, da + b\bar{c})`
 187         * `\bar{x} = (a,-b)`
 188
 189         where `\bar{x}` denotes the conjugate of `x`.
 190
 191         SETUP::
 192
 193             sage: from mjo.octonions import Octonions
 194
 195         EXAMPLES::
 196
 197             sage: O = Octonions()
 198             sage: x = sum(O.gens())
 199             sage: x.cayley_dickson()
 200             (1 + i + j + k, 1 + i + j + k)
 201
 202         """
 203         if Q is None:
 204             Q = QuaternionAlgebra(self.base_ring(), -1, -1)
 205
 206         i,j,k = Q.gens()
 207         a = (self.coefficient(0)*Q.one() +
 208              self.coefficient(1)*i +
 209              self.coefficient(2)*j +
 210              self.coefficient(3)*k )
 211         b = (self.coefficient(4)*Q.one() +
 212              self.coefficient(5)*i +
 213              self.coefficient(6)*j +
 214              self.coefficient(7)*k )
 215
 216         from sage.categories.sets_cat import cartesian_product
 217         P = cartesian_product([Q,Q])
 218         return P((a,b))
 219
 220
 221 class Octonions(CombinatorialFreeModule):
 222     r"""
 223     SETUP::
 224
 225         sage: from mjo.octonions import Octonions
 226
 227     EXAMPLES::
 228
 229         sage: Octonions()
 230         Octonion algebra with base ring Algebraic Real Field
 231         sage: Octonions(field=QQ)
 232         Octonion algebra with base ring Rational Field
 233
 234     """
 235     def __init__(self,
 236                  field=AA,
 237                  prefix="e"):
 238
 239         # Not associative, not commutative
 240         category = MagmaticAlgebras(field).FiniteDimensional()
 241         category = category.WithBasis().Unital()
 242
 243         super().__init__(field,
 244                          range(8),
 245                          element_class=Octonion,
 246                          category=category,
 247                          prefix=prefix,
 248                          bracket=False)
 249
 250         # The product of each basis element is plus/minus another
 251         # basis element that can simply be looked up on
 252         # https://en.wikipedia.org/wiki/Octonion
 253         e0, e1, e2, e3, e4, e5, e6, e7 = self.gens()
 254         self._multiplication_table = (
 255             (e0, e1, e2, e3, e4, e5, e6, e7),
 256             (e1,-e0, e3,-e2, e5,-e4,-e7, e6),
 257             (e2,-e3,-e0, e1, e6, e7,-e4,-e5),
 258             (e3, e2,-e1,-e0, e7,-e6, e5,-e4),
 259             (e4,-e5,-e6,-e7,-e0, e1, e2, e3),
 260             (e5, e4,-e7, e6,-e1,-e0,-e3, e2),
 261             (e6, e7, e4,-e5,-e2, e3,-e0,-e1),
 262             (e7,-e6, e5, e4,-e3,-e2, e1,-e0),
 263         )
 264
 265     def product_on_basis(self, i, j):
 266         return self._multiplication_table[i][j]
 267
 268     def one_basis(self):
 269         r"""
 270         Return the monomial index (basis element) corresponding to the
 271         octonion unit element.
 272
 273         SETUP::
 274
 275             sage: from mjo.octonions import Octonions
 276
 277         TESTS:
 278
 279         This gives the correct unit element::
 280
 281             sage: set_random_seed()
 282             sage: O = Octonions()
 283             sage: x = O.random_element()
 284             sage: x*O.one() == x and O.one()*x == x
 285             True
 286
 287         """
 288         return 0
 289
 290     def _repr_(self):
 291         return ("Octonion algebra with base ring %s" % self.base_ring())
 292
 293     def multiplication_table(self):
 294         """
 295         Return a visual representation of this algebra's multiplication
 296         table (on basis elements).
 297
 298         SETUP::
 299
 300             sage: from mjo.octonions import Octonions
 301
 302         EXAMPLES:
 303
 304         The multiplication table is what Wikipedia says it is::
 305
 306             sage: Octonions().multiplication_table()
 307             +----++----+-----+-----+-----+-----+-----+-----+-----+
 308             | *  || e0 | e1  | e2  | e3  | e4  | e5  | e6  | e7  |
 309             +====++====+=====+=====+=====+=====+=====+=====+=====+
 310             | e0 || e0 | e1  | e2  | e3  | e4  | e5  | e6  | e7  |
 311             +----++----+-----+-----+-----+-----+-----+-----+-----+
 312             | e1 || e1 | -e0 | e3  | -e2 | e5  | -e4 | -e7 | e6  |
 313             +----++----+-----+-----+-----+-----+-----+-----+-----+
 314             | e2 || e2 | -e3 | -e0 | e1  | e6  | e7  | -e4 | -e5 |
 315             +----++----+-----+-----+-----+-----+-----+-----+-----+
 316             | e3 || e3 | e2  | -e1 | -e0 | e7  | -e6 | e5  | -e4 |
 317             +----++----+-----+-----+-----+-----+-----+-----+-----+
 318             | e4 || e4 | -e5 | -e6 | -e7 | -e0 | e1  | e2  | e3  |
 319             +----++----+-----+-----+-----+-----+-----+-----+-----+
 320             | e5 || e5 | e4  | -e7 | e6  | -e1 | -e0 | -e3 | e2  |
 321             +----++----+-----+-----+-----+-----+-----+-----+-----+
 322             | e6 || e6 | e7  | e4  | -e5 | -e2 | e3  | -e0 | -e1 |
 323             +----++----+-----+-----+-----+-----+-----+-----+-----+
 324             | e7 || e7 | -e6 | e5  | e4  | -e3 | -e2 | e1  | -e0 |
 325             +----++----+-----+-----+-----+-----+-----+-----+-----+
 326
 327         """
 328         n = self.dimension()
 329         # Prepend the header row.
 330         M = [["*"] + list(self.gens())]
 331
 332         # And to each subsequent row, prepend an entry that belongs to
 333         # the left-side "header column."
 334         M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
 335                                     for j in range(n) ]
 336                for i in range(n) ]
 337
 338         return table(M, header_row=True, header_column=True, frame=True)
 339
 340
 341 class OctonionMatrix:
 342     r"""
 343     A pseudo-matrix class that supports octonion entries.
 344
 345     Matrices in SageMath can't have base rings that are
 346     non-commutative, much less non-associative. The "matrix" scaling,
 347     addition, and multiplication operations for this class are all
 348     wholly inefficient, but are hand-written to guarantee that they
 349     are performed in the correct order. Of course, it can't guarantee
 350     that you won't write something visually ambiguous like
 351     `A*B*C`... but you already have that problem with the
 352     non-associative octonions themselves.
 353
 354     This class is only as sophisticated as it need to be to implement
 355     the Jordan and inner-products in the space of Hermitian matrices
 356     with octonion entries, namely ``(X*Y+Y*X)/2`` and
 357     ``(X*Y).trace().real()`` for two octonion matrices ``X`` and
 358     ``Y``.
 359
 360     .. WARNING:
 361
 362         These are not matrices in the usual sense! Matrix
 363         multiplication is associative. Multiplication of octonion
 364         "matrices" cannot be, since the multiplication of the
 365         underlying octonions is not (consider two 1-by-1 matrices each
 366         containing a single octonion).
 367     """
 368     def __init__(self, entries):
 369         r"""
 370         Initialize this matrix with a list of lists in (row,column) order,
 371         just like in SageMath.
 372         """
 373         self._nrows = len(entries)
 374
 375         if self._nrows == 0:
 376             self._ncols = 0
 377         else:
 378             # We don't check that you haven't supplied two rows (or
 379             # columns) of different lengths!
 380             self._ncols = len(entries[0])
 381
 382         self._entries = entries
 383
 384     def __getitem__(self, indices):
 385         r"""
 386         SETUP::
 387
 388             sage: from mjo.octonions import Octonions, OctonionMatrix
 389
 390         EXAMPLES::
 391
 392             sage: O = Octonions(field=QQ)
 393             sage: M = OctonionMatrix([ [O.one(),  O.zero()],
 394             ....:                      [O.zero(), O.one() ] ])
 395             sage: M[0,0]
 396             e0
 397             sage: M[0,1]
 398             0
 399             sage: M[1,0]
 400             0
 401             sage: M[1,1]
 402             e0
 403         """
 404         i,j = indices
 405         return self._entries[i][j]
 406
 407     def nrows(self):
 408         r"""
 409         SETUP::
 410
 411             sage: from mjo.octonions import Octonions, OctonionMatrix
 412
 413         EXAMPLES::
 414
 415             sage: O = Octonions(field=QQ)
 416             sage: M = OctonionMatrix([ [O.one(),  O.zero()],
 417             ....:                      [O.zero(), O.one() ],
 418             ....:                      [O.zero(), O.zero()] ])
 419             sage: M.nrows()
 420             3
 421
 422         """
 423         return self._nrows
 424
 425     def ncols(self):
 426         r"""
 427
 428         SETUP::
 429
 430             sage: from mjo.octonions import Octonions, OctonionMatrix
 431
 432         EXAMPLES::
 433
 434             sage: O = Octonions(field=QQ)
 435             sage: M = OctonionMatrix([ [O.one(),  O.zero()],
 436             ....:                      [O.zero(), O.one() ],
 437             ....:                      [O.zero(), O.zero()] ])
 438             sage: M.ncols()
 439             2
 440
 441         """
 442         return self._ncols
 443
 444     def __repr__(self):
 445         return table(self._entries, frame=True)._repr_()
 446
 447     def __mul__(self,rhs):
 448         r"""
 449
 450         SETUP::
 451
 452             sage: from mjo.octonions import Octonions, OctonionMatrix
 453
 454         EXAMPLES::
 455
 456             sage: O = Octonions(QQ)
 457             sage: e1 = O.monomial(1)
 458             sage: e2 = O.monomial(2)
 459             sage: e1*e2
 460             e3
 461             sage: e2*e1
 462             -e3
 463             sage: E1 = OctonionMatrix([[e1]])
 464             sage: E2 = OctonionMatrix([[e2]])
 465             sage: E1*E2
 466             +----+
 467             | e3 |
 468             +----+
 469             sage: E2*E1
 470             +-----+
 471             | -e3 |
 472             +-----+
 473
 474         """
 475         if not self.ncols() == rhs.nrows():
 476             raise ValueError("dimension mismatch")
 477
 478         m = self.nrows()
 479         n = self.ncols()
 480         p = rhs.ncols()
 481
 482         C = lambda i,j: sum( self[i,k]*rhs[k,j] for k in range(n) )
 483         return OctonionMatrix([ [C(i,j) for j in range(m)]
 484                                 for i in range(p) ] )
 485
 486     def __rmul__(self,scalar):
 487         r"""
 488
 489         SETUP::
 490
 491             sage: from mjo.octonions import Octonions, OctonionMatrix
 492
 493         EXAMPLES::
 494
 495             sage: O = Octonions(QQ)
 496             sage: M = OctonionMatrix([[O.one(), O.zero()],
 497             ....:                     [O.zero(),O.one() ] ])
 498             sage: 2*M
 499             +------+------+
 500             | 2*e0 | 0    |
 501             +------+------+
 502             | 0    | 2*e0 |
 503             +------+------+
 504
 505         r"""
 506         # SCALAR GOES ON THE LEFT HERE
 507         return OctonionMatrix([ [scalar*self[i,j]
 508                                  for j in range(self.ncols())]
 509                                 for i in range(self.nrows()) ])
 510
 511     def __add__(self,rhs):
 512         r"""
 513         SETUP::
 514
 515             sage: from mjo.octonions import Octonions, OctonionMatrix
 516
 517         EXAMPLES::
 518
 519             sage: O = Octonions(QQ)
 520             sage: e0,e1,e2,e3,e4,e5,e6,e7 = O.gens()
 521             sage: A = OctonionMatrix([ [e0,e1],
 522             ....:                      [e2,e3] ])
 523             sage: B = OctonionMatrix([ [e4,e5],
 524             ....:                      [e6,e7] ])
 525             sage: A+B
 526             +---------+---------+
 527             | e0 + e4 | e1 + e5 |
 528             +---------+---------+
 529             | e2 + e6 | e3 + e7 |
 530             +---------+---------+
 531
 532         """
 533         if not self.ncols() == rhs.ncols():
 534             raise ValueError("column dimension mismatch")
 535         if not self.nrows() == rhs.nrows():
 536             raise ValueError("row dimension mismatch")
 537         return OctonionMatrix([ [self[i,j] + rhs[i,j]
 538                                  for j in range(self.ncols())]
 539                                 for i in range(self.nrows()) ])