Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/algebras/quatalg/quaternion_algebra_element.pyx
4156 views
1
"""
2
Elements of Quaternion Algebras
3
4
Sage allows for computation with elements of quaternion algebras over
5
a nearly arbitrary base field of characteristic not 2. Sage also has
6
very highly optimized implementation of arithmetic in rational
7
quaternion algebras and quaternion algebras over number fields.
8
"""
9
10
#*****************************************************************************
11
# Copyright (C) 2009 William Stein <[email protected]>
12
# Copyright (C) 2009 Jonathon Bober <[email protected]>
13
#
14
# Distributed under the terms of the GNU General Public License (GPL)
15
#
16
# This code is distributed in the hope that it will be useful,
17
# but WITHOUT ANY WARRANTY; without even the implied warranty of
18
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19
# General Public License for more details.
20
#
21
# The full text of the GPL is available at:
22
#
23
# http://www.gnu.org/licenses/
24
#*****************************************************************************
25
26
from sage.structure.element cimport AlgebraElement, RingElement, ModuleElement, Element
27
from sage.algebras.quatalg.quaternion_algebra_element cimport QuaternionAlgebraElement_abstract
28
from sage.rings.rational cimport Rational
29
from sage.rings.integer cimport Integer
30
from sage.rings.arith import lcm
31
from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint
32
from sage.rings.number_field.number_field_element cimport NumberFieldElement
33
from sage.rings.all import PolynomialRing
34
from sage.matrix.all import matrix
35
36
include "../../ext/gmp.pxi"
37
include "../../ext/stdsage.pxi"
38
39
include "../../libs/flint/fmpz.pxi"
40
include "../../libs/flint/fmpz_poly.pxi"
41
include "../../libs/flint/ntl_interface.pxd"
42
43
# variables for holding temporary values computed in
44
# QuaternionAlgebraElement_rational_field._mul_()
45
cdef mpz_t T1, T2, t3, t4, t5, t6, t7, t8, s1, s2, U1, U2
46
cdef fmpz_poly_t fT1, fT2, ft3, ft4, ft5, ft6, ft7, ft8, fs1, fs2, fU1, fU2
47
48
def _init_globals():
49
"""
50
Clear all global variables allocated for optimization of
51
quaternion algebra arithmetic.
52
53
Do *not* call this unless you want to leak memory.
54
55
EXAMPLES::
56
57
sage: sage.algebras.quatalg.quaternion_algebra_element._clear_globals()
58
sage: sage.algebras.quatalg.quaternion_algebra_element._init_globals()
59
"""
60
# over QQ
61
mpz_init(T1)
62
mpz_init(T2)
63
mpz_init(t3)
64
mpz_init(t4)
65
mpz_init(t5)
66
mpz_init(t6)
67
mpz_init(t7)
68
mpz_init(t8)
69
70
mpz_init(s1)
71
mpz_init(s2)
72
73
mpz_init(U1)
74
mpz_init(U2)
75
76
# Number fields
77
fmpz_poly_init(fT1)
78
fmpz_poly_init(fT2)
79
fmpz_poly_init(ft3)
80
fmpz_poly_init(ft4)
81
fmpz_poly_init(ft5)
82
fmpz_poly_init(ft6)
83
fmpz_poly_init(ft7)
84
fmpz_poly_init(ft8)
85
86
fmpz_poly_init(fs1)
87
fmpz_poly_init(fs2)
88
89
fmpz_poly_init(fU1)
90
fmpz_poly_init(fU2)
91
92
93
# Initialize module-scope global C variables.
94
_init_globals()
95
96
def _clear_globals():
97
"""
98
Clear all global variables allocated for optimization of
99
quaternion algebra arithmetic.
100
101
Do *not* call this except on exit of Sage, unless you want to see segfaults.
102
103
This is called in the function quit_sage(), which is defined in all.py.
104
105
EXAMPLES::
106
107
sage: sage.algebras.quatalg.quaternion_algebra_element._clear_globals()
108
sage: sage.algebras.quatalg.quaternion_algebra_element._init_globals()
109
"""
110
mpz_clear(T1)
111
mpz_clear(T2)
112
mpz_clear(t3)
113
mpz_clear(t4)
114
mpz_clear(t5)
115
mpz_clear(t6)
116
mpz_clear(t7)
117
mpz_clear(t8)
118
119
mpz_clear(s1)
120
mpz_clear(s2)
121
122
mpz_clear(U1)
123
mpz_clear(U2)
124
125
fmpz_poly_clear(fT1)
126
fmpz_poly_clear(fT2)
127
fmpz_poly_clear(ft3)
128
fmpz_poly_clear(ft4)
129
fmpz_poly_clear(ft5)
130
fmpz_poly_clear(ft6)
131
fmpz_poly_clear(ft7)
132
fmpz_poly_clear(ft8)
133
134
fmpz_poly_clear(fs1)
135
fmpz_poly_clear(fs2)
136
137
fmpz_poly_clear(fU1)
138
fmpz_poly_clear(fU2)
139
140
cdef to_quaternion(R, x):
141
"""
142
Internal function used implicitly by quaternion algebra creation.
143
144
INPUT::
145
146
- R -- callable
147
- x -- element or 4-tuple
148
149
Given a callable R and an x that defines a quaternion, which can be a
150
4-tuple, list of length 4, or something that coerces to R, return
151
4-tuple of elements of R.
152
153
EXAMPLES::
154
sage: Q.<i,j,kkkk> = QuaternionAlgebra(QQ,-7, 13)
155
sage: kkkk._repr_() # implicit doctest
156
'kkkk'
157
"""
158
if isinstance(x, (list, tuple)):
159
return R(x[0]), R(x[1]), R(x[2]), R(x[3])
160
else:
161
return R(x), R(0), R(0), R(0)
162
163
cdef inline print_coeff(y, i, bint atomic):
164
"""
165
Internal function used implicitly by all quaternion algebra printing.
166
167
INPUT::
168
169
- y -- coefficient
170
- i -- string (name of a generator)
171
- atomic -- boolean int; whether or not elements of base ring
172
print atomically
173
174
EXAMPLES::
175
176
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-7, 13)
177
sage: i._repr_() # implicit doctest
178
'i'
179
"""
180
if not y:
181
return ''
182
if y == 1:
183
return i
184
elif y == -1:
185
return "-%s"%i
186
y = str(y)
187
if not atomic and ('+' in y or '-' in y):
188
return '(%s)*%s'%(y, i)
189
else:
190
return '%s*%s'%(y, i)
191
192
cdef class QuaternionAlgebraElement_abstract(AlgebraElement):
193
cpdef bint is_constant(self):
194
"""
195
Return True if this quaternion is constant, i.e., has no i, j, or k term.
196
197
OUTPUT:
198
bool
199
200
EXAMPLES::
201
202
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
203
sage: A(1).is_constant()
204
True
205
sage: A(1+i).is_constant()
206
False
207
sage: A(i).is_constant()
208
False
209
"""
210
return not (self[1] or self[2] or self[3])
211
212
def __int__(self):
213
"""
214
Try to coerce this quaternion to a Python int.
215
216
EXAMPLES::
217
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
218
sage: int(A(-3))
219
-3
220
sage: int(A(-3/2))
221
-2
222
sage: int(-3 + i)
223
Traceback (most recent call last):
224
...
225
TypeError
226
"""
227
if self.is_constant():
228
return int(self[0])
229
raise TypeError
230
231
232
def __long__(self):
233
"""
234
Try to coerce this quaternion to a Python long.
235
236
EXAMPLES::
237
238
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
239
sage: long(A(-3))
240
-3L
241
sage: long(A(-3/2))
242
-2L
243
sage: long(-3 + i)
244
Traceback (most recent call last):
245
...
246
TypeError
247
"""
248
if self.is_constant():
249
return long(self[0])
250
raise TypeError
251
252
def __float__(self):
253
"""
254
Try to coerce this quaternion to a Python float.
255
256
EXAMPLES::
257
258
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
259
sage: float(A(-3/2))
260
-1.5
261
sage: float(A(-3))
262
-3.0
263
sage: float(-3 + i)
264
Traceback (most recent call last):
265
...
266
TypeError
267
"""
268
if self.is_constant():
269
return float(self[0])
270
raise TypeError
271
272
def _integer_(self, ZZ=None):
273
"""
274
Try to coerce this quaternion to an Integer.
275
276
EXAMPLES::
277
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
278
sage: Integer(A(-3)) # indirect doctest
279
-3
280
sage: type(Integer(A(-3)))
281
<type 'sage.rings.integer.Integer'>
282
sage: Integer(A(-3/2))
283
Traceback (most recent call last):
284
...
285
TypeError: no conversion of this rational to integer
286
sage: Integer(-3 + i)
287
Traceback (most recent call last):
288
...
289
TypeError
290
"""
291
if self.is_constant():
292
return Integer(self[0])
293
raise TypeError
294
295
def _rational_(self):
296
"""
297
Try to coerce this quaternion to a Rational.
298
299
EXAMPLES::
300
301
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(QQ['x']),-5,-2)
302
sage: Rational(Q(2/3)) # indirect doctest
303
2/3
304
sage: Rational(2/3 + i)
305
Traceback (most recent call last):
306
...
307
TypeError
308
"""
309
if self.is_constant():
310
return Rational(self[0])
311
raise TypeError
312
313
def __nonzero__(self):
314
"""
315
Return True if this quaternion is nonzero.
316
317
EXAMPLES::
318
319
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(QQ['x']),-5,-2)
320
sage: bool(i+j)
321
True
322
sage: (i+j).__nonzero__()
323
True
324
sage: Q(0).__nonzero__()
325
False
326
"""
327
return self[0] or self[1] or self[2] or self[3]
328
329
cdef _do_print(self, x,y,z,w):
330
"""
331
Used internally by the print function.
332
333
EXAMPLES::
334
335
sage: Q.<i,j,k> = QuaternionAlgebra(-17,-19)
336
sage: str(i+j+k-3/4) # indirect doctest
337
'-3/4 + i + j + k'
338
"""
339
cdef bint atomic = self._parent._base.is_atomic_repr()
340
v = []
341
i,j,k = self._parent.variable_names()
342
if x:
343
v.append(str(x))
344
c = print_coeff(y,i,atomic)
345
if c: v.append(c)
346
c = print_coeff(z,j,atomic)
347
if c: v.append(c)
348
c = print_coeff(w,k,atomic)
349
if c: v.append(c)
350
if len(v) == 0: return '0'
351
return ' + '.join(v).replace('+ -','- ')
352
353
def _repr_(self):
354
"""
355
Return string representation of this quaternion:
356
357
EXAMPLES::
358
359
sage: R.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(R,-5*x,-2)
360
sage: a = x + i*x^3 + j*x^2 + k*x
361
sage: a._repr_()
362
'x + x^3*i + x^2*j + x*k'
363
sage: a = x+2/3 + i*x^3 + j*(x^2-5/2) + k*x
364
sage: a._repr_()
365
'x + 2/3 + x^3*i + (x^2 - 5/2)*j + x*k'
366
sage: type(a)
367
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
368
sage: Q(0)._repr_()
369
'0'
370
"""
371
return self._do_print(self[0], self[1], self[2], self[3])
372
373
cdef int _cmp_c_impl(self, sage.structure.element.Element right) except -2:
374
"""
375
Comparing elements.
376
377
TESTS::
378
379
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(QQ['x']),-5,-2)
380
sage: a = 1/2 + 2/3*i - 3/4*j + 5/7*k
381
sage: a == a
382
True
383
sage: a == a + 1
384
False
385
sage: a < a + 1
386
True
387
388
sage: K.<x> = QQ[]; Q.<i,j,k> = QuaternionAlgebra(x, 2*x); a = x + 2*x*i + 3*j + (x-2)*k
389
sage: (1/a)*a == 1
390
True
391
sage: (1/a)*a
392
1
393
"""
394
cdef int i
395
for i in range(4):
396
c = cmp(self[i], right[i])
397
if c: return c
398
return 0
399
400
cpdef conjugate(self):
401
"""
402
Return the conjugate of the quaternion: if `\\theta = x + yi + zj + wk`,
403
return `x - yi - zj - wk`; that is, return theta.reduced_trace() - theta.
404
405
EXAMPLES::
406
407
sage: A.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
408
sage: a = 3*i - j + 2
409
sage: type(a)
410
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
411
sage: a.conjugate()
412
2 - 3*i + j
413
414
The "universal" test::
415
416
sage: K.<x,y,z,w,a,b> = QQ[]
417
sage: Q.<i,j,k> = QuaternionAlgebra(a,b)
418
sage: theta = x+y*i+z*j+w*k
419
sage: theta.conjugate()
420
x + (-y)*i + (-z)*j + (-w)*k
421
"""
422
return self.__class__(self._parent, (self[0], -self[1], -self[2], -self[3]))
423
424
cpdef reduced_trace(self):
425
"""
426
Return the reduced trace of self: if `\\theta = x + yi + zj +
427
wk`, then `\\theta` has reduced trace `2x`.
428
429
EXAMPLES::
430
431
sage: K.<x,y,z,w,a,b> = QQ[]
432
sage: Q.<i,j,k> = QuaternionAlgebra(a,b)
433
sage: theta = x+y*i+z*j+w*k
434
sage: theta.reduced_trace()
435
2*x
436
"""
437
return 2*self[0]
438
439
cpdef reduced_norm(self):
440
"""
441
Return the reduced norm of self: if `\\theta = x + yi + zj +
442
wk`, then `\\theta` has reduced norm `x^2 - ay^2 - bz^2 +
443
abw^2`.
444
445
EXAMPLES::
446
447
sage: K.<x,y,z,w,a,b> = QQ[]
448
sage: Q.<i,j,k> = QuaternionAlgebra(a,b)
449
sage: theta = x+y*i+z*j+w*k
450
sage: theta.reduced_norm()
451
w^2*a*b - y^2*a - z^2*b + x^2
452
"""
453
a,b,x,y,z,w = self._parent._a,self._parent._b,self[0],self[1],self[2],self[3]
454
return w*w*a*b - y*y*a - z*z*b + x*x
455
456
def __invert__(self):
457
"""
458
Return inverse of self.
459
460
EXAMPLES::
461
462
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-7,-13)
463
sage: theta = 1/3 - 2/3*i + 4/19*j - 17/3*k
464
sage: (1/theta) * theta
465
1
466
sage: type(theta)
467
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
468
sage: 1/Q(0)
469
Traceback (most recent call last):
470
...
471
ZeroDivisionError: rational division by zero
472
473
Note that the quaternion algebra need not be a division
474
algebra, in which case we can get a ZeroDivisionException::
475
476
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,4,9)
477
sage: theta = 2-i
478
sage: theta.reduced_norm()
479
0
480
sage: 1/theta
481
Traceback (most recent call last):
482
...
483
ZeroDivisionError: rational division by zero
484
485
The ``universal`` test:
486
487
sage: K.<x,y,z,w,a,b> = QQ[]
488
sage: Q.<i,j,k> = QuaternionAlgebra(a,b)
489
sage: theta = x+y*i+z*j+w*k
490
sage: 1/theta == theta.conjugate()/theta.reduced_norm()
491
True
492
"""
493
return self.reduced_norm().__invert__() * self.conjugate()
494
495
cpdef RingElement _div_(self, RingElement right):
496
"""
497
Return quotient of self by right.
498
499
EXAMPLES::
500
501
sage: K.<x> = QQ[]; Q.<i,j,k> = QuaternionAlgebra(x, 2*x);
502
sage: theta = x + 2*x*i + 3*j + (x-2)*k
503
sage: type(theta)
504
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
505
sage: theta._div_(theta)
506
1
507
sage: theta._div_(theta) == 1
508
True
509
"""
510
return self * right.__invert__()
511
512
def reduced_characteristic_polynomial(self, var='x'):
513
"""
514
Return the reduced characteristic polynomial of this
515
quaternion algebra element, which is `X^2 - tX + n`, where `t`
516
is the reduced trace and `n` is the reduced norm.
517
518
INPUT:
519
- var -- string (default: 'x'); indeterminate of characteristic polynomial
520
521
EXAMPLES::
522
523
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
524
sage: i.reduced_characteristic_polynomial()
525
x^2 + 1
526
sage: j.reduced_characteristic_polynomial()
527
x^2 + 2
528
sage: (i+j).reduced_characteristic_polynomial()
529
x^2 + 3
530
sage: (2+j+k).reduced_trace()
531
4
532
sage: (2+j+k).reduced_characteristic_polynomial('T')
533
T^2 - 4*T + 8
534
"""
535
536
R = PolynomialRing(self.base_ring(), var)
537
return R([self.reduced_norm(), -self.reduced_trace(), 1])
538
539
def matrix(self, action='right'):
540
"""
541
Return the matrix of right or left multiplication of self on
542
the basis for the ambient quaternion algebra. In particular,
543
if action is 'right' (the default), returns the matrix of the
544
mapping sending x to x*self.
545
546
INPUT:
547
548
- ``action`` -- (default: 'right') 'right' or 'left'.
549
550
OUTPUT:
551
552
- a matrix
553
554
EXAMPLES::
555
556
sage: Q.<i,j,k> = QuaternionAlgebra(-3,-19)
557
sage: a = 2/3 -1/2*i + 3/5*j - 4/3*k
558
sage: a.matrix()
559
[ 2/3 -1/2 3/5 -4/3]
560
[ 3/2 2/3 4 3/5]
561
[-57/5 -76/3 2/3 1/2]
562
[ 76 -57/5 -3/2 2/3]
563
sage: a.matrix() == a.matrix(action='right')
564
True
565
sage: a.matrix(action='left')
566
[ 2/3 -1/2 3/5 -4/3]
567
[ 3/2 2/3 -4 -3/5]
568
[-57/5 76/3 2/3 -1/2]
569
[ 76 57/5 3/2 2/3]
570
sage: (i*a,j*a,k*a)
571
(3/2 + 2/3*i + 4*j + 3/5*k, -57/5 - 76/3*i + 2/3*j + 1/2*k, 76 - 57/5*i - 3/2*j + 2/3*k)
572
sage: a.matrix(action='foo')
573
Traceback (most recent call last):
574
...
575
ValueError: action must be either 'left' or 'right'
576
577
We test over a more generic base field:
578
579
sage: K.<x> = QQ['x']
580
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(K),-5,-2)
581
sage: a = 1/2*x^2 + 2/3*x*i - 3/4*j + 5/7*k
582
sage: type(a)
583
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
584
sage: a.matrix()
585
[1/2*x^2 2/3*x -3/4 5/7]
586
[-10/3*x 1/2*x^2 -25/7 -3/4]
587
[ 3/2 10/7 1/2*x^2 -2/3*x]
588
[ -50/7 3/2 10/3*x 1/2*x^2]
589
"""
590
if action == 'right':
591
v = [(a*self).coefficient_tuple() for a in self._parent.basis()]
592
elif action == 'left':
593
v = [(self*a).coefficient_tuple() for a in self._parent.basis()]
594
else:
595
raise ValueError, "action must be either 'left' or 'right'"
596
return matrix(self.base_ring(), 4, v, check=False)
597
598
def coefficient_tuple(self):
599
"""
600
Return 4-tuple of coefficients of this quaternion.
601
602
EXAMPLES::
603
604
sage: K.<x> = QQ['x']
605
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(K),-5,-2)
606
sage: a = 1/2*x^2 + 2/3*x*i - 3/4*j + 5/7*k
607
sage: type(a)
608
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
609
sage: a.coefficient_tuple()
610
(1/2*x^2, 2/3*x, -3/4, 5/7)
611
"""
612
return (self[0], self[1], self[2], self[3])
613
614
def pair(self, right):
615
"""
616
Return the result of pairing self and right, which should both
617
be elements of a quaternion algebra. The pairing is
618
(x,y) = (x.conjugate()*y).reduced_trace().
619
620
INPUT:
621
622
- ``right`` -- quaternion
623
624
EXAMPLES::
625
626
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
627
sage: (1+i+j-2*k).pair(2/3+5*i-3*j+k)
628
-26/3
629
sage: x = 1+i+j-2*k; y = 2/3+5*i-3*j+k
630
sage: x.pair(y)
631
-26/3
632
sage: y.pair(x)
633
-26/3
634
sage: (x.conjugate()*y).reduced_trace()
635
-26/3
636
"""
637
return (self.conjugate() * right).reduced_trace()
638
639
cdef class QuaternionAlgebraElement_generic(QuaternionAlgebraElement_abstract):
640
"""
641
TESTS:
642
643
We test pickling::
644
645
sage: R.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(R,-5*x,-2)
646
sage: theta = x + i*x^3 + j*x^2 + k*x
647
sage: theta == loads(dumps(theta))
648
True
649
"""
650
def __init__(self, parent, v):
651
"""
652
Create a quaternion over some general base field.
653
654
EXAMPLES::
655
656
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
657
sage: sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic(Q, (x,1,-7,2/3*x^3))
658
x + i + (-7)*j + 2/3*x^3*k
659
"""
660
self._parent = parent
661
self.x, self.y, self.z, self.w = to_quaternion(parent._base, v)
662
663
def __getitem__(self, int i):
664
"""
665
EXAMPLES::
666
667
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(QQ['x']),-5,-2)
668
sage: theta = 1/2 + 2/3*i - 3/4*j + 5/7*k
669
sage: type(theta)
670
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
671
sage: list(theta)
672
[1/2, 2/3, -3/4, 5/7]
673
"""
674
if i == 0:
675
return self.x
676
elif i == 1:
677
return self.y
678
elif i == 2:
679
return self.z
680
elif i == 3:
681
return self.w
682
else:
683
raise IndexError, "quaternion element index out of range"
684
685
def __reduce__(self):
686
"""
687
Used for pickling.
688
689
TESTS::
690
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
691
sage: theta = 1/x + x*i - (x+1)*j + 2/(3*x^3+5)*k
692
sage: loads(dumps(theta)) == theta
693
True
694
sage: type(theta)
695
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
696
"""
697
return (unpickle_QuaternionAlgebraElement_generic_v0,
698
(self._parent, (self.x, self.y, self.z, self.w)))
699
700
cpdef ModuleElement _add_(self, ModuleElement _right):
701
"""
702
Return the sum of self and _right.
703
704
EXAMPLES::
705
706
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
707
sage: (x+i+j+x^3*k) + (x-i-j+ (2/3*x^3+x)*k) # indirect doctest
708
2*x + (5/3*x^3 + x)*k
709
sage: type(i)
710
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
711
"""
712
cdef QuaternionAlgebraElement_generic right = _right
713
# TODO -- make this, etc. use PY_NEW
714
return QuaternionAlgebraElement_generic(self._parent, (self.x + right.x, self.y + right.y, self.z + right.z, self.w + right.w))
715
716
cpdef ModuleElement _sub_(self, ModuleElement _right):
717
"""
718
Return the difference of self and _right.
719
720
EXAMPLES::
721
722
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
723
sage: type(i)
724
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
725
sage: (x+i+j+x^3*k)._sub_(x-i-j+ (2/3*x^3+x)*k)
726
2*i + 2*j + (1/3*x^3 - x)*k
727
"""
728
cdef QuaternionAlgebraElement_generic right = _right
729
return QuaternionAlgebraElement_generic(self._parent, (self.x - right.x, self.y - right.y, self.z - right.z, self.w - right.w))
730
731
cpdef RingElement _mul_(self, RingElement _right):
732
"""
733
Return the product of self and _right.
734
735
EXAMPLES::
736
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
737
sage: type(i)
738
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
739
sage: (x+i+j+x^3*k)._mul_(x-i-j+ (2/3*x^3+x)*k)
740
-20/3*x^6 - 10*x^4 + x^2 + 7 + (10/3*x^3 + 2*x)*i + (-25/3*x^3 - 5*x)*j + (5/3*x^4 + x^2)*k
741
"""
742
cdef QuaternionAlgebraElement_generic right = _right
743
744
a = self._parent._a
745
b = self._parent._b
746
747
x1, y1, z1, w1 = self.x, self.y, self.z, self.w
748
x2, y2, z2, w2 = right.x, right.y, right.z, right.w
749
750
#x = x1*x2 + y1*y2*a + z1*z2*b - w1*w2*a*b
751
#y = x1*y2 + y1*x2 - z1*w2*b + w1*z2*b
752
#z = x1*z2 + y1*w2 + z1*x2 - w1*y2*a
753
#w = x1*w2 + y1*z2 - z1*y2 + w1*x2
754
755
t1 = x1 * x2
756
t2 = y1 * y2
757
t3 = z1 * z2
758
t4 = w1 * w2
759
t5 = x2 * z1
760
t6 = y2 * w1
761
t7 = x1 * z2
762
t8 = y1 * w2
763
764
x = t1 + a * t2 + b * (t3 - a*t4)
765
y = (x1 + y1)*(x2 + y2) - t1 - t2 + b*( (z1 + w1)*(z2 - w2) - t3 + t4)
766
z = t5 - a*t6 + t7 + a*t8
767
w = (x2 - y2)*(z1 + w1) - t5 + t6 + (x1 + y1)*(z2 + w2) - t7 - t8
768
769
return QuaternionAlgebraElement_generic(self._parent, (x, y, z, w))
770
771
def _repr_(self):
772
"""
773
Print representation of self.
774
775
EXAMPLES::
776
777
sage: K.<x> = Frac(QQ['x']); Q.<i,j,k> = QuaternionAlgebra(K,-5,-2)
778
sage: theta = 1/x + x*i - (x+1)*j + 2/(3*x^3+5)*k
779
sage: type(theta)
780
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
781
sage: theta._repr_()
782
'1/x + x*i + (-x - 1)*j + (2/(3*x^3 + 5))*k'
783
"""
784
return self._do_print(self.x, self.y, self.z, self.w)
785
786
787
cdef class QuaternionAlgebraElement_rational_field(QuaternionAlgebraElement_abstract):
788
"""
789
TESTS:
790
791
We test pickling::
792
793
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
794
sage: i + j + k == loads(dumps(i+j+k))
795
True
796
"""
797
798
# Implementation Notes:
799
#
800
# A Quaternion algebra element (call it a) over Q are implemented as a 4-tuple of
801
# integers x, y, z, w and a denominator d, all of type mpz_t, such that
802
#
803
# theta = (1/d) * (x + y * i + z * j + w * k)
804
#
805
# (although different letters may be specified instead of i, j, and k, if desired).
806
#
807
# Inside the element we also store mpz_t integers a and b, where
808
#
809
# i^2 = a and j^2 = b
810
811
def __cinit__(self):
812
"""
813
Initialize C variables.
814
815
EXAMPLES::
816
817
sage: QuaternionAlgebra(QQ,-5,-2)([1/2,-1/3,2/3,4/5]) # implicit doctest
818
1/2 - 1/3*i + 2/3*j + 4/5*k
819
"""
820
mpz_init(self.x)
821
mpz_init(self.y)
822
mpz_init(self.z)
823
mpz_init(self.w)
824
mpz_init(self.a)
825
mpz_init(self.b)
826
mpz_init(self.d)
827
828
def __dealloc__(self):
829
mpz_clear(self.x)
830
mpz_clear(self.y)
831
mpz_clear(self.z)
832
mpz_clear(self.w)
833
mpz_clear(self.a)
834
mpz_clear(self.b)
835
mpz_clear(self.d)
836
837
def _rational_(self):
838
"""
839
Try to coerce this quaternion to a Rational.
840
841
EXAMPLES::
842
843
sage: A.<i,j,k> = QuaternionAlgebra(-1,-2)
844
sage: Rational(A(-2/3)) # indirect doctest
845
-2/3
846
sage: Rational(i)
847
Traceback (most recent call last):
848
...
849
TypeError
850
"""
851
cdef Rational x = Rational()
852
if self.is_constant():
853
mpq_set_num(x.value, self.x)
854
mpq_set_den(x.value, self.d)
855
mpq_canonicalize(x.value)
856
return x
857
raise TypeError
858
859
cpdef bint is_constant(self):
860
"""
861
Return True if this quaternion is constant, i.e., has no i, j, or k term.
862
863
OUTPUT:
864
bool
865
866
EXAMPLES::
867
868
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
869
sage: A(1/3).is_constant()
870
True
871
sage: A(-1).is_constant()
872
True
873
sage: (1+i).is_constant()
874
False
875
sage: j.is_constant()
876
False
877
"""
878
return not (mpz_sgn(self.y) or mpz_sgn(self.z) or mpz_sgn(self.w))
879
880
def __nonzero__(self):
881
"""
882
Return True if this quaternion is nonzero.
883
884
EXAMPLES::
885
886
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
887
sage: bool(1+j+k)
888
True
889
sage: A(0).__nonzero__()
890
False
891
"""
892
return bool(mpz_sgn(self.x) or mpz_sgn(self.y) or mpz_sgn(self.z) or mpz_sgn(self.w))
893
894
cdef int _cmp_c_impl(self, sage.structure.element.Element _right) except -2:
895
"""
896
Compare two quaternions. The comparison is fairly arbitrary
897
-- first the denominators are compared and if equal then each
898
of the other coefficients are compared.
899
900
TESTS::
901
902
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
903
sage: i < j
904
False
905
sage: -i < j
906
True
907
sage: i == i
908
True
909
"""
910
cdef QuaternionAlgebraElement_rational_field right = _right
911
cdef int i
912
i = mpz_cmp(self.d, right.d)
913
if i < 0: return -1
914
elif i > 0: return 1
915
i = mpz_cmp(self.x, right.x)
916
if i < 0: return -1
917
elif i > 0: return 1
918
i = mpz_cmp(self.y, right.y)
919
if i < 0: return -1
920
elif i > 0: return 1
921
i = mpz_cmp(self.z, right.z)
922
if i < 0: return -1
923
elif i > 0: return 1
924
i = mpz_cmp(self.w, right.w)
925
if i < 0: return -1
926
elif i > 0: return 1
927
return 0
928
929
def __init__(self, parent, v):
930
"""
931
Setup element data from parent and coordinates.
932
933
EXAMPLES::
934
935
sage: A.<i,j,k>=QuaternionAlgebra(-4,-5)
936
sage: A(2/3)
937
2/3
938
sage: type(A(2/3))
939
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
940
941
sage: A([-1/2,-10/3,-2/3,-4/5]) # implicit doctest
942
-1/2 - 10/3*i - 2/3*j - 4/5*k
943
sage: A(vector([1,2/3,3/4,4/5]))
944
1 + 2/3*i + 3/4*j + 4/5*k
945
946
::
947
948
sage: QA = QuaternionAlgebra(QQ, -1, -1)
949
sage: foo = QA(3.0); foo
950
3
951
sage: parent(foo)
952
Quaternion Algebra (-1, -1) with base ring Rational Field
953
sage: foo[0]
954
3
955
sage: parent(foo[0])
956
Rational Field
957
"""
958
self._parent = parent
959
960
# cache a and b
961
mpz_set(self.a, (<Integer>parent._a).value)
962
mpz_set(self.b, (<Integer>parent._b).value)
963
964
cdef Rational x, y, z, w
965
cdef mpz_t lcm
966
cdef mpq_t lcm_rat
967
if not isinstance(v, (list,tuple)):
968
try:
969
x = Rational(v)
970
mpz_set(self.x, mpq_numref(x.value))
971
mpz_set_si(self.y, 0)
972
mpz_set_si(self.z, 0)
973
mpz_set_si(self.w, 0)
974
mpz_set(self.d, mpq_denref(x.value))
975
return
976
except TypeError:
977
pass
978
v = list(v)
979
# Now v is definitely a list or tuple, and we convert each
980
# entry to a rational, then clear denominators, etc.
981
x = Rational(v[0])
982
y = Rational(v[1])
983
z = Rational(v[2])
984
w = Rational(v[3])
985
mpz_init(lcm)
986
mpz_lcm(lcm, mpq_denref(x.value), mpq_denref(y.value))
987
mpz_lcm(lcm, lcm, mpq_denref(z.value))
988
mpz_lcm(lcm, lcm, mpq_denref(w.value))
989
if mpz_cmp_si(lcm, 1):
990
mpz_init_set(mpq_numref(lcm_rat), lcm)
991
mpz_init_set_si(mpq_denref(lcm_rat), 1)
992
mpq_mul(x.value, x.value, lcm_rat)
993
mpq_mul(y.value, y.value, lcm_rat)
994
mpq_mul(z.value, z.value, lcm_rat)
995
mpq_mul(w.value, w.value, lcm_rat)
996
mpq_clear(lcm_rat)
997
mpz_set(self.x, mpq_numref(x.value))
998
mpz_set(self.y, mpq_numref(y.value))
999
mpz_set(self.z, mpq_numref(z.value))
1000
mpz_set(self.w, mpq_numref(w.value))
1001
mpz_set(self.d, lcm)
1002
mpz_clear(lcm)
1003
1004
1005
def __getitem__(self, int i):
1006
"""
1007
TESTS::
1008
1009
sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
1010
sage: theta = 1/2 + 2/3*i - 3/4*j + 5/7*k
1011
sage: type(theta)
1012
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
1013
sage: list(theta)
1014
[1/2, 2/3, -3/4, 5/7]
1015
"""
1016
cdef Rational r = Rational()
1017
if i == 0:
1018
mpq_set_num(r.value, self.x)
1019
elif i == 1:
1020
mpq_set_num(r.value, self.y)
1021
elif i == 2:
1022
mpq_set_num(r.value, self.z)
1023
elif i == 3:
1024
mpq_set_num(r.value, self.w)
1025
else:
1026
raise IndexError, "quaternion element index out of range"
1027
mpq_set_den(r.value, self.d)
1028
mpq_canonicalize(r.value)
1029
return r
1030
1031
def __reduce__(self):
1032
"""
1033
Used for pickling.
1034
1035
TESTS::
1036
1037
sage: K.<x> = QQ[]
1038
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(K),-5,-19)
1039
sage: theta = 1/2 + 2/3*i - 3/4*j + 5/7*k
1040
sage: type(theta)
1041
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
1042
sage: loads(dumps(theta)) == theta
1043
True
1044
1045
"""
1046
return (unpickle_QuaternionAlgebraElement_rational_field_v0,
1047
(self._parent, (self[0], self[1], self[2], self[3])))
1048
1049
cpdef ModuleElement _add_(self, ModuleElement _right):
1050
"""
1051
EXAMPLES::
1052
1053
sage: Q.<i,j,k> = QuaternionAlgebra(15)
1054
sage: type(i)
1055
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
1056
sage: (2/3 + 3/4*i + 5/6*j + 7/8*k)._add_(-2/3 - 3/4*i + 5/6*j + 7/8*k)
1057
5/3*j + 7/4*k
1058
"""
1059
1060
# Given two quaternion algebra elements
1061
# theta = (1/d1)*(x1 + y1 * i + z1 * j + w1 * k)
1062
# nu = (1/d2)*(x2 + y2 * i + z2 * j + w2 * k)
1063
#
1064
# we compute their sum as
1065
#
1066
# (theta + nu) = (1/d3)*(x3 + y3 * i + z3 * j + w3 * k)
1067
#
1068
# with d3 = d1 * d2
1069
# x3 = d1 * x2 + d2 * x1
1070
# y3 = d1 * y2 + d2 * y1
1071
# z3 = d1 * z2 + d2 * z1
1072
# w3 = d1 * w2 + d2 * w1
1073
#
1074
# and then we reduce the sum by dividing everything
1075
# by the gcd of d3, x3, y3, z3, and w3
1076
1077
cdef QuaternionAlgebraElement_rational_field right = _right
1078
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1079
result._parent = self._parent
1080
1081
mpz_mul(U1, self.x, right.d) # U1 = x1 * d2
1082
mpz_mul(U2, right.x, self.d) # U2 = x2 * d1
1083
mpz_add(result.x, U1, U2) # x3 = x1 * d2 + x2 * d1
1084
1085
mpz_mul(U1, self.y, right.d) # U1 = y1 * d2
1086
mpz_mul(U2, right.y, self.d) # U2 = y2 * d1
1087
mpz_add(result.y, U1, U2) # x3 = y1 * d2 + y2 * d1
1088
1089
mpz_mul(U1, self.z, right.d) # U1 = z1 * d2
1090
mpz_mul(U2, right.z, self.d) # U2 = z2 * d1
1091
mpz_add(result.z, U1, U2) # z3 = z1 * d2 + z2 * d1
1092
1093
mpz_mul(U1, self.w, right.d) # U1 = w1 * d2
1094
mpz_mul(U2, right.w, self.d) # U2 = w2 * d1
1095
mpz_add(result.w, U1, U2) # w3 = w1 * d2 + w2 * d1
1096
1097
mpz_mul(result.d, self.d, right.d) # d3 = d1 * d2
1098
1099
result.canonicalize()
1100
1101
mpz_set(result.a, self.a)
1102
mpz_set(result.b, self.b)
1103
return result
1104
1105
cpdef ModuleElement _sub_(self, ModuleElement _right):
1106
"""
1107
EXAMPLES::
1108
1109
sage: Q.<i,j,k> = QuaternionAlgebra(15)
1110
sage: type(i)
1111
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
1112
sage: (2/3 + 3/4*i + 5/6*j + 7/8*k)._sub_(-2/3 - 3/4*i + 5/6*j + 7/8*k)
1113
4/3 + 3/2*i
1114
"""
1115
cdef QuaternionAlgebraElement_rational_field right = _right
1116
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1117
result._parent = self._parent
1118
1119
# Implementation Note: To obtain _sub_, we simply replace every occurrence of
1120
# "add" in _add_ with "sub"; that is, we s/add/sub to get _sub_
1121
1122
mpz_mul(U1, self.x, right.d)
1123
mpz_mul(U2, right.x, self.d)
1124
mpz_sub(result.x, U1, U2)
1125
1126
mpz_mul(U1, self.y, right.d)
1127
mpz_mul(U2, right.y, self.d)
1128
mpz_sub(result.y, U1, U2)
1129
1130
mpz_mul(U1, self.z, right.d)
1131
mpz_mul(U2, right.z, self.d)
1132
mpz_sub(result.z, U1, U2)
1133
1134
mpz_mul(U1, self.w, right.d)
1135
mpz_mul(U2, right.w, self.d)
1136
mpz_sub(result.w, U1, U2)
1137
1138
mpz_mul(result.d, self.d, right.d)
1139
1140
result.canonicalize()
1141
1142
mpz_set(result.a, self.a)
1143
mpz_set(result.b, self.b)
1144
return result
1145
1146
cpdef RingElement _mul_(self, RingElement _right):
1147
"""
1148
EXAMPLES::
1149
1150
sage: Q.<i,j,k> = QuaternionAlgebra(15)
1151
sage: type(i)
1152
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
1153
sage: (2/3 + 3/4*i + 5/6*j + 7/8*k)._mul_(-2/3 - 3/4*i + 5/6*j + 7/8*k)
1154
9331/576 - i - 63/16*j + 5/4*k
1155
"""
1156
1157
# We use the following formula for multiplication:
1158
#
1159
# Given two quaternion algebra elements
1160
#
1161
# theta = (1/d1)*(x1 + y1 * i + z1 * j + w1 * k)
1162
# nu = (1/d2)*(x2 + y2 * i + z2 * j + w2 * k)
1163
#
1164
# we compute their product as
1165
#
1166
# theta*nu = (1/d3)*(x3 + y3 * i + z3 * j + w3 * k)
1167
#
1168
# where
1169
# x3 = t1 + a * t2 + b * (t3 - a*t4)
1170
# y3 = s1*(x2 + y2) - t1 - t2 + b*( s2*(z2 - w2) - t3 + t4)
1171
# z3 = t5 - a*t6 + t7 + a*t8
1172
# w3 = (x2 - y2)*s2 - t5 + t6 + s1*(z2 + w2) - t7 - t8
1173
#
1174
# and where
1175
# t1 = x1 * x2
1176
# t2 = y1 * y2
1177
# t3 = z1 * z2
1178
# t4 = w1 * w2
1179
# t5 = x2 * z1
1180
# t6 = y2 * w1
1181
# t7 = x1 * z2
1182
# t8 = y1 * w2
1183
#
1184
# s1 = x1 + y1
1185
# s2 = z1 + w1
1186
#
1187
# This takes more integer addition operations but fewer integer multiplication
1188
# operations than the "straightforward" multiplication method.
1189
#
1190
# There might be a way to optimize this formula further.
1191
1192
1193
cdef QuaternionAlgebraElement_rational_field right = _right
1194
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1195
result._parent = self._parent
1196
1197
mpz_set(result.a, self.a)
1198
mpz_set(result.b, self.b)
1199
1200
mpz_mul(T1, self.x, right.x) # t1 = x1 * x2
1201
mpz_mul(T2, self.y, right.y) # t2 = y1 * y2
1202
mpz_mul(t3, self.z, right.z) # t3 = z1 * z2
1203
mpz_mul(t4, self.w, right.w) # t4 = w1 * w2
1204
mpz_mul(t5, right.x, self.z) # t5 = x2 * z1
1205
mpz_mul(t6, right.y, self.w) # t6 = y2 * w1
1206
mpz_mul(t7, self.x, right.z) # t7 = x1 * z2
1207
mpz_mul(t8, self.y, right.w) # t8 = y1 * w2
1208
1209
mpz_add(s1, self.x, self.y) # s1 = x1 + y1
1210
mpz_add(s2, self.z, self.w) # s2 = z1 + w1
1211
1212
#------------------
1213
1214
mpz_mul(U1, self.a, t4)
1215
mpz_sub(U1, t3, U1)
1216
mpz_mul(U1, U1, self.b)
1217
mpz_mul(U2, self.a, T2)
1218
mpz_add(result.x, T1, U2)
1219
mpz_add(result.x, result.x, U1)
1220
1221
#------------------
1222
1223
mpz_sub(U1, right.z, right.w)
1224
mpz_mul(U1, U1, s2)
1225
mpz_sub(U1, U1, t3)
1226
mpz_add(U1, U1, t4)
1227
mpz_mul(U1, U1, self.b)
1228
mpz_sub(U1, U1, T2)
1229
mpz_sub(U1, U1, T1)
1230
mpz_add(U2, right.x, right.y)
1231
mpz_mul(U2, s1, U2)
1232
mpz_add(result.y, U1, U2)
1233
1234
#------------------
1235
1236
mpz_mul(U1, self.a, t8)
1237
mpz_add(U1, U1, t7)
1238
mpz_mul(U2, self.a, t6)
1239
mpz_sub(U1, U1, U2)
1240
mpz_add(result.z, U1, t5)
1241
1242
#------------------
1243
1244
mpz_add(U1, right.z, right.w)
1245
mpz_mul(U1, U1, s1)
1246
mpz_sub(U1, U1, t7)
1247
mpz_sub(U1, U1, t8)
1248
mpz_add(U1, U1, t6)
1249
mpz_sub(U1, U1, t5)
1250
mpz_sub(U2, right.x, right.y)
1251
mpz_mul(U2, U2, s2)
1252
mpz_add(result.w, U1, U2)
1253
1254
1255
1256
mpz_mul(result.d, self.d, right.d)
1257
1258
result.canonicalize()
1259
1260
return result
1261
1262
cpdef reduced_norm(self):
1263
"""
1264
Return the reduced norm of self. Given a quaternion
1265
`x+yi+zj+wk`, this is `x^2 - ay^2 - bz^2 + abw^2`.
1266
1267
EXAMPLES::
1268
1269
sage: K.<i,j,k> = QuaternionAlgebra(QQ, -5, -2)
1270
sage: i.reduced_norm()
1271
5
1272
sage: j.reduced_norm()
1273
2
1274
sage: a = 1/3 + 1/5*i + 1/7*j + k
1275
sage: a.reduced_norm()
1276
22826/2205
1277
"""
1278
mpz_mul(U1, self.x, self.x) # U1 = x*x
1279
mpz_mul(U2, self.b, self.z) # U2 = b*x
1280
mpz_mul(U2, U2, self.z) # U2 = b*z*z
1281
mpz_sub(U2, U1, U2) # U2 = -b*z*z + x*x
1282
1283
mpz_mul(U1, self.y, self.a) # U1 = a*y
1284
mpz_mul(U1, U1, self.y) # U1 = a*y*y
1285
mpz_sub(U2, U2, U1) # U2 = -a*y*y - b*z*z + x*x
1286
1287
mpz_mul(U1, self.w, self.w) # U1 = w*w
1288
mpz_mul(U1, U1, self.a) # U1 = w*w*a
1289
mpz_mul(U1, U1, self.b) # U1 = w*w*a*b
1290
1291
mpz_add(U1, U1, U2) # U1 = w*w*a*b - a*y*y - b*z*z + x*x
1292
1293
mpz_mul(U2, self.d, self.d)
1294
1295
cdef Rational result = PY_NEW(Rational)
1296
mpq_set_num(result.value, U1)
1297
mpq_set_den(result.value, U2)
1298
mpq_canonicalize(result.value)
1299
1300
return result
1301
1302
cpdef conjugate(self):
1303
"""
1304
Return the conjugate of this quaternion.
1305
1306
EXAMPLES::
1307
1308
sage: A.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
1309
sage: a = 3*i - j + 2
1310
sage: type(a)
1311
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
1312
sage: a.conjugate()
1313
2 - 3*i + j
1314
sage: b = 1 + 1/3*i + 1/5*j - 1/7*k
1315
sage: b.conjugate()
1316
1 - 1/3*i - 1/5*j + 1/7*k
1317
"""
1318
1319
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1320
result._parent = self._parent
1321
1322
mpz_set(result.a, self.a)
1323
mpz_set(result.b, self.b)
1324
mpz_set(result.d, self.d)
1325
1326
mpz_set(result.x, self.x)
1327
mpz_mul_si(result.y, self.y, -1)
1328
mpz_mul_si(result.z, self.z, -1)
1329
mpz_mul_si(result.w, self.w, -1)
1330
1331
return result
1332
1333
cpdef reduced_trace(self):
1334
"""
1335
Return the reduced trace of self, which is `2x` if self is `x+iy+zj+wk`.
1336
1337
EXAMPLES::
1338
1339
sage: K.<i,j,k> = QuaternionAlgebra(QQ, -5, -2)
1340
sage: i.reduced_trace()
1341
0
1342
sage: j.reduced_trace()
1343
0
1344
sage: a = 1/3 + 1/5*i + 1/7*j + k
1345
sage: a.reduced_trace()
1346
2/3
1347
"""
1348
#return 2*self[0]
1349
1350
mpz_mul_si(U1, self.x, 2)
1351
cdef Rational result = PY_NEW(Rational)
1352
mpq_set_num(result.value, U1)
1353
mpq_set_den(result.value, self.d)
1354
mpq_canonicalize(result.value)
1355
return result
1356
1357
cdef inline canonicalize(self):
1358
"""
1359
Put the representation of this quaternion element into
1360
smallest form. For `a = (1/d)(x + yi + zj + wk)` we
1361
divide `a`, `x`, `y`, `z`, and `w` by the gcd of all of them.
1362
1363
TESTS::
1364
sage: K.<i,j,k> = QuaternionAlgebra(QQ, -10, -7)
1365
sage: (1/4 + 1/2 * i + 1/7 * j + 1/28 * k)*14*i # implicit doctest
1366
-70 + 7/2*i + 5*j - 2*k
1367
"""
1368
1369
# Note: this function changes the module-level global variable
1370
# U1, so it isn't always safe to use this in the middle of
1371
# another function. Normally this function is called
1372
# at the end of an arithmetic routine, so this is fine.
1373
1374
# Implementation-wise, we compute the GCD's one at a time,
1375
# and quit if it ever becomes one
1376
1377
1378
mpz_gcd(U1, self.d, self.x)
1379
if mpz_cmp_ui(U1, 1) != 0:
1380
mpz_gcd(U1, U1, self.y)
1381
if mpz_cmp_ui(U1, 1) != 0:
1382
mpz_gcd(U1, U1, self.z)
1383
if mpz_cmp_ui(U1, 1) != 0:
1384
mpz_gcd(U1, U1, self.w)
1385
if mpz_cmp_ui(U1, 1) != 0:
1386
# at this point U1 actually contains the gcd of all the terms, and we divide
1387
mpz_divexact(self.d, self.d, U1)
1388
mpz_divexact(self.x, self.x, U1)
1389
mpz_divexact(self.y, self.y, U1)
1390
mpz_divexact(self.z, self.z, U1)
1391
mpz_divexact(self.w, self.w, U1)
1392
1393
def denominator(self):
1394
"""
1395
Return the lowest common multiple of the denominators of the coefficients
1396
of i, j and k for this quaternion.
1397
1398
EXAMPLES::
1399
1400
sage: A = QuaternionAlgebra(QQ, -1, -1)
1401
sage: A.<i,j,k> = QuaternionAlgebra(QQ, -1, -1)
1402
sage: a = (1/2) + (1/5)*i + (5/12)*j + (1/13)*k
1403
sage: a
1404
1/2 + 1/5*i + 5/12*j + 1/13*k
1405
sage: a.denominator()
1406
780
1407
sage: lcm([2, 5, 12, 13])
1408
780
1409
sage: (a * a).denominator()
1410
608400
1411
sage: (a + a).denominator()
1412
390
1413
"""
1414
1415
cdef Integer d = Integer()
1416
mpz_set(d.value, self.d)
1417
return d
1418
1419
def denominator_and_integer_coefficient_tuple(self):
1420
"""
1421
Return 5-tuple d, x, y, z, w, where this rational quaternion
1422
is equal to `(x + yi + zj + wk)/d` and x, y, z, w do not share
1423
a common factor with d.
1424
1425
OUTPUT:
1426
5-tuple of Integers
1427
1428
EXAMPLES::
1429
1430
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
1431
sage: (2 + 3*i + 4/3*j - 5*k).denominator_and_integer_coefficient_tuple()
1432
(3, 6, 9, 4, -15)
1433
"""
1434
cdef Integer d = Integer()
1435
cdef Integer y = Integer()
1436
cdef Integer x = Integer()
1437
cdef Integer z = Integer()
1438
cdef Integer w = Integer()
1439
1440
mpz_set(d.value, self.d)
1441
mpz_set(x.value, self.x)
1442
mpz_set(y.value, self.y)
1443
mpz_set(z.value, self.z)
1444
mpz_set(w.value, self.w)
1445
1446
return (d, x, y, z, w)
1447
1448
def integer_coefficient_tuple(self):
1449
"""
1450
Returns integer part of this quaternion, ignoring the common denominator.
1451
1452
OUTPUT:
1453
4-tuple of Integers
1454
1455
EXAMPLES::
1456
1457
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
1458
sage: (2 + 3*i + 4/3*j - 5*k).integer_coefficient_tuple()
1459
(6, 9, 4, -15)
1460
"""
1461
cdef Integer y = Integer()
1462
cdef Integer x = Integer()
1463
cdef Integer z = Integer()
1464
cdef Integer w = Integer()
1465
1466
mpz_set(x.value, self.x)
1467
mpz_set(y.value, self.y)
1468
mpz_set(z.value, self.z)
1469
mpz_set(w.value, self.w)
1470
1471
return (x, y, z, w)
1472
1473
1474
1475
def coefficient_tuple(self):
1476
"""
1477
Return 4-tuple of rational numbers which are the coefficients of this quaternion.
1478
1479
EXAMPLES::
1480
1481
sage: A.<i,j,k>=QuaternionAlgebra(-1,-2)
1482
sage: (2/3 + 3/5*i + 4/3*j - 5/7*k).coefficient_tuple()
1483
(2/3, 3/5, 4/3, -5/7)
1484
"""
1485
cdef Rational x = Rational()
1486
cdef Rational y = Rational()
1487
cdef Rational z = Rational()
1488
cdef Rational w = Rational()
1489
1490
mpq_set_num(x.value, self.x)
1491
mpq_set_num(y.value, self.y)
1492
mpq_set_num(z.value, self.z)
1493
mpq_set_num(w.value, self.w)
1494
1495
mpq_set_den(x.value, self.d)
1496
mpq_set_den(y.value, self.d)
1497
mpq_set_den(z.value, self.d)
1498
mpq_set_den(w.value, self.d)
1499
1500
mpq_canonicalize(x.value)
1501
mpq_canonicalize(y.value)
1502
mpq_canonicalize(z.value)
1503
mpq_canonicalize(w.value)
1504
return (x, y, z, w)
1505
1506
def _multiply_by_integer(self, Integer n):
1507
"""
1508
Return the product of self times the integer n.
1509
1510
EXAMPLES::
1511
sage: A = QuaternionAlgebra(7)
1512
sage: a = A.random_element()
1513
sage: 5*a == a._multiply_by_integer(5)
1514
True
1515
"""
1516
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1517
result._parent = self._parent
1518
1519
mpz_set(result.a, self.a)
1520
mpz_set(result.b, self.b)
1521
1522
if mpz_divisible_p(self.d, n.value) != 0:
1523
mpz_divexact(result.d, self.d, n.value)
1524
mpz_set(result.x, self.x)
1525
mpz_set(result.y, self.y)
1526
mpz_set(result.z, self.z)
1527
mpz_set(result.w, self.w)
1528
return result
1529
if mpz_divisible_p(n.value, self.d):
1530
mpz_divexact(T1, n.value, self.d)
1531
else:
1532
mpz_set(T1, n.value)
1533
1534
mpz_set(result.d, self.d)
1535
mpz_mul(result.x, self.x, T1)
1536
mpz_mul(result.y, self.y, T1)
1537
mpz_mul(result.z, self.z, T1)
1538
mpz_mul(result.w, self.w, T1)
1539
1540
return result
1541
1542
def _divide_by_integer(self, Integer n):
1543
"""
1544
Return the quotient of self by the integer n.
1545
1546
EXAMPLES::
1547
sage: A = QuaternionAlgebra(7)
1548
sage: a = A.random_element()
1549
sage: a/5 == a._divide_by_integer(5)
1550
True
1551
sage: a._divide_by_integer(0)
1552
Traceback (most recent call last):
1553
...
1554
ZeroDivisionError
1555
"""
1556
if mpz_sgn(n.value) == 0:
1557
raise ZeroDivisionError
1558
1559
cdef QuaternionAlgebraElement_rational_field result = <QuaternionAlgebraElement_rational_field> PY_NEW(QuaternionAlgebraElement_rational_field)
1560
result._parent = self._parent
1561
1562
mpz_set(result.a, self.a)
1563
mpz_set(result.b, self.b)
1564
1565
mpz_mul(result.d, self.d, n.value)
1566
mpz_set(result.x, self.x)
1567
mpz_set(result.y, self.y)
1568
mpz_set(result.z, self.z)
1569
mpz_set(result.w, self.w)
1570
1571
result.canonicalize()
1572
1573
return result
1574
1575
1576
1577
cdef class QuaternionAlgebraElement_number_field(QuaternionAlgebraElement_abstract):
1578
def __cinit__(self):
1579
"""
1580
Allocate memory for this quaternion over a number field.
1581
"""
1582
fmpz_poly_init(self.x)
1583
fmpz_poly_init(self.y)
1584
fmpz_poly_init(self.z)
1585
fmpz_poly_init(self.w)
1586
fmpz_poly_init(self.a)
1587
fmpz_poly_init(self.b)
1588
fmpz_poly_init(self.modulus)
1589
mpz_init(self.d)
1590
1591
def __dealloc__(self):
1592
"""
1593
Free memory used by this quaternion over a number field.
1594
"""
1595
fmpz_poly_clear(self.x)
1596
fmpz_poly_clear(self.y)
1597
fmpz_poly_clear(self.z)
1598
fmpz_poly_clear(self.w)
1599
fmpz_poly_clear(self.a)
1600
fmpz_poly_clear(self.b)
1601
fmpz_poly_clear(self.modulus)
1602
mpz_clear(self.d)
1603
1604
def __init__(self, parent, v):
1605
"""
1606
EXAMPLES::
1607
1608
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K,-a,a+1)
1609
sage: Q([a,-2/3,a^2-1/2,a*2]) # implicit doctest
1610
a + (-2/3)*i + (a^2 - 1/2)*j + 2*a*k
1611
"""
1612
self._parent = parent
1613
x, y, z, w = to_quaternion(parent._base, v)
1614
cdef NumberFieldElement a = <NumberFieldElement>(parent._base(parent._a))
1615
cdef NumberFieldElement b = <NumberFieldElement>(parent._base(parent._b))
1616
ZZX_to_fmpz_poly(self.x, (<NumberFieldElement>x).__numerator)
1617
ZZX_to_fmpz_poly(self.y, (<NumberFieldElement>y).__numerator)
1618
ZZX_to_fmpz_poly(self.z, (<NumberFieldElement>z).__numerator)
1619
ZZX_to_fmpz_poly(self.w, (<NumberFieldElement>w).__numerator)
1620
1621
ZZ_to_mpz(&T1, &(<NumberFieldElement>x).__denominator)
1622
ZZ_to_mpz(&T2, &(<NumberFieldElement>y).__denominator)
1623
ZZ_to_mpz(&t3, &(<NumberFieldElement>z).__denominator)
1624
ZZ_to_mpz(&t4, &(<NumberFieldElement>w).__denominator)
1625
1626
mpz_lcm(self.d, T1, T2)
1627
mpz_lcm(self.d, self.d, t3)
1628
mpz_lcm(self.d, self.d, t4)
1629
1630
mpz_divexact(T1, self.d, T1)
1631
mpz_divexact(T2, self.d, T2)
1632
mpz_divexact(t3, self.d, t3)
1633
mpz_divexact(t4, self.d, t4)
1634
1635
fmpz_poly_scalar_mul_mpz(self.x, self.x, T1)
1636
fmpz_poly_scalar_mul_mpz(self.y, self.y, T2)
1637
fmpz_poly_scalar_mul_mpz(self.z, self.z, t3)
1638
fmpz_poly_scalar_mul_mpz(self.w, self.w, t4)
1639
1640
ZZX_to_fmpz_poly(self.a, a.__numerator) # we will assume that the denominator of a and b are 1
1641
ZZX_to_fmpz_poly(self.b, b.__numerator)
1642
1643
ZZX_to_fmpz_poly(self.modulus, (<NumberFieldElement>x).__fld_numerator.x) # and same for the modulus
1644
1645
def __getitem__(self, int i):
1646
"""
1647
EXAMPLES::
1648
1649
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K,-a,a+1)
1650
sage: Q([a,-2/3,a^2-1/2,a*2])
1651
a + (-2/3)*i + (a^2 - 1/2)*j + 2*a*k
1652
sage: x = Q([a,-2/3,a^2-1/2,a*2])
1653
sage: type(x)
1654
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_number_field'>
1655
sage: x[0]
1656
a
1657
sage: x[1]
1658
-2/3
1659
sage: x[2]
1660
a^2 - 1/2
1661
sage: x[3]
1662
2*a
1663
sage: list(x)
1664
[a, -2/3, a^2 - 1/2, 2*a]
1665
"""
1666
# general number -- this code assumes that the number field is not quadratic!!
1667
1668
cdef NumberFieldElement el = <NumberFieldElement>(self._parent.base_ring().an_element())
1669
cdef NumberFieldElement item = el._new()
1670
1671
if i == 0:
1672
fmpz_poly_to_ZZX(item.__numerator, self.x)
1673
elif i == 1:
1674
fmpz_poly_to_ZZX(item.__numerator, self.y)
1675
elif i == 2:
1676
fmpz_poly_to_ZZX(item.__numerator, self.z)
1677
elif i == 3:
1678
fmpz_poly_to_ZZX(item.__numerator, self.w)
1679
else:
1680
raise IndexError, "quaternion element index out of range"
1681
1682
mpz_to_ZZ(&item.__denominator, &self.d)
1683
1684
return item
1685
1686
1687
def __reduce__(self):
1688
"""
1689
EXAMPLES::
1690
1691
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K, -3, a)
1692
sage: z = (i+j+k+a)^2; z
1693
a^2 + 4*a - 3 + 2*a*i + 2*a*j + 2*a*k
1694
sage: f, t = z.__reduce__()
1695
sage: f(*t)
1696
a^2 + 4*a - 3 + 2*a*i + 2*a*j + 2*a*k
1697
sage: loads(dumps(z)) == z
1698
True
1699
"""
1700
return (unpickle_QuaternionAlgebraElement_number_field_v0,
1701
(self._parent, (self[0], self[1], self[2], self[3])))
1702
1703
cpdef ModuleElement _add_(self, ModuleElement _right):
1704
"""
1705
Add self and _right:
1706
1707
EXAMPLES::
1708
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K, -3, a)
1709
sage: z = a + i + (2/3)*a^3*j + (1+a)*k; w = a - i - (2/3)*a^3*j + (1/3+a)*k
1710
sage: type(z)
1711
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_number_field'>
1712
sage: z._add_(w)
1713
2*a + (2*a + 4/3)*k
1714
"""
1715
1716
# Given two quaternion algebra elements
1717
# a = (1/d1)*(x1 + y1 * i + z1 * j + w1 * k)
1718
# b = (1/d2)*(x2 + y2 * i + z2 * j + w2 * k)
1719
#
1720
# we compute their sum as
1721
#
1722
# (a + b) = (1/d3)*(x3 + y3 * i + z3 * j + w3 * k)
1723
#
1724
# with d3 = d1 * d2
1725
# x3 = d1 * x2 + d2 * x1
1726
# y3 = d1 * y2 + d2 * y1
1727
# z3 = d1 * z2 + d2 * z1
1728
# w3 = d1 * w2 + d2 * w1
1729
#
1730
# and then we reduce the sum by calling canonicalize().
1731
1732
# Note: We are assuming in this routine that the modulus is monic. This shouldn't
1733
# currently be an issue because it is impossible to create a number field with
1734
# a modulus that is not monic.
1735
1736
cdef QuaternionAlgebraElement_number_field right = _right
1737
cdef QuaternionAlgebraElement_number_field result = <QuaternionAlgebraElement_number_field> PY_NEW(QuaternionAlgebraElement_number_field)
1738
1739
fmpz_poly_set(result.a, self.a)
1740
fmpz_poly_set(result.b, self.b)
1741
fmpz_poly_set(result.modulus, self.modulus)
1742
result._parent = self._parent
1743
1744
fmpz_poly_scalar_mul_mpz(fU1, self.x, right.d)
1745
fmpz_poly_scalar_mul_mpz(fU2, right.x, self.d)
1746
fmpz_poly_add(result.x, fU1, fU2)
1747
1748
fmpz_poly_scalar_mul_mpz(fU1, self.y, right.d)
1749
fmpz_poly_scalar_mul_mpz(fU2, right.y, self.d)
1750
fmpz_poly_add(result.y, fU1, fU2)
1751
1752
fmpz_poly_scalar_mul_mpz(fU1, self.w, right.d)
1753
fmpz_poly_scalar_mul_mpz(fU2, right.w, self.d)
1754
fmpz_poly_add(result.w, fU1, fU2)
1755
1756
fmpz_poly_scalar_mul_mpz(fU1, self.z, right.d)
1757
fmpz_poly_scalar_mul_mpz(fU2, right.z, self.d)
1758
fmpz_poly_add(result.z, fU1, fU2)
1759
1760
mpz_mul(result.d, self.d, right.d)
1761
1762
self.canonicalize()
1763
1764
return result
1765
1766
1767
1768
1769
1770
cpdef ModuleElement _sub_(self, ModuleElement _right):
1771
"""
1772
Subtract _right from self.
1773
1774
EXAMPLES::
1775
1776
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K, -3, a)
1777
sage: z = a + i + (2/3)*a^3*j + (1+a)*k; w = a - i - (2/3)*a^3*j + (1/3+a)*k
1778
sage: type(z)
1779
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_number_field'>
1780
sage: z._sub_(w)
1781
2*i + 8/3*j + 2/3*k
1782
1783
"""
1784
# Implementation Note: To obtain _sub_, we simply replace every occurrence of
1785
# "add" in _add_ with "sub"; that is, we s/add/sub to get _sub_
1786
1787
# Note: We are assuming in this routine that the modulus is monic. This shouldn't
1788
# currently be an issue because it is impossible to create a number field with
1789
# a modulus that is not monic.
1790
cdef QuaternionAlgebraElement_number_field right = _right
1791
cdef QuaternionAlgebraElement_number_field result = <QuaternionAlgebraElement_number_field> PY_NEW(QuaternionAlgebraElement_number_field)
1792
1793
1794
fmpz_poly_set(result.a, self.a)
1795
fmpz_poly_set(result.b, self.b)
1796
fmpz_poly_set(result.modulus, self.modulus)
1797
result._parent = self._parent
1798
1799
fmpz_poly_scalar_mul_mpz(fU1, self.x, right.d)
1800
fmpz_poly_scalar_mul_mpz(fU2, right.x, self.d)
1801
fmpz_poly_sub(result.x, fU1, fU2)
1802
1803
fmpz_poly_scalar_mul_mpz(fU1, self.y, right.d)
1804
fmpz_poly_scalar_mul_mpz(fU2, right.y, self.d)
1805
fmpz_poly_sub(result.y, fU1, fU2)
1806
1807
fmpz_poly_scalar_mul_mpz(fU1, self.w, right.d)
1808
fmpz_poly_scalar_mul_mpz(fU2, right.w, self.d)
1809
fmpz_poly_sub(result.w, fU1, fU2)
1810
1811
fmpz_poly_scalar_mul_mpz(fU1, self.z, right.d)
1812
fmpz_poly_scalar_mul_mpz(fU2, right.z, self.d)
1813
fmpz_poly_sub(result.z, fU1, fU2)
1814
1815
mpz_mul(result.d, self.d, right.d)
1816
1817
self.canonicalize()
1818
1819
return result
1820
1821
cpdef RingElement _mul_(self, RingElement _right):
1822
"""
1823
Multiply self and _right.
1824
1825
EXAMPLES::
1826
1827
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K, -3, a)
1828
sage: z = a + i + (2/3)*a^3*j + (1+a)*k; w = a - i - (2/3)*a^3*j + (1/3+a)*k
1829
sage: type(z)
1830
<type 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_number_field'>
1831
sage: z._mul_(w)
1832
5*a^2 - 7/9*a + 9 + (-8/3*a^2 - 16/9*a)*i + (-6*a - 4)*j + (2*a^2 + 4/3*a)*k
1833
"""
1834
1835
# We use the following formula for multiplication:
1836
#
1837
# Given two quaternion algebra elements
1838
#
1839
# a = (1/d1)*(x1 + y1 * i + z1 * j + w1 * k)
1840
# b = (1/d2)*(x2 + y2 * i + z2 * j + w2 * k)
1841
#
1842
# we compute their product as
1843
#
1844
# ab = (1/d3)*(x3 + y3 * i + z3 * j + w3 * k)
1845
#
1846
# where
1847
# x3 = t1 + a * t2 + b * (t3 - a*t4)
1848
# y3 = s1*(x2 + y2) - t1 - t2 + b*( s2*(z2 - w2) - t3 + t4)
1849
# z3 = t5 - a*t6 + t7 + a*t8
1850
# w3 = (x2 - y2)*s2 - t5 + t6 + s1*(z2 + w2) - t7 - t8
1851
#
1852
# and where
1853
# t1 = x1 * x2
1854
# t2 = y1 * y2
1855
# t3 = z1 * z2
1856
# t4 = w1 * w2
1857
# t5 = x2 * z1
1858
# t6 = y2 * w1
1859
# t7 = x1 * z2
1860
# t8 = y1 * w2
1861
#
1862
# s1 = x1 + y1
1863
# s2 = z1 + w1
1864
#
1865
# This takes more polynomial addition operations but fewer polynomial multiplication
1866
# operations than the "straightforward" multiplication method.
1867
#
1868
# There might be a way to optimize this formula further.
1869
1870
cdef QuaternionAlgebraElement_number_field right = _right
1871
cdef QuaternionAlgebraElement_number_field result = <QuaternionAlgebraElement_number_field> PY_NEW(QuaternionAlgebraElement_number_field)
1872
1873
mpz_set_si(result.d, 1)
1874
1875
fmpz_poly_set(result.a, self.a)
1876
fmpz_poly_set(result.b, self.b)
1877
fmpz_poly_set(result.modulus, self.modulus)
1878
result._parent = self._parent
1879
1880
fmpz_poly_mul(fT1, self.x, right.x) # t1 = x1 * x2
1881
fmpz_poly_mul(fT2, self.y, right.y) # t2 = y1 * y2
1882
fmpz_poly_mul(ft3, self.z, right.z) # t3 = x1 * x2
1883
fmpz_poly_mul(ft4, self.w, right.w) # t4 = w1 * w2
1884
fmpz_poly_mul(ft5, right.x, self.z) # t5 = x2 * z1
1885
fmpz_poly_mul(ft6, right.y, self.w) # t6 = y2 * w1
1886
fmpz_poly_mul(ft7, self.x, right.z) # t7 = x1 * z2
1887
fmpz_poly_mul(ft8, self.y, right.w) # t8 = y1 * w2
1888
1889
fmpz_poly_add(fs1, self.x, self.y) # s1 = x1 + y1
1890
fmpz_poly_add(fs2, self.z, self.w) # s2 = z1 + w
1891
1892
#------------------
1893
1894
fmpz_poly_mul(fU1, self.a, ft4)
1895
fmpz_poly_sub(fU1, ft3, fU1)
1896
fmpz_poly_mul(fU1, fU1, self.b)
1897
fmpz_poly_mul(fU2, self.a, fT2)
1898
fmpz_poly_add(result.x, fT1, fU2)
1899
fmpz_poly_add(result.x, result.x, fU1)
1900
1901
#------------------
1902
1903
fmpz_poly_sub(fU1, right.z, right.w)
1904
fmpz_poly_mul(fU1, fU1, fs2)
1905
fmpz_poly_sub(fU1, fU1, ft3)
1906
fmpz_poly_add(fU1, fU1, ft4)
1907
fmpz_poly_mul(fU1, fU1, self.b)
1908
fmpz_poly_sub(fU1, fU1, fT2)
1909
fmpz_poly_sub(fU1, fU1, fT1)
1910
fmpz_poly_add(fU2, right.x, right.y)
1911
fmpz_poly_mul(fU2, fs1, fU2)
1912
fmpz_poly_add(result.y, fU1, fU2)
1913
1914
#------------------
1915
1916
fmpz_poly_mul(fU1, self.a, ft8)
1917
fmpz_poly_add(fU1, fU1, ft7)
1918
fmpz_poly_mul(fU2, self.a, ft6)
1919
fmpz_poly_sub(fU1, fU1, fU2)
1920
fmpz_poly_add(result.z, fU1, ft5)
1921
1922
#------------------
1923
1924
fmpz_poly_add(fU1, right.z, right.w)
1925
fmpz_poly_mul(fU1, fU1, fs1)
1926
fmpz_poly_sub(fU1, fU1, ft7)
1927
fmpz_poly_sub(fU1, fU1, ft8)
1928
fmpz_poly_add(fU1, fU1, ft6)
1929
fmpz_poly_sub(fU1, fU1, ft5)
1930
fmpz_poly_sub(fU2, right.x, right.y)
1931
fmpz_poly_mul(fU2, fU2, fs2)
1932
fmpz_poly_add(result.w, fU1, fU2)
1933
1934
# At this point we have essentially computed the product, but we still
1935
# need to reduce modulo the modulus, which is what the following 12 lines do.
1936
#
1937
# When this was written, the version of flint in Sage had problems with
1938
# fpmz_poly_divrem(). This should be fixed in the newest version of
1939
# flint, which also should have some new functions which should do
1940
# this faster (Bill Hart sent an email to Bober and William about this).
1941
#
1942
# This should be fixed in the near future. (I don't know how much
1943
# faster it will be when it is updated, but the following code is
1944
# currently quite a bottleneck.
1945
1946
fmpz_poly_div(fT1, result.x, result.modulus)
1947
fmpz_poly_mul(fT1, fT1, result.modulus)
1948
fmpz_poly_sub(result.x, result.x, fT1)
1949
1950
fmpz_poly_div(fT1, result.y, result.modulus)
1951
fmpz_poly_mul(fT1, fT1, result.modulus)
1952
fmpz_poly_sub(result.y, result.y, fT1)
1953
1954
fmpz_poly_div(fT1, result.z, result.modulus)
1955
fmpz_poly_mul(fT1, fT1, result.modulus)
1956
fmpz_poly_sub(result.z, result.z, fT1)
1957
1958
fmpz_poly_div(fT1, result.w, result.modulus)
1959
fmpz_poly_mul(fT1, fT1, result.modulus)
1960
fmpz_poly_sub(result.w, result.w, fT1)
1961
1962
mpz_mul(result.d, self.d, right.d)
1963
1964
self.canonicalize()
1965
1966
return result
1967
1968
cdef inline canonicalize(self):
1969
"""
1970
Put the representation of this quaternion element into
1971
smallest form. For a = `(1/d)(x + yi + zj + wk)` we
1972
divide `a`, `x`, `y`, `z`, and `w` by the gcd of all of them.
1973
1974
TESTS::
1975
1976
sage: F = QQ[3^(1/3)]
1977
sage: a = F.gen()
1978
sage: K.<i,j,k> = QuaternionAlgebra(F, -10 + a, -7 - a)
1979
sage: ((1/4 + 1/2 * i + a^3/7 * j + a/28 * k)*14*i)^3 # implicit doctest
1980
34503/2*a^2 + 132195/2*a + 791399/4 + (203/8*a^2 - 10591*a + 169225/4)*i + (-84695/4*a^2 + 483413/8*a + 18591/4)*j + (-87/2*a^2 + 18156*a - 72525)*k
1981
"""
1982
1983
# Note: this function changes the module-level global variables
1984
# U1 and U2, so it isn't always safe to use this in the middle of
1985
# another function. Normally this function is called
1986
# at the end of an arithmetic routine, so this is fine.
1987
1988
# Implementation-wise, we compute the GCD's one at a time,
1989
# and quit if it ever becomes one
1990
1991
cdef fmpz_t content = fmpz_init(fmpz_poly_max_limbs(self.x)) # TODO: think about how big this should be (probably the size of d)
1992
# Note that we have to allocate this here, and not
1993
# as a global variable, because fmpz_t's do not
1994
# self allocate memory
1995
fmpz_poly_content(content, self.x)
1996
fmpz_to_mpz(U1, content)
1997
mpz_gcd(U1, self.d, U1)
1998
if mpz_cmp_ui(U1, 1) != 0:
1999
fmpz_poly_content(content, self.y)
2000
fmpz_to_mpz(U2, content)
2001
mpz_gcd(U1, U1, U2)
2002
if mpz_cmp_ui(U1, 1) != 0:
2003
fmpz_poly_content(content, self.z)
2004
fmpz_to_mpz(U2, content)
2005
mpz_gcd(U1, U1, U2)
2006
if mpz_cmp_ui(U1, 1) != 0:
2007
fmpz_poly_content(content, self.w)
2008
fmpz_to_mpz(U2, content)
2009
mpz_gcd(U1, U1, U2)
2010
if mpz_cmp_ui(U1, 1) != 0:
2011
fmpz_poly_scalar_div_mpz(self.x, self.x, U1)
2012
fmpz_poly_scalar_div_mpz(self.y, self.y, U1)
2013
fmpz_poly_scalar_div_mpz(self.z, self.z, U1)
2014
fmpz_poly_scalar_div_mpz(self.w, self.w, U1)
2015
mpz_divexact(self.d, self.d, U1)
2016
2017
fmpz_clear(content)
2018
2019
2020
2021
2022
#######################################################################
2023
# Versioned unpickle functions
2024
#######################################################################
2025
2026
def unpickle_QuaternionAlgebraElement_generic_v0(*args):
2027
"""
2028
EXAMPLES::
2029
2030
sage: K.<X> = QQ[]
2031
sage: Q.<i,j,k> = QuaternionAlgebra(Frac(K), -5,-19); z = 2/3 + i*X - X^2*j + X^3*k
2032
sage: f, t = z.__reduce__()
2033
sage: sage.algebras.quatalg.quaternion_algebra_element.unpickle_QuaternionAlgebraElement_generic_v0(*t)
2034
2/3 + X*i + (-X^2)*j + X^3*k
2035
sage: sage.algebras.quatalg.quaternion_algebra_element.unpickle_QuaternionAlgebraElement_generic_v0(*t) == z
2036
True
2037
"""
2038
return QuaternionAlgebraElement_generic(*args)
2039
2040
def unpickle_QuaternionAlgebraElement_rational_field_v0(*args):
2041
"""
2042
EXAMPLES::
2043
2044
sage: Q.<i,j,k> = QuaternionAlgebra(-5,-19); a = 2/3 + i*5/7 - j*2/5 +19/2
2045
sage: f, t = a.__reduce__()
2046
sage: sage.algebras.quatalg.quaternion_algebra_element.unpickle_QuaternionAlgebraElement_rational_field_v0(*t)
2047
61/6 + 5/7*i - 2/5*j
2048
"""
2049
return QuaternionAlgebraElement_rational_field(*args)
2050
2051
def unpickle_QuaternionAlgebraElement_number_field_v0(*args):
2052
"""
2053
EXAMPLES::
2054
2055
sage: K.<a> = QQ[2^(1/3)]; Q.<i,j,k> = QuaternionAlgebra(K, -3, a); z = i + j
2056
sage: f, t = z.__reduce__()
2057
sage: sage.algebras.quatalg.quaternion_algebra_element.unpickle_QuaternionAlgebraElement_number_field_v0(*t)
2058
i + j
2059
sage: sage.algebras.quatalg.quaternion_algebra_element.unpickle_QuaternionAlgebraElement_number_field_v0(*t) == z
2060
True
2061
"""
2062
return QuaternionAlgebraElement_number_field(*args)
2063
2064