Ticket #13215: trac_13215_skew_polynomials.patch
File trac_13215_skew_polynomials.patch, 254.4 KB (added by , 8 years ago) |
---|
-
module_list.py
# HG changeset patch # User Xavier Caruso <xavier.caruso@normalesup.org> # Date 1350912684 -7200 # Node ID 055730e426275f1a961973a8e22fe55776088242 # Parent 4977d9a6c35a4c78910a846057e861f3e17e248d imported patch skew_polynomial.patch diff --git a/module_list.py b/module_list.py
a b 1795 1795 Extension('sage.rings.polynomial.symmetric_reduction', 1796 1796 sources = ['sage/rings/polynomial/symmetric_reduction.pyx']), 1797 1797 1798 Extension('sage.rings.polynomial.skew_polynomial_element', 1799 sources = ['sage/rings/polynomial/skew_polynomial_element.pyx'], 1800 depends = numpy_depends), 1801 1802 Extension('sage.rings.polynomial.skew_polynomial_finite_field', 1803 sources = ['sage/rings/polynomial/skew_polynomial_finite_field.pyx'], 1804 depends = numpy_depends), 1805 1798 1806 ################################ 1799 1807 ## 1800 1808 ## sage.schemes -
sage/rings/polynomial/all.py
diff --git a/sage/rings/polynomial/all.py b/sage/rings/polynomial/all.py
a b 45 45 # Infinite Polynomial Rings 46 46 from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing 47 47 48 # Skew Polynomial Rings 49 from sage.rings.polynomial.skew_polynomial_ring_constructor import SkewPolynomialRing 50 from sage.rings.polynomial.skew_polynomial_ring import is_SkewPolynomialRing 51 from sage.rings.polynomial.skew_polynomial_element import is_SkewPolynomial 48 52 # Evaluation of cyclotomic polynomials 49 53 from sage.rings.polynomial.cyclotomic import cyclotomic_value -
new file sage/rings/polynomial/skew_polynomial_element.pxd
diff --git a/sage/rings/polynomial/skew_polynomial_element.pxd b/sage/rings/polynomial/skew_polynomial_element.pxd new file mode 100644
- + 1 include "../../ext/interrupt.pxi" 2 include "../../ext/cdefs.pxi" 3 include '../../ext/stdsage.pxi' 4 5 6 from sage.rings.integer cimport Integer 7 8 from sage.structure.element import Element, AlgebraElement 9 from sage.structure.element cimport Element, AlgebraElement, ModuleElement 10 from sage.structure.parent cimport Parent 11 from polynomial_compiled import CompiledPolynomialFunction 12 from polynomial_compiled cimport CompiledPolynomialFunction 13 14 from sage.rings.morphism cimport RingHomomorphism 15 from sage.structure.element cimport RingElement 16 17 from sage.rings.polynomial.polynomial_element cimport Polynomial_generic_dense 18 19 20 cdef class CenterSkewPolynomial_generic_dense(Polynomial_generic_dense): 21 pass 22 23 24 cdef class SkewPolynomial(AlgebraElement): 25 cdef _is_gen 26 27 cdef long _hash_c(self) 28 cdef list _list_c(self) 29 cdef SkewPolynomial _new_c(self,list coeffs,Parent P,char check=*) 30 cpdef SkewPolynomial _new_constant_poly(self,RingElement a,Parent P,char check=*) 31 cpdef ModuleElement _neg_(self) 32 33 34 35 cdef class SkewPolynomial_generic_dense(SkewPolynomial): 36 cdef list __coeffs 37 cdef void __normalize(self) 38 cpdef _pow_(self,right,modulus=*,side=*) 39 40 # Inplace functions 41 cdef void _inplace_add(self, SkewPolynomial_generic_dense right) 42 cdef void _inplace_sub(self, SkewPolynomial_generic_dense right) 43 cdef void _inplace_rmul(self, SkewPolynomial_generic_dense right) 44 cdef void _inplace_lmul(self, SkewPolynomial_generic_dense right) 45 cdef void _inplace_pow(self, Py_ssize_t n) 46 47 48 cdef class SkewPolynomialBaseringInjection(RingHomomorphism): 49 cdef RingElement _an_element 50 cdef object _new_constant_poly_ -
new file sage/rings/polynomial/skew_polynomial_element.pyx
diff --git a/sage/rings/polynomial/skew_polynomial_element.pyx b/sage/rings/polynomial/skew_polynomial_element.pyx new file mode 100644
- + 1 r""" 2 This module implements elements in skew polynomial rings. 3 4 DEFINITION:: 5 6 Let `R` be a commutative ring equipped with an endomorphism `\sigma`. 7 8 The skew polynomial ring over `(R, \sigma)` is the ring `S = `R[X,\sigma]` 9 is the usual ring of polynomials over `R` equipped with the skew 10 multiplication defined by the rule `X*a = \sigma(a)*X` for all `a` 11 in `R`. 12 13 EXAMPLES:: 14 15 We illustrate some functionalities implemented in this class. 16 17 We create the skew polynomial ring:: 18 19 sage: R.<t> = ZZ[] 20 sage: sigma = R.hom([t+1]) 21 sage: S.<x> = R['x',sigma]; S 22 Skew Polynomial Ring in x over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 23 24 and some elements in it:: 25 26 sage: a = t + x + 1; a 27 x + t + 1 28 sage: b = S([t^2,t+1,1]); b 29 x^2 + (t + 1)*x + t^2 30 sage: c = S.random_element(degree=3,monic=True); c 31 x^3 + (-95*t^2 + t + 2)*x^2 + (-t^2 + t)*x + 2*t - 8 32 33 Ring operations are supported:: 34 35 sage: a + b 36 x^2 + (t + 2)*x + t^2 + t + 1 37 sage: a - b 38 -x^2 - t*x - t^2 + t + 1 39 40 sage: a * b 41 x^3 + (2*t + 3)*x^2 + (2*t^2 + 4*t + 2)*x + t^3 + t^2 42 sage: b * a 43 x^3 + (2*t + 4)*x^2 + (2*t^2 + 3*t + 2)*x + t^3 + t^2 44 sage: a * b == b * a 45 False 46 47 sage: b^2 48 x^4 + (2*t + 4)*x^3 + (3*t^2 + 7*t + 6)*x^2 + (2*t^3 + 4*t^2 + 3*t + 1)*x + t^4 49 sage: b^2 == b*b 50 True 51 52 53 Sage implements also some arithmetics over skew polynomial rings. 54 You will find below a short panorama. 55 56 DEFINITION: 57 58 Let `a` and `b` be two skew polynomials in the same ring `S`. 59 The *left (resp. right) euclidean division* of `a` by `b` is 60 a couple `(q,r)` of elements in `S` such that 61 62 - `a = q*b + r` (resp. `a = b*q + r`) 63 64 - the degree of `r` is less than the degree of `b` 65 66 `q` (resp. `r`) is called the *quotient* (resp. the remainder) 67 of this euclidean division. 68 69 PROPERTIES: 70 71 Keeping the previous notations, if the leading coefficient of `b` 72 is a unit (e.g. if `b` is monic) then the quotient and the remainder 73 in the *right* euclidean division exist and are unique. 74 75 The same result holds for the *left* euclidean division if in addition 76 the twist map defining the skew polynomial ring is invertible. 77 78 EXAMPLES:: 79 80 sage: q,r = c.quo_rem(b) # default side is right 81 sage: q 82 x - 95*t^2 83 sage: r 84 (95*t^3 + 93*t^2 - t - 1)*x + 95*t^4 + 2*t - 8 85 sage: c == q*b + r 86 True 87 88 The operators ``//`` and ``%`` give respectively the quotient 89 and the remainder of the right euclidean division:: 90 91 sage: q == c // b 92 True 93 sage: r == c % b 94 True 95 96 If we want left euclidean division, we need to specify 97 ``side=Left``. Nonetheless, it won't work over our current 98 `S` because Sage can't invert the twist map:: 99 100 sage: q,r = c.quo_rem(b,side=Left) 101 Traceback (most recent call last): 102 ... 103 NotImplementedError 104 105 Here is a working example over a finite field:: 106 107 sage: k.<t> = GF(5^3) 108 sage: Frob = k.frobenius_endomorphism() 109 sage: S.<x> = k['x',Frob] 110 sage: a = x^4 + (4*t + 1)*x^3 + (t^2 + 3*t + 3)*x^2 + (3*t^2 + 2*t + 2)*x + 3*t^2 + 3*t + 1 111 sage: b = (2*t^2 + 3)*x^2 + (3*t^2 + 1)*x + 4*t + 2 112 sage: q,r = a.quo_rem(b,side=Left) 113 sage: q 114 (4*t^2 + t + 1)*x^2 + (2*t^2 + 2*t + 2)*x + 2*t^2 + 4*t + 3 115 sage: r 116 (t + 2)*x + 3*t^2 + 2*t + 4 117 sage: a == b*q + r 118 True 119 120 121 Once we have euclidean divisions, we have for free gcd and lcm 122 (at least if the base ring is a field). 123 This class provides an implementation of gcd and lcm:: 124 125 sage: a = (x + t) * (x + t^2)^2 126 sage: b = (x + t) * (t*x + t + 1) * (x + t^2) 127 sage: a.gcd(b) # default side is right 128 x + t^2 129 sage: a.gcd(b,side=Left) 130 x + t 131 132 For lcm, the default side is left but be very careful: by 133 convention, a left (resp. right) lcm is common multiple on 134 the right (resp. left):: 135 136 sage: c = a.lcm(b); c # default side is left 137 x^5 + (4*t^2 + t + 3)*x^4 + (3*t^2 + 4*t)*x^3 + 2*t^2*x^2 + (2*t^2 + t)*x + 4*t^2 + 4 138 sage: c.is_divisible_by(a) 139 True 140 sage: c.is_divisible_by(b) 141 True 142 143 sage: d = a.lcm(b,side=Right); d 144 x^5 + (t^2 + 1)*x^4 + (3*t^2 + 3*t + 3)*x^3 + (3*t^2 + t + 2)*x^2 + (4*t^2 + 3*t)*x + 4*t + 4 145 sage: d.is_divisible_by(a,side=Left) 146 True 147 sage: d.is_divisible_by(b,side=Left) 148 True 149 150 AUTHOR: 151 152 - Xavier Caruso (2012-06-29) 153 """ 154 155 ############################################################################# 156 # Copyright (C) 2012 Xavier Caruso <xavier.caruso@normalesup.org> 157 # 158 # Distributed under the terms of the GNU General Public License (GPL) 159 # 160 # http://www.gnu.org/licenses/ 161 #**************************************************************************** 162 163 164 include "../../ext/stdsage.pxi" 165 166 import operator, copy, re 167 168 import skew_polynomial_ring 169 import sage.rings.infinity as infinity 170 from sage.misc.latex import latex 171 from sage.structure.factorization import Factorization 172 173 from sage.categories.homset import Hom 174 175 from sage.structure.element import RingElement 176 from sage.structure.element cimport Element, RingElement, ModuleElement 177 178 from sage.rings.ring import Field 179 180 from sage.structure.parent_gens cimport ParentWithGens 181 182 from sage.rings.integer cimport Integer 183 from sage.categories.map cimport Map 184 from sage.rings.morphism cimport RingHomomorphism 185 from sage.rings.polynomial.polynomial_element cimport Polynomial_generic_dense 186 187 from sage.structure.side import Left, Right 188 189 def is_SkewPolynomial(a): 190 """ 191 Return True if `a` is a skew polynomial (over some base). 192 193 INPUT: 194 195 - ``a`` -- an object 196 197 EXAMPLES:: 198 199 sage: from sage.rings.polynomial.skew_polynomial_element import is_SkewPolynomial 200 sage: R.<t> = ZZ[] 201 sage: sigma = R.hom([t+1]) 202 sage: S.<x> = R['x',sigma] 203 sage: a = x^2 204 sage: is_SkewPolynomial(a) 205 True 206 """ 207 return PY_TYPE_CHECK(a, SkewPolynomial) 208 209 210 cdef class CenterSkewPolynomial_generic_dense(Polynomial_generic_dense): 211 """ 212 A class for elements in the center of a skew polynomial ring. 213 """ 214 pass 215 216 217 cdef class SkewPolynomial(AlgebraElement): 218 """ 219 A skew polynomial. 220 """ 221 def __init__(self,parent,is_gen=False,construct=False): 222 """ 223 The following examples illustrate the creation of elements of 224 skew polynomial rings. 225 226 sage: R.<t> = ZZ[] 227 sage: sigma = R.hom([t+1]) 228 sage: S.<x> = R['x',sigma] 229 sage: P = x+t; P 230 x + t 231 sage: Q = S([1,t,t+2]); Q 232 (t + 2)*x^2 + t*x + 1 233 """ 234 AlgebraElement.__init__(self, parent) 235 self._is_gen = is_gen 236 237 # you may have to replicate this boilerplate code in derived classes if you override 238 # __richcmp__. The python documentation at http://docs.python.org/api/type-structs.html 239 # explains how __richcmp__, __hash__, and __cmp__ are tied together. 240 def __hash__(self): 241 return self._hash_c() 242 243 cdef long _hash_c(self): 244 """ 245 This hash incorporates the name of the variable. 246 """ 247 cdef long result = 0 # store it in a c-int and just let the overflowing additions wrap 248 cdef long result_mon 249 cdef long c_hash 250 cdef long var_name_hash 251 cdef int i 252 for i from 0<= i <= self.degree(): 253 if i == 1: 254 # we delay the hashing until now to not waste it one a constant poly 255 var_name_hash = hash((<ParentWithGens>self._parent)._names[0]) 256 # I'm assuming (incorrectly) that hashes of zero indicate that the element is 0. 257 # This assumption is not true, but I think it is true enough for the purposes and it 258 # it allows us to write fast code that omits terms with 0 coefficients. This is 259 # important if we want to maintain the '==' relationship with sparse polys. 260 c_hash = hash(self[i]) 261 if c_hash != 0: 262 if i == 0: 263 result += c_hash 264 else: 265 # Hash (self[i], generator, i) as a tuple according to the algorithm. 266 result_mon = c_hash 267 result_mon = (1000003 * result_mon) ^ var_name_hash 268 result_mon = (1000003 * result_mon) ^ i 269 result += result_mon 270 if result == -1: 271 return -2 272 return result 273 274 275 # Comparison 276 # ---------- 277 278 def __richcmp__(left, right, int op): 279 return (<Element>left)._richcmp(right,op) 280 281 282 cdef int _cmp_c_impl(self, Element other) except -2: 283 """ 284 Compare the two skew polynomials self and other. 285 286 We order polynomials first by degree, then in dictionary order 287 starting with the coefficient of largest degree. 288 """ 289 cdef list x = (<SkewPolynomial>self)._list_c() 290 cdef list y = (<SkewPolynomial>other)._list_c() 291 cdef Py_ssize_t dx = len(x), dy = len(y) 292 cdef Py_ssize_t i 293 cdef int c 294 c = cmp(dx,dy) 295 if c: return c 296 for i from dx > i >= 0: 297 c = cmp(x[i],y[i]) 298 if c: return c 299 return 0 300 301 302 # Some c functions 303 # ---------------- 304 305 cdef SkewPolynomial _new_c(self,list coeffs,Parent P,char check=0): 306 """ 307 Fast creation of a new skew polynomial 308 309 .. NOTE:: 310 311 Override this function in classes which inherit 312 from SkewPolynomial. 313 """ 314 return P(list) 315 316 317 cpdef SkewPolynomial _new_constant_poly(self,RingElement a,Parent P,char check=0): 318 """ 319 Fast creation of a new constant skew polynomial 320 """ 321 if a: 322 return self._new_c([a],P,check) 323 else: 324 return self._new_c([],P) 325 326 327 cdef list _list_c(self): 328 """ 329 Return the list of the underlying elements of this 330 skew polynomial. 331 332 .. WARNING:: 333 334 It is *not* a copy; do not modify this list! 335 """ 336 raise NotImplementedError 337 338 339 # A skew polynomial as a list of coefficients 340 # ------------------------------------------- 341 342 def list(self): 343 """ 344 Return a new copy of the list of the underlying elements of self. 345 346 EXAMPLES:: 347 348 sage: R.<t> = QQ[] 349 sage: sigma = R.hom([t+1]) 350 sage: S.<x> = R['x',sigma] 351 sage: a = 1 + x^4 + (t+1)*x^2 + t^2 352 sage: l = a.list(); l 353 [t^2 + 1, 0, t + 1, 0, 1] 354 355 Note that v is a list, it is mutable, and each call to the list 356 method returns a new list:: 357 358 sage: type(l) 359 <type 'list'> 360 sage: l[0] = 5 361 sage: a.list() 362 [t^2 + 1, 0, t + 1, 0, 1] 363 """ 364 return list((<SkewPolynomial>self)._list_c()) 365 366 367 def __getitem__(self,n): 368 """ 369 INPUT:: 370 371 - ``n`` -- an integer 372 373 OUTPUT:: 374 375 - the `n`-th coefficient of self 376 377 .. NOTE:: 378 379 Coefficients are on the left (ie before the variable) 380 381 EXAMPLES:: 382 383 sage: R.<t> = QQ[] 384 sage: sigma = R.hom([t+1]) 385 sage: S.<x> = R['x',sigma] 386 sage: a = t*x^2 + (t + 3/7)*x + t^2 387 sage: a[1] 388 t + 3/7 389 sage: a[3] 390 0 391 """ 392 try: 393 return (<SkewPolynomial>self)._list_c()[n] 394 except IndexError: 395 return self.base_ring()(0) 396 397 398 def __getslice__(self, Py_ssize_t i, Py_ssize_t j): 399 """ 400 EXAMPLES:: 401 sage: R.<t> = QQ[] 402 sage: sigma = R.hom([t+1]) 403 sage: S.<x> = R['x',sigma] 404 sage: a = t*x^2 + (t + 3/7)*x + t^2 405 sage: a[1:] 406 t*x^2 + (t + 3/7)*x 407 sage: a[:1] 408 t^2 409 sage: a[3:] 410 0 411 """ 412 if i <= 0: 413 i = 0 414 zeros = [] 415 elif i > 0: 416 zeros = [self._parent.base_ring()(0)] * i 417 return self._new_c(zeros + self._list_c()[i:j], self._parent, 1) 418 419 420 def __setitem__(self, n, value): 421 """ 422 Set the `n`-th coefficient of this skew polynomial. This always 423 raises an IndexError, since in Sage polynomials are immutable. 424 425 INPUT: 426 427 - ``n`` - an integer 428 429 - ``value`` - value to set the `n`-th coefficient to 430 431 OUTPUT: an IndexError is always raised. 432 433 EXAMPLES:: 434 435 sage: k.<t> = GF(5^3) 436 sage: Frob = k.frobenius_endomorphism() 437 sage: S.<x> = k['x',Frob] 438 sage: a = x + t 439 sage: a[1] = t + 1 440 Traceback (most recent call last): 441 ... 442 IndexError: skew polynomials are immutable 443 """ 444 raise IndexError, "skew polynomials are immutable" 445 446 447 # Basic operations 448 # ---------------- 449 450 def degree(self): 451 """ 452 Return the degree of this skew polynomial. The zero 453 skew polynomial has degree -1. 454 455 EXAMPLES:: 456 457 sage: R.<t> = ZZ[] 458 sage: sigma = R.hom([t+1]) 459 sage: S.<x> = R['x',sigma] 460 sage: a = x^2 + t*x^3 + t^2*x + 1 461 sage: a.degree() 462 3 463 464 By convention, the degree of 0 is -1:: 465 466 sage: S(0).degree() 467 -1 468 """ 469 return len((<SkewPolynomial>self)._list_c())-1 470 471 472 def valuation(self): 473 """ 474 Return the valuation of this skew polynomial. 475 The zero skew polynomial has valuation +Infinity. 476 477 EXAMPLES:: 478 479 sage: R.<t> = ZZ[] 480 sage: sigma = R.hom([t+1]) 481 sage: S.<x> = R['x',sigma] 482 sage: a = x^2 + t*x^3 + t^2*x 483 sage: a.valuation() 484 1 485 486 By convention, the valuation of 0 is +Infinity:: 487 488 sage: S(0).valuation() 489 +Infinity 490 """ 491 cdef list x = (<SkewPolynomial>self)._list_c() 492 if len(x) == 0: 493 return infinity.infinity 494 cdef Py_ssize_t v = 0 495 while x[v].is_zero(): 496 v += 1 497 return v 498 499 500 cpdef ModuleElement _add_(self, ModuleElement right): 501 """ 502 Compute self + right 503 504 INPUT: 505 506 - right -- a skew polynomial over the same base 507 508 TESTS:: 509 510 sage: R.<t> = QQ[] 511 sage: sigma = R.hom([t+1]) 512 sage: S.<x> = R['x',sigma] 513 sage: a = S.random_element(monic=True); a 514 x^2 + (-12*t^2 + 1/2*t - 1/95)*x - 1/2*t^2 - 4 515 sage: b = -S.random_element(monic=True); b 516 -x^2 + (5/2*t - 2/3)*x + 1/4*t^2 - 1/2*t + 1 517 sage: c = a+b; c 518 (-12*t^2 + 3*t - 193/285)*x - 1/4*t^2 - 1/2*t - 3 519 sage: c.degree() 520 1 521 """ 522 cdef Py_ssize_t i, min 523 cdef list x = (<SkewPolynomial>self)._list_c() 524 cdef list y = (<SkewPolynomial>right)._list_c() 525 cdef Py_ssize_t dx = len(x), dy = len(y) 526 527 if dx > dy: 528 return self._new_c([x[i] + y[i] for i from 0 <= i < dy] + x[dy:], self._parent, 0) 529 elif dx < dy: 530 return self._new_c([x[i] + y[i] for i from 0 <= i < dx] + y[dx:], self._parent, 0) 531 else: 532 return self._new_c([x[i] + y[i] for i from 0 <= i < dx], self._parent, 1) 533 534 535 cpdef ModuleElement _sub_(self, ModuleElement right): 536 """ 537 Compute self - right 538 539 INPUT: 540 541 - right -- a skew polynomial over the same base 542 543 TESTS:: 544 545 sage: R.<t> = QQ[] 546 sage: sigma = R.hom([t+1]) 547 sage: S.<x> = R['x',sigma] 548 sage: a = S.random_element(monic=True); a 549 x^2 + (-12*t^2 + 1/2*t - 1/95)*x - 1/2*t^2 - 4 550 sage: b = S.random_element(monic=True); b 551 x^2 + (-5/2*t + 2/3)*x - 1/4*t^2 + 1/2*t - 1 552 sage: c = a-b; c 553 (-12*t^2 + 3*t - 193/285)*x - 1/4*t^2 - 1/2*t - 3 554 sage: c.degree() 555 1 556 """ 557 cdef Py_ssize_t i, min 558 cdef list x = (<SkewPolynomial>self)._list_c() 559 cdef list y = (<SkewPolynomial>right)._list_c() 560 cdef Py_ssize_t dx = len(x), dy = len(y) 561 cdef RingElement c 562 563 if dx > dy: 564 return self._new_c([x[i] - y[i] for i from 0 <= i < dy] + x[dy:], self._parent, 0) 565 elif dx < dy: 566 return self._new_c([x[i] - y[i] for i from 0 <= i < dx] + [ -c for c in y[dx:] ], self._parent, 0) 567 else: 568 return self._new_c([x[i] - y[i] for i from 0 <= i < dx], self._parent, 1) 569 570 571 cpdef ModuleElement _neg_(self): 572 """ 573 Return the opposite of self 574 575 TESTS:: 576 577 sage: R.<t> = QQ[] 578 sage: sigma = R.hom([t+1]) 579 sage: S.<x> = R['x',sigma] 580 sage: a = t*x^2 + x - 3 581 sage: -a 582 -t*x^2 - x + 3 583 """ 584 return self._new_c([-x for x in (<SkewPolynomial>self)._list_c()], self._parent, 0) 585 586 587 cpdef ModuleElement _lmul_(self, RingElement right): 588 """ 589 Compute self * right (in this order) 590 591 INPUT: 592 593 - right -- an element of the base ring 594 595 TESTS:: 596 597 sage: R.<t> = QQ[] 598 sage: sigma = R.hom([t+1]) 599 sage: S.<x> = R['x',sigma] 600 sage: a = x + t 601 sage: b = t 602 sage: a * b 603 (t + 1)*x + t^2 604 sage: a * b == b * a 605 False 606 """ 607 if right == 0: 608 return self._parent(0) 609 cdef list x = (<SkewPolynomial>self)._list_c() 610 cdef Py_ssize_t i 611 map = self._parent._map 612 return self._new_c([ (map**i)(right)*x[i] for i from 0 <= i < len(x) ], self._parent, 0) 613 614 615 cpdef ModuleElement _rmul_(self, RingElement left): 616 """ 617 Compute left * self (in this order) 618 619 INPUT: 620 621 - left -- an element of the base ring 622 623 TESTS:: 624 625 sage: R.<t> = QQ[] 626 sage: sigma = R.hom([t+1]) 627 sage: S.<x> = R['x',sigma] 628 sage: a = t 629 sage: b = x + t 630 sage: a * b 631 t*x + t^2 632 sage: a * b == b * a 633 False 634 """ 635 if left == 0: 636 return self.parent().zero_element() 637 cdef list x = (<SkewPolynomial>self)._list_c() 638 cdef Py_ssize_t i 639 map = self._parent._map 640 return self._new_c([ left*x[i] for i from 0 <= i < len(x) ], self._parent, 0) 641 642 643 cpdef RingElement _mul_(self, RingElement right): 644 """ 645 Compute self * right (in this order) 646 647 INPUT: 648 649 - right -- a skew polynomial in the same ring 650 651 TESTS:: 652 653 sage: R.<t> = QQ[] 654 sage: sigma = R.hom([t+1]) 655 sage: S.<x> = R['x',sigma] 656 sage: a = x^2 + t; a 657 x^2 + t 658 sage: b = x^2 + (t + 1)*x; b 659 x^2 + (t + 1)*x 660 sage: a * b 661 x^4 + (t + 3)*x^3 + t*x^2 + (t^2 + t)*x 662 sage: a * b == b * a 663 False 664 """ 665 cdef list x = (<SkewPolynomial>self)._list_c() 666 cdef list y = (<SkewPolynomial>right)._list_c() 667 cdef Py_ssize_t i, k, start, end 668 cdef Py_ssize_t dx = len(x)-1, dy = len(y)-1 669 parent = self._parent 670 if dx == -1: 671 return self 672 elif dy == -1: 673 return right 674 elif dx == 0: 675 c = x[0] 676 return self._new_c([c*a for a in y], parent, 0) 677 cdef list coeffs = [] 678 for k from 0 <= k <= dx+dy: 679 start = 0 if k <= dy else k-dy # max(0, k-dy) 680 end = k if k <= dx else dx # min(k, dx) 681 sum = x[start] * parent.twist_map(start)(y[k-start]) 682 for i from start < i <= end: 683 sum += x[i] * parent.twist_map(i)(y[k-i]) 684 coeffs.append(sum) 685 return self._new_c(coeffs, parent, 0) 686 687 688 def square(self): 689 """ 690 Return the square of self 691 692 TESTS:: 693 694 sage: R.<t> = QQ[] 695 sage: sigma = R.hom([t+1]) 696 sage: S.<x> = R['x',sigma] 697 sage: a = x + t; a 698 x + t 699 sage: a.square() 700 x^2 + (2*t + 1)*x + t^2 701 sage: a.square() == a*a 702 True 703 """ 704 return self * self 705 706 707 def conjugate(self,n): 708 """ 709 INPUT: 710 711 - `n` -- an integer 712 713 OUTPUT: 714 715 - this skew polynomial conjugated by x^n (where x is 716 the variable) 717 718 .. NOTE:: 719 720 This conjugate is obtained from the skew polynomial by 721 applying the `n`-th iterate of the twist map to each of 722 its coefficients. 723 724 EXAMPLES:: 725 726 sage: R.<t> = QQ[] 727 sage: sigma = R.hom([t+1]) 728 sage: S.<x> = R['x',sigma] 729 sage: a = t*x^3 + (t^2 + 1)*x^2 + 2*t 730 sage: b = a.conjugate(2); b 731 (t + 2)*x^3 + (t^2 + 4*t + 5)*x^2 + 2*t + 4 732 sage: x^2*a == b*x^2 733 True 734 735 In principle, negative values for `n` are allowed... but Sage 736 needs to be able to invert the twist map:: 737 738 sage: b = a.conjugate(-1) 739 Traceback (most recent call last): 740 ... 741 NotImplementedError 742 743 Here is a working example:: 744 745 sage: k.<t> = GF(5^3) 746 sage: Frob = k.frobenius_endomorphism() 747 sage: T.<y> = k['y',Frob] 748 sage: u = T.random_element(); u 749 (2*t^2 + 3)*y^2 + (4*t^2 + t + 4)*y + 2*t^2 + 2 750 sage: v = u.conjugate(-1); v 751 (3*t^2 + t)*y^2 + (4*t^2 + 2*t + 4)*y + 3*t^2 + t + 4 752 sage: u*y == y*v 753 True 754 """ 755 return self._new_c([ self._parent.twist_map(n)(x) for x in (<SkewPolynomial>self)._list_c() ], self._parent, 0) 756 757 758 # Other useful mathematical functions 759 # ----------------------------------- 760 761 def constant_coefficient(self): 762 """ 763 Return the constant coefficient (i.e. the coefficient of degree 764 `0`) of this skew polynomial. 765 766 OUTPUT: element of base ring 767 768 EXAMPLES:: 769 770 sage: R.<t> = ZZ[] 771 sage: sigma = R.hom([t+1]) 772 sage: S.<x> = R['x',sigma] 773 sage: a = x + t^2 + 2 774 sage: a.constant_coefficient() 775 t^2 + 2 776 """ 777 cdef x = (<SkewPolynomial>self)._list_c() 778 if len(x) == 0: 779 return self.base_ring()(0) 780 else: 781 return x[0] 782 783 784 def leading_coefficient(self): 785 """ 786 Return the leading coefficient of this skew polynomial. 787 788 OUTPUT: element of the base ring 789 790 EXAMPLES:: 791 792 sage: R.<t> = ZZ[] 793 sage: sigma = R.hom([t+1]) 794 sage: S.<x> = R['x',sigma] 795 sage: a = x + (t+1)*x^5 + t^2*x^3 - x^5 796 sage: a.leading_coefficient() 797 t 798 """ 799 cdef x = (<SkewPolynomial>self)._list_c() 800 if len(x) == 0: 801 raise ValueError("") 802 return x[-1] 803 804 805 def __nonzero__(self): 806 return self.degree() >= 0 807 808 809 def is_unit(self): 810 """ 811 Return True if this skew polynomial is a unit. 812 813 Not yet implemented. 814 """ 815 raise NotImplementedError 816 817 818 def is_monic(self): 819 """ 820 Returns True if this skew polynomial is monic. The zero polynomial 821 is by definition not monic. 822 823 EXAMPLES:: 824 825 sage: R.<t> = ZZ[] 826 sage: sigma = R.hom([t+1]) 827 sage: S.<x> = R['x',sigma] 828 sage: a = x + t 829 sage: a.is_monic() 830 True 831 sage: a = 0*x 832 sage: a.is_monic() 833 False 834 sage: a = t*x^3 + x^4 + (t+1)*x^2 835 sage: a.is_monic() 836 True 837 sage: a = (t^2 + 2*t)*x^2 + x^3 + t^10*x^5 838 sage: a.is_monic() 839 False 840 """ 841 return not self.is_zero() and self[self.degree()] == 1 842 843 844 def lmonic(self): 845 """ 846 Return the unique monic skew polynomial `a` of the same 847 degree which divides this skew polynomial on the left 848 849 .. Note:: 850 851 This skew polynomial is self dividing on the 852 *right* by the `n`-th iterative (`n` is the degree of 853 self) of the inverse of the twist map applied to the 854 leading coefficient. 855 856 EXAMPLES:: 857 858 sage: k.<t> = GF(5^3) 859 sage: Frob = k.frobenius_endomorphism() 860 sage: S.<x> = k['x',Frob] 861 sage: a = (3*t^2 + 3*t + 2)*x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 862 sage: b = a.lmonic(); b 863 x^3 + (4*t^2 + 3*t)*x^2 + (4*t + 2)*x + 2*t^2 + 4*t + 3 864 865 Check list:: 866 867 sage: b.degree() == a.degree() 868 True 869 sage: b.is_divisible_by(a,side=Left) 870 True 871 sage: twist = S.twist_map(-a.degree()) 872 sage: a == b * twist(a.leading_coefficient()) 873 True 874 875 Note that `b` does not divise `a` on the right:: 876 877 sage: b.is_divisible_by(a,side=Right) 878 False 879 880 This function does not work if the leading coefficient is not a 881 unit:: 882 883 sage: R.<t> = QQ[] 884 sage: sigma = R.hom([t+1]) 885 sage: S.<x> = R['x',sigma] 886 sage: a = t*x 887 sage: a.lmonic() 888 Traceback (most recent call last): 889 ... 890 NotImplementedError: the leading coefficient is not a unit 891 """ 892 try: 893 a = self.base_ring()(~self.leading_coefficient()) 894 except (ZeroDivisionError, TypeError): 895 raise NotImplementedError("the leading coefficient is not a unit") 896 return self*self._parent.twist_map(-self.degree())(a) 897 898 899 def rmonic(self): 900 """ 901 Return the unique monic skew polynomial `a` of the same 902 degree which divides this skew polynomial on the right 903 904 .. Note:: 905 906 This skew polynomial is self dividing on the *left* 907 by its leading coefficient. 908 909 EXAMPLES:: 910 911 sage: k.<t> = GF(5^3) 912 sage: Frob = k.frobenius_endomorphism() 913 sage: S.<x> = k['x',Frob] 914 sage: a = (3*t^2 + 3*t + 2)*x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 915 sage: b = a.rmonic(); b 916 x^3 + (2*t^2 + 3*t + 4)*x^2 + (3*t^2 + 4*t + 1)*x + 2*t^2 + 4*t + 3 917 918 Check list:: 919 920 sage: b.degree() == a.degree() 921 True 922 sage: b.is_divisible_by(a,side=Right) 923 True 924 sage: a == a.leading_coefficient() * b 925 True 926 927 Note that `b` does not divise `a` on the right:: 928 929 sage: b.is_divisible_by(a,side=Left) 930 False 931 932 This function does not work if the leading coefficient is not a 933 unit:: 934 935 sage: R.<t> = QQ[] 936 sage: sigma = R.hom([t+1]) 937 sage: S.<x> = R['x',sigma] 938 sage: a = t*x 939 sage: a.rmonic() 940 Traceback (most recent call last): 941 ... 942 NotImplementedError: the leading coefficient is not a unit 943 """ 944 try: 945 a = self.base_ring()(~self.leading_coefficient()) 946 except (ZeroDivisionError, TypeError): 947 raise NotImplementedError("the leading coefficient is not a unit") 948 return a*self 949 950 951 def monic(self, side=Right): 952 """ 953 INPUT: 954 955 - ``side`` -- ``Left`` or ``Right`` (default: Right) 956 957 OUTPUT: 958 959 - the unique monic skew polynomial `a` of the same 960 degree which divides this skew polynomial on ``side`` 961 962 .. Note:: 963 964 if ``side`` is ``Right``, this skew polynomial is self 965 dividing on the *left* by its leading coefficient; if 966 ``side`` is ``Left``, it is self dividing on the 967 *right* by the `n`-th iterative (`n` is the degree of 968 self) of the inverse of the twist map applied to the 969 leading coefficient. 970 971 EXAMPLES:: 972 973 sage: k.<t> = GF(5^3) 974 sage: Frob = k.frobenius_endomorphism() 975 sage: S.<x> = k['x',Frob] 976 sage: a = (3*t^2 + 3*t + 2)*x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 977 sage: b = a.monic(); b 978 x^3 + (2*t^2 + 3*t + 4)*x^2 + (3*t^2 + 4*t + 1)*x + 2*t^2 + 4*t + 3 979 980 Check list:: 981 982 sage: b.degree() == a.degree() 983 True 984 sage: b.is_divisible_by(a,side=Right) 985 True 986 sage: a == a.leading_coefficient() * b 987 True 988 989 Note that `b` does not divise `a` on the left:: 990 991 sage: b.is_divisible_by(a,side=Left) 992 False 993 994 The same on the left:: 995 996 sage: b = a.monic(side=Left); b 997 x^3 + (4*t^2 + 3*t)*x^2 + (4*t + 2)*x + 2*t^2 + 4*t + 3 998 sage: b.degree() == a.degree() 999 True 1000 sage: b.is_divisible_by(a,side=Left) 1001 True 1002 sage: twist = S.twist_map(-a.degree()) 1003 sage: a == b * twist(a.leading_coefficient()) 1004 True 1005 1006 This function does not work if the leading coefficient is not a 1007 unit:: 1008 1009 sage: R.<t> = QQ[] 1010 sage: sigma = R.hom([t+1]) 1011 sage: S.<x> = R['x',sigma] 1012 sage: a = t*x 1013 sage: a.monic() 1014 Traceback (most recent call last): 1015 ... 1016 NotImplementedError: the leading coefficient is not a unit 1017 """ 1018 if self.is_monic(): 1019 return self 1020 if side == Right: 1021 return self.rmonic() 1022 else: 1023 return self.lmonic() 1024 1025 1026 # Divisibility 1027 # ------------ 1028 1029 1030 def lquo_rem(self, other): 1031 """ 1032 INPUT: 1033 1034 - ``other`` -- a skew polynomial ring over the same 1035 base ring 1036 1037 OUTPUT: 1038 1039 - the quotient and the remainder of the left euclidean 1040 division of this skew polynomial by ``other`` 1041 1042 .. NOTE:: 1043 1044 Doesn't work if the leading coefficient of the divisor 1045 is not a unit or if Sage can't invert the twist map. 1046 1047 EXAMPLES:: 1048 1049 sage: k.<t> = GF(5^3) 1050 sage: Frob = k.frobenius_endomorphism() 1051 sage: S.<x> = k['x',Frob] 1052 sage: a = (3*t^2 + 3*t + 2)*x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 1053 sage: b = (3*t^2 + 4*t + 2)*x^2 + (2*t^2 + 4*t + 3)*x + 2*t^2 + t + 1 1054 sage: q,r = a.lquo_rem(b) 1055 sage: a == b*q + r 1056 True 1057 1058 In the following example, Sage does not know the inverse 1059 of the twist map:: 1060 1061 sage: R.<t> = ZZ[] 1062 sage: sigma = R.hom([t+1]) 1063 sage: S.<x> = R['x',sigma] 1064 sage: a = (-2*t^2 - t + 1)*x^3 + (-t^2 + t)*x^2 + (-12*t - 2)*x - t^2 - 95*t + 1 1065 sage: b = x^2 + (5*t - 6)*x - 4*t^2 + 4*t - 1 1066 sage: a.lquo_rem(b) 1067 Traceback (most recent call last): 1068 ... 1069 NotImplementedError 1070 """ 1071 cdef list a = list((<SkewPolynomial>self)._list_c()) 1072 cdef list b = (<SkewPolynomial>other)._list_c() 1073 cdef Py_ssize_t i, j 1074 cdef Py_ssize_t da = self.degree(), db = other.degree() 1075 if db < 0: 1076 raise ZeroDivisionError 1077 if da < db: 1078 return self._new_c([],self._parent), self 1079 try: 1080 inv = self.base_ring()(~b[db]) 1081 except (ZeroDivisionError, TypeError): 1082 raise NotImplementedError("the leading coefficient of the divisor is not invertible") 1083 cdef list q = [ ] 1084 parent = self._parent 1085 for i from da-db >= i >= 0: 1086 c = parent.twist_map(-db)(inv*a[i+db]) 1087 for j from 0 <= j < db: 1088 a[i+j] -= b[j] * parent.twist_map(j)(c) 1089 q.append(c) 1090 q.reverse() 1091 return self._new_c(q,parent), self._new_c(a[:db],parent,1) 1092 1093 1094 def rquo_rem(self, other): 1095 """ 1096 INPUT: 1097 1098 - ``other`` -- a skew polynomial ring over the same 1099 base ring 1100 1101 OUTPUT: 1102 1103 - the quotient and the remainder of the left euclidean 1104 division of this skew polynomial by ``other`` 1105 1106 .. NOTE:: 1107 1108 Doesn't work if the leading coefficient of the divisor 1109 is not a unit. 1110 1111 EXAMPLES:: 1112 1113 sage: R.<t> = ZZ[] 1114 sage: sigma = R.hom([t+1]) 1115 sage: S.<x> = R['x',sigma] 1116 sage: a = S.random_element(degree=4); a 1117 t^2*x^4 + (-12*t^2 - 2*t - 1)*x^3 + (-95*t^2 + t + 2)*x^2 + (-t^2 + t)*x + 2*t - 8 1118 sage: b = S.random_element(monic=True); b 1119 x^2 + (4*t^2 - t - 2)*x - t^2 + t - 1 1120 sage: q,r = a.rquo_rem(b) 1121 sage: a == q*b + r 1122 True 1123 1124 The leading coefficient of the divisor need to be invertible:: 1125 1126 sage: c = S.random_element(); c 1127 (-4*t^2 + t)*x^2 - 2*t^2*x + 5*t^2 - 6*t - 4 1128 sage: a.rquo_rem(c) 1129 Traceback (most recent call last): 1130 ... 1131 NotImplementedError: the leading coefficient of the divisor is not invertible 1132 """ 1133 cdef list a = list((<SkewPolynomial>self)._list_c()) 1134 cdef list b = (<SkewPolynomial>other)._list_c() 1135 cdef Py_ssize_t i, j 1136 cdef Py_ssize_t da = self.degree(), db = other.degree() 1137 parent = self._parent 1138 if db < 0: 1139 raise ZeroDivisionError 1140 if da < db: 1141 return self._new_c([],parent), self 1142 try: 1143 inv = self.base_ring()(~b[db]) 1144 except (ZeroDivisionError, TypeError): 1145 raise NotImplementedError("the leading coefficient of the divisor is not invertible") 1146 cdef list q = [ ] 1147 parent = self._parent 1148 for i from da-db >= i >= 0: 1149 c = parent.twist_map(i)(inv) * a[i+db] 1150 for j from 0 <= j < db: 1151 a[i+j] -= c * parent.twist_map(i)(b[j]) 1152 q.append(c) 1153 q.reverse() 1154 return self._new_c(q,parent), self._new_c(a[:db],parent,1) 1155 1156 1157 def quo_rem(self,other,side=Right): 1158 r""" 1159 INPUT: 1160 1161 - ``other`` -- a skew polynomial ring over the same 1162 base ring 1163 1164 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1165 1166 OUTPUT: 1167 1168 - the quotient and the remainder of the euclidean 1169 division of this skew polynomial by ``other`` 1170 1171 .. NOTE:: 1172 1173 Doesn't work if the leading coefficient of the divisor 1174 is not a unit. 1175 1176 EXAMPLES:: 1177 1178 sage: R.<t> = ZZ[] 1179 sage: sigma = R.hom([t+1]) 1180 sage: S.<x> = R['x',sigma] 1181 sage: a = S.random_element(degree=4); a 1182 t^2*x^4 + (-12*t^2 - 2*t - 1)*x^3 + (-95*t^2 + t + 2)*x^2 + (-t^2 + t)*x + 2*t - 8 1183 sage: b = S.random_element(monic=True); b 1184 x^2 + (4*t^2 - t - 2)*x - t^2 + t - 1 1185 sage: q,r = a.quo_rem(b) # right euclidean division 1186 sage: a == q*b + r 1187 True 1188 1189 The left euclidean division doesn't work over this `S` because 1190 Sage cannot invert `\sigma`:: 1191 1192 sage: q,r = a.quo_rem(b,side=Left) 1193 Traceback (most recent call last): 1194 ... 1195 NotImplementedError 1196 1197 In any case, the leading coefficient of the divisor need to be 1198 invertible:: 1199 1200 sage: c = S.random_element(); c 1201 (-4*t^2 + t)*x^2 - 2*t^2*x + 5*t^2 - 6*t - 4 1202 sage: a.quo_rem(c) # right euclidean division 1203 Traceback (most recent call last): 1204 ... 1205 NotImplementedError: the leading coefficient of the divisor is not invertible 1206 1207 When the base ring is a finite field, everything works fine:: 1208 1209 sage: k.<t> = GF(5^3) 1210 sage: Frob = k.frobenius_endomorphism() 1211 sage: S.<x> = k['x',Frob] 1212 sage: a = 4*x^3 + (t^2 + 2*t + 4)*x^2 + (4*t^2 + 2*t + 4)*x + t^2 + 2 1213 sage: b = (3*t + 2)*x^2 + 2*t*x + 3*t^2 1214 sage: q,r = a.quo_rem(b) # right euclidean division 1215 sage: a == q*b + r 1216 True 1217 sage: q,r = a.quo_rem(b,side=Left) 1218 sage: a == b*q + r 1219 True 1220 """ 1221 if side == Right: 1222 return self.rquo_rem(other) 1223 else: 1224 return self.lquo_rem(other) 1225 1226 1227 def rem(self,other,side=Right): 1228 r""" 1229 INPUT: 1230 1231 - ``other`` -- a skew polynomial ring over the same 1232 base ring 1233 1234 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1235 1236 OUTPUT: 1237 1238 - the remainder of the euclidean division of this 1239 skew polynomial by ``other`` 1240 1241 .. NOTE:: 1242 1243 Doesn't work if the leading coefficient of the divisor 1244 is not a unit. 1245 1246 EXAMPLES:: 1247 1248 sage: R.<t> = ZZ[] 1249 sage: sigma = R.hom([t+1]) 1250 sage: S.<x> = R['x',sigma] 1251 sage: b = x^2 + 2*t*x + 2 1252 sage: a = (x+t)*b + t*x + 1 1253 sage: a.rem(b) # right euclidean division 1254 t*x + 1 1255 1256 The left euclidean division doesn't work over this `S` because 1257 Sage cannot invert `\sigma`:: 1258 1259 sage: a.rem(b,side=Left) 1260 Traceback (most recent call last): 1261 ... 1262 NotImplementedError 1263 1264 In any case, the leading coefficient of the divisor need to be 1265 invertible:: 1266 1267 sage: (a*t).rem(b*t) 1268 Traceback (most recent call last): 1269 ... 1270 NotImplementedError: the leading coefficient of the divisor is not invertible 1271 1272 When the base ring is a finite field, everything works fine:: 1273 1274 sage: k.<t> = GF(5^3) 1275 sage: Frob = k.frobenius_endomorphism() 1276 sage: S.<x> = k['x',Frob] 1277 sage: b = x^2 + 2*t*x + 2 1278 sage: a = (x+t)*b + t*x + 1 1279 sage: a.rem(b) # right euclidean division 1280 t*x + 1 1281 sage: a.rem(b,side=Left) 1282 (3*t^2 + 2)*x + 2*t^2 1283 """ 1284 _,r = self.quo_rem(other,side=side) 1285 return r 1286 1287 1288 def __mod__(self,other): 1289 """ 1290 Return the remainder in the *right* euclidean division of 1291 this skew polynomial by ``other`` 1292 1293 TESTS:: 1294 1295 sage: R.<t> = ZZ[] 1296 sage: sigma = R.hom([t+1]) 1297 sage: S.<x> = R['x',sigma] 1298 sage: b = x^2 + 2*t*x + 2 1299 sage: a = (x+t)*b + t*x + 1 1300 sage: a % b 1301 t*x + 1 1302 1303 sage: (a*t).rem(b*t) 1304 Traceback (most recent call last): 1305 ... 1306 NotImplementedError: the leading coefficient of the divisor is not invertible 1307 """ 1308 _,r = self.rquo_rem(other) 1309 return r 1310 1311 1312 def __floordiv__(self,right): 1313 """ 1314 Return the quotient of the right euclidean division of self by right. 1315 1316 The algorithm fails if the leading coefficient of the divisor (right) 1317 is not invertible. 1318 1319 .. SEEALSO:: 1320 1321 1322 TESTS:: 1323 1324 sage: R.<t> = QQ[] 1325 sage: sigma = R.hom([t+1]) 1326 sage: S.<x> = R['x',sigma] 1327 sage: b = x^2 + t 1328 sage: a = (x^2 + t*x + 1)*b + t^3*x 1329 sage: a // b 1330 x^2 + t*x + 1 1331 1332 sage: (t*a) // (t*b) 1333 Traceback (most recent call last): 1334 ... 1335 NotImplementedError: the leading coefficient of the divisor is not invertible 1336 1337 """ 1338 q,_ = self.rquo_rem(right) 1339 return q 1340 1341 1342 cpdef RingElement _div_(self, RingElement right): 1343 """ 1344 Not Implemented (since localization of Ore rings is 1345 not yet implemented, see trac #13215). 1346 1347 Use the operator `//` even for exact division. 1348 1349 EXAMPLES:: 1350 1351 sage: R.<t> = QQ[] 1352 sage: sigma = R.hom([t+1]) 1353 sage: S.<x> = R['x',sigma] 1354 sage: a = x^5 + (t + 2)*x^2 + t^2 1355 sage: b = x^3 + 4*t 1356 sage: c = a*b 1357 1358 sage: c / b 1359 Traceback (most recent call last): 1360 ... 1361 NotImplementedError: Please use `//` (even for exact division) 1362 1363 sage: c // b == a 1364 True 1365 """ 1366 raise NotImplementedError("Please use `//` (even for exact division)") 1367 1368 1369 def is_divisible_by(self,other,side=Right): 1370 """ 1371 INPUT: 1372 1373 - ``other`` -- a skew polynomial over the same base 1374 1375 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1376 1377 OUTPUT: 1378 1379 Return True iff self is divisible by other on ``side`` 1380 1381 EXAMPLES:: 1382 1383 sage: k.<t> = GF(5^3) 1384 sage: Frob = k.frobenius_endomorphism() 1385 sage: S.<x> = k['x',Frob] 1386 sage: a = x^2 + t*x + t^2 + 3 1387 sage: b = x^3 + (t + 1)*x^2 + 1 1388 sage: c = a*b 1389 sage: c.is_divisible_by(a) # on the right 1390 False 1391 sage: c.is_divisible_by(b) 1392 True 1393 sage: c.is_divisible_by(a,side=Left) 1394 True 1395 sage: c.is_divisible_by(b,side=Left) 1396 False 1397 1398 Divisibility by 0 has no sense:: 1399 1400 sage: c.is_divisible_by(S(0)) 1401 Traceback (most recent call last): 1402 ... 1403 ZeroDivisionError 1404 1405 This function does not work if the leading coefficient of the divisor 1406 is not a unit:: 1407 1408 sage: R.<t> = QQ[] 1409 sage: sigma = R.hom([t+1]) 1410 sage: S.<x> = R['x',sigma] 1411 sage: a = x^2 + 2*x + t 1412 sage: b = (t+1)*x + t^2 1413 sage: c = a*b 1414 sage: c.is_divisible_by(b) 1415 Traceback (most recent call last): 1416 ... 1417 NotImplementedError: the leading coefficient of the divisor is not invertible 1418 """ 1419 return self.rem(other,side=side).is_zero() 1420 1421 1422 def divides(self,other,side=Right): 1423 """ 1424 INPUT: 1425 1426 - ``other`` -- a skew polynomial over the same base 1427 1428 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1429 1430 OUTPUT: 1431 1432 Return True iff self divides other on ``side`` 1433 1434 EXAMPLES:: 1435 1436 sage: k.<t> = GF(5^3) 1437 sage: Frob = k.frobenius_endomorphism() 1438 sage: S.<x> = k['x',Frob] 1439 sage: a = x^2 + t*x + t^2 + 3 1440 sage: b = x^3 + (t + 1)*x^2 + 1 1441 sage: c = a*b 1442 sage: a.divides(c) # on the right 1443 False 1444 sage: b.divides(c) 1445 True 1446 sage: a.divides(c,side=Left) 1447 True 1448 sage: b.divides(c,side=Left) 1449 False 1450 1451 Divisibility by 0 has no sense:: 1452 1453 sage: S(0).divides(c) 1454 Traceback (most recent call last): 1455 ... 1456 ZeroDivisionError 1457 1458 This function does not work if the leading coefficient of the divisor 1459 is not a unit:: 1460 1461 sage: R.<t> = QQ[] 1462 sage: sigma = R.hom([t+1]) 1463 sage: S.<x> = R['x',sigma] 1464 sage: a = x^2 + 2*x + t 1465 sage: b = (t+1)*x + t^2 1466 sage: c = a*b 1467 sage: b.divides(c) 1468 Traceback (most recent call last): 1469 ... 1470 NotImplementedError: the leading coefficient of the divisor is not invertible 1471 """ 1472 return other.rem(self,side=side).is_zero() 1473 1474 1475 # greastest commun divisor 1476 # ------------------------ 1477 1478 def lxgcd(self,other,monic=True): 1479 """ 1480 INPUT: 1481 1482 - ``other`` -- an other skew polynomial over the same 1483 base 1484 1485 - ``monic`` -- boolean (default: True) 1486 1487 OUTPUT: 1488 1489 - The left gcd of self and other, that is a skew polynomial 1490 `g` with the following property: any skew polynomial is 1491 divisible on the left by `g` iff it is divisible on the left 1492 by both ``self`` and ``other``. 1493 If monic is ``True``, `g` is in addition monic. (With this 1494 extra condition, it is uniquely determined.) 1495 1496 - Two skew polynomials `u` and `v` such that: 1497 1498 .. MATH:: 1499 1500 g = self * u + other * v 1501 1502 .. NOTE:: 1503 1504 Works only if two following conditions are fulfilled 1505 (otherwise left gcd do not exist in general): 1506 1) the base ring is a field and 2) the twist map on 1507 this field is bijective. 1508 1509 EXAMPLES:: 1510 1511 sage: k.<t> = GF(5^3) 1512 sage: Frob = k.frobenius_endomorphism() 1513 sage: S.<x> = k['x',Frob] 1514 sage: a = (x + t) * (x^2 + t*x + 1) 1515 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1516 sage: g,u,v = a.lxgcd(b); g 1517 x + t 1518 sage: a*u + b*v == g 1519 True 1520 1521 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1522 1523 sage: g,u,v = a.lxgcd(b,monic=False); g 1524 2*t*x + 4*t + 2 1525 sage: a*u + b*v == g 1526 True 1527 1528 The base ring need to be a field:: 1529 1530 sage: R.<t> = QQ[] 1531 sage: sigma = R.hom([t+1]) 1532 sage: S.<x> = R['x',sigma] 1533 sage: a = (x + t) * (x^2 + t*x + 1) 1534 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1535 sage: a.lxgcd(b) 1536 Traceback (most recent call last): 1537 ... 1538 TypeError: the base ring must be a field 1539 1540 And the twist map need to be bijective:: 1541 1542 sage: FR = R.fraction_field() 1543 sage: f = FR.hom([FR(t)^2]) 1544 sage: S.<x> = FR['x',f] 1545 sage: a = (x + t) * (x^2 + t*x + 1) 1546 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1547 sage: a.lxgcd(b) 1548 Traceback (most recent call last): 1549 ... 1550 NotImplementedError 1551 """ 1552 if not isinstance(self.base_ring(),Field): 1553 raise TypeError("the base ring must be a field") 1554 G = self 1555 U = self._parent(1) 1556 if other.is_zero(): 1557 V = self._parent(0) 1558 else: 1559 V1 = self._parent(0) 1560 V3 = other 1561 while not V3.is_zero(): 1562 Q,R = G.lquo_rem(V3) 1563 T = U - V1*Q 1564 U = V1 1565 G = V3 1566 V1 = T 1567 V3 = R 1568 V,_ = (G-self*U).lquo_rem(other) 1569 if monic: 1570 lc = ~G.leading_coefficient() 1571 lc = self._parent.twist_map(-G.degree())(lc) 1572 G = G*lc 1573 U = U*lc 1574 V = V*lc 1575 return G,U,V 1576 1577 1578 def rxgcd(self,other,monic=True): 1579 """ 1580 INPUT: 1581 1582 - ``other`` -- an other skew polynomial over the same 1583 base 1584 1585 - ``monic`` -- boolean (default: True) 1586 1587 OUTPUT: 1588 1589 - The right gcd of self and other, that is a skew polynomial 1590 `g` with the following property: any skew polynomial is 1591 divisible on the right by `g` iff it is divisible on the right 1592 by both ``self`` and ``other``. 1593 If monic is ``True``, `g` is in addition monic. (With this 1594 extra condition, it is uniquely determined.) 1595 1596 - Two skew polynomials `u` and `v` such that: 1597 1598 .. MATH:: 1599 1600 g = u * self + v * other 1601 1602 .. NOTE:: 1603 1604 Works only if the base ring is a field (otherwise right 1605 gcd do not exist in general). 1606 1607 EXAMPLES:: 1608 1609 sage: k.<t> = GF(5^3) 1610 sage: Frob = k.frobenius_endomorphism() 1611 sage: S.<x> = k['x',Frob] 1612 sage: a = (x^2 + t*x + 1) * (x + t) 1613 sage: b = 2 * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1614 sage: g,u,v = a.rxgcd(b); g 1615 x + t 1616 sage: u*a + v*b == g 1617 True 1618 1619 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1620 1621 sage: g,u,v = a.rxgcd(b,monic=False); g 1622 (4*t^2 + 4*t + 1)*x + 4*t^2 + 4*t + 3 1623 sage: u*a + v*b == g 1624 True 1625 1626 The base ring need to be a field:: 1627 1628 sage: R.<t> = QQ[] 1629 sage: sigma = R.hom([t+1]) 1630 sage: S.<x> = R['x',sigma] 1631 sage: a = (x^2 + t*x + 1) * (x + t) 1632 sage: b = 2 * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1633 sage: a.rxgcd(b) 1634 Traceback (most recent call last): 1635 ... 1636 TypeError: the base ring must be a field 1637 """ 1638 if not isinstance(self.base_ring(),Field): 1639 raise TypeError("the base ring must be a field") 1640 G = self 1641 U = self._parent(1) 1642 if other.is_zero(): 1643 V = self._parent(0) 1644 else: 1645 V1 = self._parent(0) 1646 V3 = other 1647 while not V3.is_zero(): 1648 Q, R = G.rquo_rem(V3) 1649 T = U - Q*V1 1650 U = V1 1651 G = V3 1652 V1 = T 1653 V3 = R 1654 V,_ = (G-U*self).rquo_rem(other) 1655 if monic: 1656 lc = ~G.leading_coefficient() 1657 G = lc*G 1658 U = lc*U 1659 V = lc*V 1660 return G,U,V 1661 1662 1663 def xgcd(self,other,side=Right,monic=True): 1664 """ 1665 INPUT: 1666 1667 - ``other`` -- an other skew polynomial over the same 1668 base 1669 1670 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1671 1672 - ``monic`` -- boolean (default: True) 1673 1674 OUTPUT: 1675 1676 - The left (resp. right) gcd of self and other, that is a 1677 skew polynomial `g` with the following property: any skew 1678 polynomial is divisible on the left (resp. right) by `g` 1679 iff it is divisible on the left (resp. right) by both 1680 ``self`` and ``other``. 1681 If monic is ``True``, `g` is in addition monic. (With this 1682 extra condition, it is uniquely determined.) 1683 1684 - Two skew polynomials `u` and `v` such that: 1685 1686 - (if side=Left) `g = self * u + other * v` 1687 1688 - (if side=Right) `g = u * self + v * other` 1689 1690 .. NOTE:: 1691 1692 Works only if the base ring is a field and, when 1693 ``side`` is ``Left`` if the twist map on this field 1694 is bijective (otherwise gcd do not exist in general). 1695 1696 EXAMPLES:: 1697 1698 sage: k.<t> = GF(5^3) 1699 sage: Frob = k.frobenius_endomorphism() 1700 sage: S.<x> = k['x',Frob] 1701 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1702 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1703 sage: g,u,v = a.xgcd(b); g # side = Right 1704 x + t 1705 sage: u*a + v*b == g 1706 True 1707 sage: g,u,v = a.xgcd(b,side=Left); g 1708 x + 2*t 1709 sage: a*u + b*v == g 1710 True 1711 1712 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1713 1714 sage: g,u,v = a.xgcd(b,monic=False); g 1715 (2*t^2 + 2*t + 4)*x + 2*t^2 + 3*t + 4 1716 sage: u*a + v*b == g 1717 True 1718 sage: g,u,v = a.xgcd(b,side=Left,monic=False); g 1719 (2*t^2 + 3*t + 4)*x + 2*t^2 + t + 3 1720 sage: a*u + b*v == g 1721 True 1722 1723 The base ring need to be a field:: 1724 1725 sage: R.<t> = QQ[] 1726 sage: sigma = R.hom([t+1]) 1727 sage: S.<x> = R['x',sigma] 1728 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1729 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1730 sage: a.xgcd(b) 1731 Traceback (most recent call last): 1732 ... 1733 TypeError: the base ring must be a field 1734 sage: a.xgcd(b,side=Left) 1735 Traceback (most recent call last): 1736 ... 1737 TypeError: the base ring must be a field 1738 1739 And the twist map need to be bijective:: 1740 1741 sage: FR = R.fraction_field() 1742 sage: f = FR.hom([FR(t)^2]) 1743 sage: S.<x> = FR['x',f] 1744 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1745 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1746 sage: g,u,v = a.xgcd(b); g 1747 x + t 1748 sage: g,u,v = a.xgcd(b,side=Left); g 1749 Traceback (most recent call last): 1750 ... 1751 NotImplementedError 1752 """ 1753 if side == Right: 1754 return self.rxgcd(other,monic=monic) 1755 else: 1756 return self.lxgcd(other,monic=monic) 1757 1758 1759 def rgcd(self,other,monic=True): 1760 """ 1761 INPUT: 1762 1763 - ``other`` -- an other skew polynomial over the same 1764 base 1765 1766 - ``monic`` -- boolean (default: True) 1767 1768 OUTPUT: 1769 1770 - The right gcd of self and other, that is a skew polynomial 1771 `g` with the following property: any skew polynomial is 1772 divisible on the right by `g` iff it is divisible on the right 1773 by both ``self`` and ``other``. 1774 If monic is ``True``, `g` is in addition monic. (With this 1775 extra condition, it is uniquely determined.) 1776 1777 .. NOTE:: 1778 1779 Works only if the base ring is a field (otherwise right 1780 gcd do not exist in general). 1781 1782 EXAMPLES:: 1783 1784 sage: k.<t> = GF(5^3) 1785 sage: Frob = k.frobenius_endomorphism() 1786 sage: S.<x> = k['x',Frob] 1787 sage: a = (x^2 + t*x + 1) * (x + t) 1788 sage: b = 2 * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1789 sage: a.rgcd(b) 1790 x + t 1791 1792 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1793 1794 sage: a.rgcd(b,monic=False) 1795 (4*t^2 + 4*t + 1)*x + 4*t^2 + 4*t + 3 1796 1797 The base ring need to be a field:: 1798 1799 sage: R.<t> = QQ[] 1800 sage: sigma = R.hom([t+1]) 1801 sage: S.<x> = R['x',sigma] 1802 sage: a = (x^2 + t*x + 1) * (x + t) 1803 sage: b = 2 * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1804 sage: a.rgcd(b) 1805 Traceback (most recent call last): 1806 ... 1807 TypeError: the base ring must be a field 1808 """ 1809 if not isinstance(self.base_ring(),Field): 1810 raise TypeError("the base ring must be a field") 1811 if other.is_zero(): 1812 return self 1813 A = self 1814 B = other 1815 while not B.is_zero(): 1816 A,B = B, A % B 1817 if monic: 1818 A = A.rmonic() 1819 return A 1820 1821 1822 def lgcd(self,other,monic=True): 1823 """ 1824 INPUT: 1825 1826 - ``other`` -- an other skew polynomial over the same 1827 base 1828 1829 - ``monic`` -- boolean (default: True) 1830 1831 OUTPUT: 1832 1833 - The left gcd of self and other, that is a skew polynomial 1834 `g` with the following property: any skew polynomial is 1835 divisible on the left by `g` iff it is divisible on the left 1836 by both ``self`` and ``other``. 1837 If monic is ``True``, `g` is in addition monic. (With this 1838 extra condition, it is uniquely determined.) 1839 1840 .. NOTE:: 1841 1842 Works only if two following conditions are fulfilled 1843 (otherwise left gcd do not exist in general): 1844 1) the base ring is a field and 2) the twist map on 1845 this field is bijective. 1846 1847 EXAMPLES:: 1848 1849 sage: k.<t> = GF(5^3) 1850 sage: Frob = k.frobenius_endomorphism() 1851 sage: S.<x> = k['x',Frob] 1852 sage: a = (x + t) * (x^2 + t*x + 1) 1853 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1854 sage: a.lgcd(b) 1855 x + t 1856 1857 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1858 1859 sage: a.lgcd(b,monic=False) 1860 2*t*x + 4*t + 2 1861 1862 The base ring need to be a field:: 1863 1864 sage: R.<t> = QQ[] 1865 sage: sigma = R.hom([t+1]) 1866 sage: S.<x> = R['x',sigma] 1867 sage: a = (x + t) * (x^2 + t*x + 1) 1868 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1869 sage: a.lgcd(b) 1870 Traceback (most recent call last): 1871 ... 1872 TypeError: the base ring must be a field 1873 1874 And the twist map need to be bijective:: 1875 1876 sage: FR = R.fraction_field() 1877 sage: f = FR.hom([FR(t)^2]) 1878 sage: S.<x> = FR['x',f] 1879 sage: a = (x + t) * (x^2 + t*x + 1) 1880 sage: b = 2 * (x + t) * (x^3 + (t+1)*x^2 + t^2) 1881 sage: a.lgcd(b) 1882 Traceback (most recent call last): 1883 ... 1884 NotImplementedError 1885 """ 1886 if not isinstance(self.base_ring(),Field): 1887 raise TypeError("the base ring must be a field") 1888 if other.is_zero(): 1889 return self 1890 A = self 1891 B = other 1892 while not B.is_zero(): 1893 A,B = B, A.rem(B,side=Left) 1894 if monic: 1895 A = A.lmonic() 1896 return A 1897 1898 1899 def gcd(self,other,side=Right,monic=True): 1900 """ 1901 INPUT: 1902 1903 - ``other`` -- an other skew polynomial over the same 1904 base 1905 1906 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1907 1908 - ``monic`` -- boolean (default: True) 1909 1910 OUTPUT: 1911 1912 - The left (resp. right) gcd of self and other, that is a 1913 skew polynomial `g` with the following property: any skew 1914 polynomial is divisible on the left (resp. right) by `g` 1915 iff it is divisible on the left (resp. right) by both 1916 ``self`` and ``other``. 1917 If monic is ``True``, `g` is in addition monic. (With this 1918 extra condition, it is uniquely determined.) 1919 1920 .. NOTE:: 1921 1922 Works only if the base ring is a field and, when 1923 ``side=Left`` if, in addition, the twist map on this 1924 field is bijective. 1925 1926 EXAMPLES:: 1927 1928 sage: k.<t> = GF(5^3) 1929 sage: Frob = k.frobenius_endomorphism() 1930 sage: S.<x> = k['x',Frob] 1931 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1932 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1933 sage: a.gcd(b); # side = Right 1934 x + t 1935 sage: a.gcd(b,side=Left) 1936 x + 2*t 1937 1938 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 1939 1940 sage: a.gcd(b,monic=False) 1941 (2*t^2 + 2*t + 4)*x + 2*t^2 + 3*t + 4 1942 sage: a.gcd(b,side=Left,monic=False) 1943 (2*t^2 + 3*t + 4)*x + 2*t^2 + t + 3 1944 1945 The base ring need to be a field:: 1946 1947 sage: R.<t> = QQ[] 1948 sage: sigma = R.hom([t+1]) 1949 sage: S.<x> = R['x',sigma] 1950 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1951 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1952 sage: a.gcd(b) 1953 Traceback (most recent call last): 1954 ... 1955 TypeError: the base ring must be a field 1956 sage: a.gcd(b,side=Left) 1957 Traceback (most recent call last): 1958 ... 1959 TypeError: the base ring must be a field 1960 1961 And, when ``side=Left``, the twist map need to be bijective:: 1962 1963 sage: FR = R.fraction_field() 1964 sage: f = FR.hom([FR(t)^2]) 1965 sage: S.<x> = FR['x',f] 1966 sage: a = (x + 2*t) * (x^2 + t*x + 1) * (x + t) 1967 sage: b = 2 * (x + 2*t) * (x^3 + (t+1)*x^2 + t^2) * (x + t) 1968 sage: a.gcd(b) 1969 x + t 1970 sage: a.gcd(b,side=Left) 1971 Traceback (most recent call last): 1972 ... 1973 NotImplementedError 1974 """ 1975 if side == Right: 1976 return self.rgcd(other,monic=monic) 1977 else: 1978 return self.lgcd(other,monic=monic) 1979 1980 1981 # Lowest common multiple 1982 # ---------------------- 1983 1984 def llcm(self,other,monic=True): 1985 """ 1986 INPUT: 1987 1988 - ``other`` -- an other skew polynomial over the same 1989 base 1990 1991 - ``monic`` -- boolean (default: True) 1992 1993 OUTPUT: 1994 1995 - The left lcm of self and other, that is a skew polynomial 1996 `g` with the following property: any skew polynomial divides 1997 `g` on the *right* iff it divides both ``self`` and ``other`` 1998 on the *right*. 1999 If monic is ``True``, `g` is in addition monic. (With this 2000 extra condition, it is uniquely determined.) 2001 2002 .. NOTE:: 2003 2004 Works only if the base ring is a field (otherwise left 2005 lcm do not exist in general). 2006 2007 EXAMPLES:: 2008 2009 sage: k.<t> = GF(5^3) 2010 sage: Frob = k.frobenius_endomorphism() 2011 sage: S.<x> = k['x',Frob] 2012 sage: a = (x + t^2) * (x + t) 2013 sage: b = 2 * (x^2 + t + 1) * (x * t) 2014 sage: c = a.llcm(b); c 2015 x^5 + (2*t^2 + t + 4)*x^4 + (3*t^2 + 4)*x^3 + (3*t^2 + 3*t + 2)*x^2 + (t^2 + t + 2)*x 2016 sage: c.is_divisible_by(a,side=Right) 2017 True 2018 sage: c.is_divisible_by(b,side=Right) 2019 True 2020 sage: a.degree() + b.degree() == c.degree() + a.rgcd(b).degree() 2021 True 2022 2023 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 2024 2025 sage: a.llcm(b,monic=False) 2026 (t^2 + t)*x^5 + (4*t^2 + 4*t + 1)*x^4 + (t + 1)*x^3 + (t^2 + 2)*x^2 + (3*t + 4)*x 2027 2028 The base ring need to be a field:: 2029 2030 sage: R.<t> = QQ[] 2031 sage: sigma = R.hom([t+1]) 2032 sage: S.<x> = R['x',sigma] 2033 sage: a = (x + t^2) * (x + t) 2034 sage: b = 2 * (x^2 + t + 1) * (x * t) 2035 sage: a.llcm(b) 2036 Traceback (most recent call last): 2037 ... 2038 TypeError: the base ring must be a field 2039 """ 2040 if not isinstance(self.base_ring(),Field): 2041 raise TypeError("the base ring must be a field") 2042 if self.is_zero() or other.is_zero(): 2043 raise ZeroDivisionError 2044 U = self._parent(1) 2045 G = self 2046 V1 = self._parent(0) 2047 V3 = other 2048 while not V3.is_zero(): 2049 Q, R = G.rquo_rem(V3) 2050 T = U - Q*V1 2051 U = V1 2052 G = V3 2053 V1 = T 2054 V3 = R 2055 V1 = V1*self 2056 if monic: 2057 V1 = V1.rmonic() 2058 return V1 2059 2060 2061 def rlcm(self,other,monic=True): 2062 """ 2063 INPUT: 2064 2065 - ``other`` -- an other skew polynomial over the same 2066 base 2067 2068 - ``monic`` -- boolean (default: True) 2069 2070 OUTPUT: 2071 2072 - The right lcm of self and other, that is a skew polynomial 2073 `g` with the following property: any skew polynomial divides 2074 `g` on the *left* iff it divides both ``self`` and ``other`` 2075 on the *left*. 2076 If monic is ``True``, `g` is in addition monic. (With this 2077 extra condition, it is uniquely determined.) 2078 2079 .. NOTE:: 2080 2081 Works only if two following conditions are fulfilled 2082 (otherwise right lcm do not exist in general): 2083 1) the base ring is a field and 2) the twist map on 2084 this field is bijective. 2085 2086 EXAMPLES:: 2087 2088 sage: k.<t> = GF(5^3) 2089 sage: Frob = k.frobenius_endomorphism() 2090 sage: S.<x> = k['x',Frob] 2091 sage: a = (x + t) * (x + t^2) 2092 sage: b = 2 * (x + t) * (x^2 + t + 1) 2093 sage: c = a.rlcm(b); c 2094 x^4 + (2*t^2 + t + 2)*x^3 + (3*t^2 + 4*t + 1)*x^2 + (3*t^2 + 4*t + 1)*x + t^2 + 4 2095 sage: c.is_divisible_by(a,side=Left) 2096 True 2097 sage: c.is_divisible_by(b,side=Left) 2098 True 2099 sage: a.degree() + b.degree() == c.degree() + a.lgcd(b).degree() 2100 True 2101 2102 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 2103 2104 sage: a.rlcm(b,monic=False) 2105 2*t*x^4 + (3*t + 1)*x^3 + (4*t^2 + 4*t + 3)*x^2 + (3*t^2 + 4*t + 2)*x + 3*t^2 + 2*t + 3 2106 2107 The base ring need to be a field:: 2108 2109 sage: R.<t> = QQ[] 2110 sage: sigma = R.hom([t+1]) 2111 sage: S.<x> = R['x',sigma] 2112 sage: a = (x + t) * (x + t^2) 2113 sage: b = 2 * (x + t) * (x^2 + t + 1) 2114 sage: a.rlcm(b) 2115 Traceback (most recent call last): 2116 ... 2117 TypeError: the base ring must be a field 2118 2119 And the twist map need to be bijective:: 2120 2121 sage: FR = R.fraction_field() 2122 sage: f = FR.hom([FR(t)^2]) 2123 sage: S.<x> = FR['x',f] 2124 sage: a = (x + t) * (x + t^2) 2125 sage: b = 2 * (x + t) * (x^2 + t + 1) 2126 sage: a.rlcm(b) 2127 Traceback (most recent call last): 2128 ... 2129 NotImplementedError 2130 """ 2131 if not isinstance(self.base_ring(),Field): 2132 raise TypeError("the base ring must be a field") 2133 if self.is_zero() or other.is_zero(): 2134 raise ZeroDivisionError 2135 R = self.parent() 2136 U = R.one_element() 2137 G = self 2138 V1 = R.zero_element() 2139 V3 = other 2140 while not V3.is_zero(): 2141 Q, R = G.lquo_rem(V3) 2142 T = U - V1*Q 2143 U = V1 2144 G = V3 2145 V1 = T 2146 V3 = R 2147 V1 = self*V1 2148 if monic: 2149 V1 = V1.lmonic() 2150 return V1 2151 2152 2153 def lcm(self,other,side=Left,monic=True): 2154 """ 2155 INPUT: 2156 2157 - ``other`` -- an other skew polynomial over the same 2158 base 2159 2160 - ``side`` -- ``Left`` or ``Right`` (default: Left) 2161 2162 - ``monic`` -- boolean (default: True) 2163 2164 OUTPUT: 2165 2166 - The left (resp. right) lcm of self and other, that is a 2167 skew polynomial `g` with the following property: any skew 2168 polynomial divides `g` on the right (resp. left) iff it 2169 divides both ``self`` and ``other`` on the rignt (resp. 2170 left) 2171 If monic is ``True``, `g` is in addition monic. (With this 2172 extra condition, it is uniquely determined.) 2173 2174 .. NOTE:: 2175 2176 Works only if the base ring is a field and, when 2177 ``side=Right`` if, in addition, the twist map on this 2178 field is bijective. 2179 2180 EXAMPLES:: 2181 2182 sage: k.<t> = GF(5^3) 2183 sage: Frob = k.frobenius_endomorphism() 2184 sage: S.<x> = k['x',Frob] 2185 sage: a = (x + 2*t) * (x + t^2) * (x + t) 2186 sage: b = 2 * (x + 2*t) * (x + t + 1) * (x + t) 2187 2188 sage: c = a.lcm(b); c # side = Left 2189 x^5 + (2*t^2 + 4*t)*x^4 + (2*t^2 + 3*t + 4)*x^3 + (t^2 + t + 3)*x^2 + (2*t^2 + 4*t + 3)*x + t^2 + t + 1 2190 sage: c.is_divisible_by(a,side=Right) 2191 True 2192 sage: c.is_divisible_by(b,side=Right) 2193 True 2194 2195 sage: c = a.lcm(b,side=Right); c 2196 x^5 + (3*t^2 + 2*t + 1)*x^4 + (4*t^2 + t + 2)*x^3 + (4*t^2 + 2*t + 3)*x^2 + 4*t^2*x + t^2 + 2 2197 sage: c.is_divisible_by(a,side=Left) 2198 True 2199 sage: c.is_divisible_by(b,side=Left) 2200 True 2201 2202 Specifying ``monic=False``, we *can* get a nonmonic gcd:: 2203 2204 sage: a.lcm(b,monic=False) 2205 (3*t^2 + 3*t + 4)*x^5 + (2*t^2 + 4*t + 1)*x^4 + (t^2 + t + 1)*x^3 + (2*t^2 + t + 4)*x^2 + (t^2 + 3*t + 3)*x + t^2 + 1 2206 sage: a.lcm(b,side=Right,monic=False) 2207 (4*t + 4)*x^5 + (3*t^2 + t + 1)*x^4 + (t^2 + 4)*x^3 + (4*t^2 + 2*t + 4)*x^2 + (3*t^2 + t)*x + 2*t^2 + 2 2208 2209 The base ring need to be a field:: 2210 2211 sage: R.<t> = QQ[] 2212 sage: sigma = R.hom([t+1]) 2213 sage: S.<x> = R['x',sigma] 2214 sage: a = (x + 2*t) * (x + t^2) * (x + t) 2215 sage: b = 2 * (x + 2*t) * (x + t + 1) * (x + t) 2216 sage: a.lcm(b) 2217 Traceback (most recent call last): 2218 ... 2219 TypeError: the base ring must be a field 2220 sage: a.lcm(b,side=Right) 2221 Traceback (most recent call last): 2222 ... 2223 TypeError: the base ring must be a field 2224 2225 And the twist map need to be bijective:: 2226 2227 sage: FR = R.fraction_field() 2228 sage: f = FR.hom([FR(t)^2]) 2229 sage: S.<x> = FR['x',f] 2230 sage: a = (x + 2*t) * (x + t^2) * (x + t) 2231 sage: b = 2 * (x + 2*t) * (x + t + 1) * (x + t) 2232 sage: a.lcm(b,side=Right) 2233 Traceback (most recent call last): 2234 ... 2235 NotImplementedError 2236 """ 2237 if side == Right: 2238 return self.rlcm(other,monic=monic) 2239 else: 2240 return self.llcm(other,monic=monic) 2241 2242 2243 # Printing 2244 # -------- 2245 2246 def _repr_(self): 2247 """ 2248 Return string representation of this skew polynomial. 2249 2250 EXAMPLES:: 2251 2252 sage: R.<t> = QQ[] 2253 sage: sigma = R.hom([t+1]) 2254 sage: S.<x> = R['x',sigma] 2255 sage: a = t^2 + 1/2*x*t; 2256 sage: a._repr_() 2257 '(1/2*t + 1/2)*x + t^2' 2258 """ 2259 return self._repr() 2260 2261 2262 def _repr(self,name=None): 2263 """ 2264 Return string representation of this skew polynomial. 2265 2266 INPUT: 2267 2268 - ``name`` -- the name of the variable (default: the 2269 name given when the skew polynomial ring was created) 2270 2271 EXAMPLES:: 2272 2273 sage: R.<t> = QQ[] 2274 sage: sigma = R.hom([t+1]) 2275 sage: S.<x> = R['x',sigma] 2276 sage: a = t^2 + 1/2*x*t; 2277 sage: a._repr() 2278 '(1/2*t + 1/2)*x + t^2' 2279 sage: a._repr(name='y') 2280 '(1/2*t + 1/2)*y + t^2' 2281 """ 2282 s = " " 2283 m = self.degree() + 1 2284 if name is None: 2285 name = self.parent().variable_name() 2286 atomic_repr = self.parent().base_ring()._repr_option('element_is_atomic') 2287 coeffs = self.list() 2288 for n in reversed(xrange(m)): 2289 x = coeffs[n] 2290 if x: 2291 if n != m-1: 2292 s += " + " 2293 x = y = repr(x) 2294 if y.find("-") == 0: 2295 y = y[1:] 2296 if not atomic_repr and n > 0 and (y.find("+") != -1 or y.find("-") != -1): 2297 x = "(%s)"%x 2298 if n > 1: 2299 var = "*%s^%s"%(name,n) 2300 elif n==1: 2301 var = "*%s"%name 2302 else: 2303 var = "" 2304 s += "%s%s"%(x,var) 2305 s = s.replace(" + -", " - ") 2306 s = re.sub(r' 1(\.0+)?\*',' ', s) 2307 s = re.sub(r' -1(\.0+)?\*',' -', s) 2308 if s == " ": 2309 return "0" 2310 return s[1:] 2311 2312 2313 def _latex_(self,name=None): 2314 """ 2315 Return a latex representation of this skew polynomial. 2316 2317 INPUT: 2318 2319 - ``name`` -- the name of the variable (default: the 2320 name given when the skew polynomial ring was created) 2321 2322 EXAMPLES:: 2323 2324 sage: R.<t> = QQ[] 2325 sage: sigma = R.hom([t+1]) 2326 sage: S.<x> = R['x',sigma] 2327 sage: a = t^2 + 1/2*x*t; 2328 sage: a._latex_() 2329 '\\left(\\frac{1}{2} t + \\frac{1}{2}\\right) x + t^{2}' 2330 sage: a._latex_(name='y') 2331 '\\left(\\frac{1}{2} t + \\frac{1}{2}\\right) y + t^{2}' 2332 """ 2333 s = " " 2334 coeffs = self.list() 2335 m = len(coeffs) 2336 if name is None: 2337 name = self.parent().latex_variable_names()[0] 2338 atomic_repr = self.parent().base_ring()._repr_option('element_is_atomic') 2339 for n in reversed(xrange(m)): 2340 x = coeffs[n] 2341 x = y = latex(x) 2342 if x != '0': 2343 if n != m-1: 2344 s += " + " 2345 if y.find("-") == 0: 2346 y = y[1:] 2347 if not atomic_repr and n > 0 and (y.find("+") != -1 or y.find("-") != -1): 2348 x = "\\left(%s\\right)"%x 2349 if n > 1: 2350 var = "|%s^{%s}"%(name,n) 2351 elif n==1: 2352 var = "|%s"%name 2353 else: 2354 var = "" 2355 s += "%s %s"%(x,var) 2356 s = s.replace(" + -", " - ") 2357 s = re.sub(" 1(\.0+)? \|"," ", s) 2358 s = re.sub(" -1(\.0+)? \|", " -", s) 2359 s = s.replace("|","") 2360 if s==" ": 2361 return "0" 2362 return s[1:].lstrip().rstrip() 2363 2364 2365 # Misc 2366 # ---- 2367 2368 def _is_atomic(self): 2369 return (self.degree() == self.valuation() and 2370 self.leading_coefficient()._is_atomic()) 2371 2372 2373 def base_ring(self): 2374 """ 2375 Return the base ring of this skew polynomial. 2376 2377 EXAMPLES:: 2378 2379 sage: R.<t> = ZZ[] 2380 sage: sigma = R.hom([t+1]) 2381 sage: S.<x> = R['x',sigma] 2382 sage: a = S.random_element() 2383 sage: a.base_ring() 2384 Univariate Polynomial Ring in t over Integer Ring 2385 sage: a.base_ring() is R 2386 True 2387 """ 2388 return self.parent().base_ring() 2389 2390 2391 def shift(self, n): 2392 """ 2393 Return this skew polynomial multiplied on the right by the 2394 power `x^n`. If `n` is negative, terms below `x^n` will be 2395 discarded. 2396 2397 EXAMPLES:: 2398 2399 sage: R.<t> = QQ[] 2400 sage: sigma = R.hom([t+1]) 2401 sage: S.<x> = R['x',sigma] 2402 sage: a = x^5 + t^4*x^4 + t^2*x^2 + t^10 2403 sage: a.shift(0) 2404 x^5 + t^4*x^4 + t^2*x^2 + t^10 2405 sage: a.shift(-1) 2406 x^4 + t^4*x^3 + t^2*x 2407 sage: a.shift(-5) 2408 1 2409 sage: a.shift(2) 2410 x^7 + t^4*x^6 + t^2*x^4 + t^10*x^2 2411 2412 One can also use the infix shift operator:: 2413 2414 sage: a >> 2 2415 x^3 + t^4*x^2 + t^2 2416 sage: a << 2 2417 x^7 + t^4*x^6 + t^2*x^4 + t^10*x^2 2418 """ 2419 if n == 0 or self.degree() < 0: 2420 return self 2421 if n > 0: 2422 return self._parent(n*[self.base_ring().zero_element()] + self.list(), check=False) 2423 if n < 0: 2424 if n > self.degree(): 2425 return self._parent([]) 2426 else: 2427 return self._parent(self.list()[-n:], check=False) 2428 2429 2430 def __lshift__(self, k): 2431 return self.shift(k) 2432 2433 2434 def __rshift__(self, k): 2435 return self.shift(-k) 2436 2437 2438 def change_variable_name(self, var): 2439 """ 2440 Return a new polynomial over the same base ring but in a different 2441 variable. 2442 2443 INPUT: 2444 2445 - ``var`` -- the name of the new variable 2446 2447 EXAMPLES:: 2448 2449 sage: R.<t> = ZZ[] 2450 sage: sigma = R.hom([t+1]) 2451 sage: S.<x> = R['x',sigma] 2452 sage: a = x^3 + (2*t + 1)*x + t^2 + 3*t + 5 2453 sage: b = a.change_variable_name('y'); b 2454 y^3 + (2*t + 1)*y + t^2 + 3*t + 5 2455 2456 Remark that a new parent is created at the same time:: 2457 2458 sage: b.parent() 2459 Skew Polynomial Ring in y over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 2460 """ 2461 parent = self._parent 2462 R = parent.base_ring()[var,parent.twist_map()] 2463 return R(self.list()) 2464 2465 2466 def __copy__(self): 2467 """ 2468 Return a copy of this skew polynomial. 2469 """ 2470 return self 2471 2472 2473 def dict(self): 2474 """ 2475 Return a sparse dictionary representation of this skew 2476 polynomial. 2477 2478 EXAMPLES:: 2479 2480 sage: R.<t> = ZZ[] 2481 sage: sigma = R.hom([t+1]) 2482 sage: S.<x> = R['x',sigma] 2483 sage: a = x^2012 + t*x^1006 + t^3 + 2*t 2484 sage: a.dict() 2485 {0: t^3 + 2*t, 2012: 1, 1006: t} 2486 """ 2487 X = {} 2488 Y = self.list() 2489 for i in xrange(len(Y)): 2490 c = Y[i] 2491 if c: 2492 X[i] = c 2493 return X 2494 2495 2496 def is_monomial(self): 2497 """ 2498 Returns True if self is a monomial, i.e., a power of the generator. 2499 2500 EXAMPLES:: 2501 2502 sage: R.<t> = ZZ[] 2503 sage: sigma = R.hom([t+1]) 2504 sage: S.<x> = R['x',sigma] 2505 sage: x.is_monomial() 2506 True 2507 sage: (x+1).is_monomial() 2508 False 2509 sage: (x^2).is_monomial() 2510 True 2511 sage: S(1).is_monomial() 2512 True 2513 2514 The coefficient must be 1:: 2515 2516 sage: (2*x^5).is_monomial() 2517 False 2518 sage: S(t).is_monomial() 2519 False 2520 2521 To allow a non-1 leading coefficient, use is_term():: 2522 2523 sage: (2*x^5).is_term() 2524 True 2525 sage: S(t).is_term() 2526 True 2527 """ 2528 return len(self.exponents()) == 1 and self.leading_coefficient() == 1 2529 2530 2531 def is_term(self): 2532 """ 2533 Return True if self is an element of the base ring times a 2534 power of the generator. 2535 2536 EXAMPLES:: 2537 2538 sage: R.<t> = ZZ[] 2539 sage: sigma = R.hom([t+1]) 2540 sage: S.<x> = R['x',sigma] 2541 sage: x.is_term() 2542 True 2543 sage: R(1).is_term() 2544 True 2545 sage: (3*x^5).is_term() 2546 True 2547 sage: (1+3*x^5).is_term() 2548 False 2549 2550 To require that the coefficient is 1, use is_monomial() instead:: 2551 2552 sage: (3*x^5).is_monomial() 2553 False 2554 """ 2555 return len(self.exponents()) == 1 2556 2557 2558 def is_gen(self): 2559 return self._is_gen 2560 2561 2562 def coefficients(self): 2563 """ 2564 Return the nonzero coefficients of the monomials 2565 appearing in self. 2566 2567 EXAMPLES:: 2568 2569 sage: R.<t> = QQ[] 2570 sage: sigma = R.hom([t+1]) 2571 sage: S.<x> = R['x',sigma] 2572 sage: a = 1 + x^4 + (t+1)*x^2 + t^2 2573 sage: a.coefficients() 2574 [t^2 + 1, t + 1, 1] 2575 """ 2576 return [c for c in self.list() if not c.is_zero()] 2577 2578 2579 def exponents(self): 2580 """ 2581 Return the exponents of the monomials appearing in self. 2582 2583 EXAMPLES:: 2584 2585 sage: R.<t> = QQ[] 2586 sage: sigma = R.hom([t+1]) 2587 sage: S.<x> = R['x',sigma] 2588 sage: a = 1 + x^4 + (t+1)*x^2 + t^2 2589 sage: a.exponents() 2590 [0, 2, 4] 2591 """ 2592 l = self.list() 2593 return [i for i in range(len(l)) if not l[i].is_zero()] 2594 2595 2596 def prec(self): 2597 """ 2598 Return the precision of this polynomial. This is always infinity, 2599 since polynomials are of infinite precision by definition (there is 2600 no big-oh). 2601 2602 EXAMPLES:: 2603 2604 sage: R.<t> = QQ[] 2605 sage: sigma = R.hom([t+1]) 2606 sage: S.<x> = R['x',sigma] 2607 sage: x.prec() 2608 +Infinity 2609 """ 2610 return infinity.infinity 2611 2612 2613 def padded_list(self,n=None): 2614 """ 2615 Return list of coefficients of self up to (but not including) 2616 `q^n`. 2617 2618 Includes 0's in the list on the right so that the list has length 2619 `n`. 2620 2621 INPUT: 2622 2623 2624 - ``n`` - (default: None); if given, an integer that 2625 is at least 0 2626 2627 2628 EXAMPLES:: 2629 2630 sage: R.<t> = QQ[] 2631 sage: sigma = R.hom([t+1]) 2632 sage: S.<x> = R['x',sigma] 2633 sage: a = 1 + t*x^3 + t^2*x^5 2634 sage: a.padded_list() 2635 [1, 0, 0, t, 0, t^2] 2636 sage: a.padded_list(10) 2637 [1, 0, 0, t, 0, t^2, 0, 0, 0, 0] 2638 sage: len(a.padded_list(10)) 2639 10 2640 sage: a.padded_list(3) 2641 [1, 0, 0] 2642 sage: a.padded_list(0) 2643 [] 2644 sage: a.padded_list(-1) 2645 Traceback (most recent call last): 2646 ... 2647 ValueError: n must be at least 0 2648 """ 2649 v = self.list() 2650 if n is None: 2651 return v 2652 if n < 0: 2653 raise ValueError("n must be at least 0") 2654 if len(v) < n: 2655 z = self._parent.base_ring().zero_element() 2656 return v + [z]*(n - len(v)) 2657 else: 2658 return v[:int(n)] 2659 2660 2661 def coeffs(self): 2662 r""" 2663 Returns ``self.list()``. 2664 2665 (It is potentially slightly faster to use 2666 ``self.list()`` directly.) 2667 2668 EXAMPLES:: 2669 2670 sage: R.<t> = QQ[] 2671 sage: sigma = R.hom([t+1]) 2672 sage: S.<x> = R['x',sigma] 2673 sage: a = 1 + x^4 + (t+1)*x^2 + t^2 2674 sage: a.coeffs() 2675 [t^2 + 1, 0, t + 1, 0, 1] 2676 """ 2677 return self.list() 2678 2679 2680 def variable_name(self): 2681 """ 2682 Return name of variable used in this polynomial as a string. 2683 2684 EXAMPLES:: 2685 2686 sage: R.<t> = QQ[] 2687 sage: sigma = R.hom([t+1]) 2688 sage: S.<x> = R['x',sigma] 2689 sage: a = x + t 2690 sage: a.variable_name() 2691 'x' 2692 """ 2693 return self.parent().variable_name() 2694 2695 2696 cdef class SkewPolynomial_generic_dense(SkewPolynomial): 2697 def __init__(self, parent, x=None, int check=1, is_gen=False, int construct=0, **kwds): 2698 """ 2699 TESTS:: 2700 2701 sage: R.<t> = QQ[] 2702 sage: sigma = R.hom([t+1]) 2703 sage: S.<x> = R['x',sigma] 2704 2705 We create a skew polynomial from a list:: 2706 2707 sage: S([t,1]) 2708 x + t 2709 2710 from another skew polynomial:: 2711 2712 sage: S(x^2 + t) 2713 x^2 + t 2714 2715 from a constant:: 2716 2717 sage: x = S(t^2 + 1); x 2718 t^2 + 1 2719 sage: x.parent() is S 2720 True 2721 """ 2722 SkewPolynomial.__init__(self, parent, is_gen=is_gen) 2723 if x is None: 2724 self.__coeffs = [] 2725 return 2726 2727 R = parent.base_ring() 2728 if PY_TYPE_CHECK(x, list): 2729 if check: 2730 self.__coeffs = [R(t) for t in x] 2731 self.__normalize() 2732 else: 2733 self.__coeffs = x 2734 return 2735 2736 if PY_TYPE_CHECK(x, SkewPolynomial): 2737 if (<Element>x)._parent is self._parent: 2738 x = list(x.list()) 2739 elif R.has_coerce_map_from((<Element>x)._parent):# is R or (<Element>x)._parent == R: 2740 try: 2741 if x.is_zero(): 2742 self.__coeffs = [] 2743 return 2744 except (AttributeError, TypeError): 2745 pass 2746 x = [x] 2747 else: 2748 self.__coeffs = [R(a, **kwds) for a in x.list()] 2749 if check: 2750 self.__normalize() 2751 return 2752 2753 elif PY_TYPE_CHECK(x, int) and x == 0: 2754 self.__coeffs = [] 2755 return 2756 2757 elif isinstance(x, dict): 2758 x = self._dict_to_list(x, R.zero_element()) 2759 2760 elif not isinstance(x, list): 2761 # We trust that the element constructors do not send x=0 2762 # if x: 2763 x = [x] # constant polynomials 2764 # else: 2765 # x = [] # zero polynomial 2766 if check: 2767 self.__coeffs = [R(z, **kwds) for z in x] 2768 self.__normalize() 2769 else: 2770 self.__coeffs = x 2771 2772 2773 cdef list _list_c(self): 2774 """ 2775 Return the list of the underlying elements of self. 2776 2777 .. WARNING:: 2778 2779 It is a priori not a copy; do not modify this list! 2780 """ 2781 return self.__coeffs 2782 2783 2784 cdef SkewPolynomial _new_c(self, list coeffs, Parent P, char check=0): 2785 """ 2786 Fast creation of a new skew polynomial 2787 """ 2788 cdef SkewPolynomial_generic_dense f = <SkewPolynomial_generic_dense>PY_NEW_SAME_TYPE(self) 2789 f._parent = P 2790 f.__coeffs = coeffs 2791 if check: 2792 f.__normalize() 2793 return f 2794 2795 2796 cdef void __normalize(self): 2797 x = self.__coeffs 2798 cdef Py_ssize_t n = len(x) - 1 2799 while n >= 0 and not x[n]: 2800 del x[n] 2801 n -= 1 2802 2803 2804 # Basic operations in place 2805 # ------------------------- 2806 2807 cdef void _inplace_add(self, SkewPolynomial_generic_dense right): 2808 """ 2809 Replace self by self+right (only for internal use). 2810 """ 2811 cdef Py_ssize_t i, min 2812 x = (<SkewPolynomial_generic_dense>self).__coeffs 2813 y = (<SkewPolynomial_generic_dense>right).__coeffs 2814 if len(x) > len(y): 2815 for i from 0 <= i < len(y): 2816 x[i] += y[i] 2817 else: 2818 x += y[len(x):] 2819 for i from 0 <= i < len(x): 2820 x[i] += y[i] 2821 x += y[len(x):] 2822 if len(x) == len(y): 2823 self.__normalize() 2824 2825 cdef void _inplace_sub(self, SkewPolynomial_generic_dense right): 2826 """ 2827 Replace self by self-right (only for internal use). 2828 """ 2829 cdef Py_ssize_t i, min 2830 cdef list x = (<SkewPolynomial_generic_dense>self).__coeffs 2831 cdef list y = (<SkewPolynomial_generic_dense>right).__coeffs 2832 if len(x) >= len(y): 2833 for i from 0 <= i < len(y): 2834 x[i] -= y[i] 2835 else: 2836 for i from 0 <= i < len(x): 2837 x[i] -= y[i] 2838 x += [-c for c in y[len(x):]] 2839 if len(x) == len(y): 2840 self.__normalize() 2841 2842 cdef void _inplace_rmul(self, SkewPolynomial_generic_dense right): 2843 """ 2844 Replace self by self*right (only for internal use). 2845 """ 2846 cdef list x = (<SkewPolynomial_generic_dense>self).__coeffs 2847 cdef list y = (<SkewPolynomial_generic_dense>right).__coeffs 2848 cdef Py_ssize_t i, k, start, end 2849 cdef Py_ssize_t d1 = len(x)-1, d2 = len(y)-1 2850 parent = self._parent 2851 if d2 == -1: 2852 (<SkewPolynomial_generic_dense>self).__coeffs = [ ] 2853 elif d1 >= 0: 2854 for k from d1 < k <= d1+d2: 2855 start = 0 if k <= d2 else k-d2 # max(0, k-d2) 2856 sum = x[start] * parent.twist_map(start)(y[k-start]) 2857 for i from start < i <= d1: 2858 sum += x[i] * parent.twist_map(i)(y[k-i]) 2859 x.append(sum) 2860 for k from d1 >= k >= 0: 2861 start = 0 if k <= d2 else k-d2 # max(0, k-d2) 2862 end = k if k <= d1 else d1 # min(k, d1) 2863 sum = x[start] * parent.twist_map(start)(y[k-start]) 2864 for i from start < i <= end: 2865 sum += x[i] * parent.twist_map(i)(y[k-i]) 2866 x[k] = sum 2867 2868 cdef void _inplace_lmul(self, SkewPolynomial_generic_dense left): 2869 """ 2870 Replace self by left*self (only for internal use). 2871 """ 2872 cdef list x = (<SkewPolynomial_generic_dense>self).__coeffs 2873 cdef list y = (<SkewPolynomial_generic_dense>left).__coeffs 2874 cdef Py_ssize_t i, k, start, end 2875 cdef Py_ssize_t d1 = len(x)-1, d2 = len(y)-1 2876 parent = self._parent 2877 if d2 == -1: 2878 (<SkewPolynomial_generic_dense>self).__coeffs = [ ] 2879 elif d1 >= 0: 2880 for k from d1 < k <= d1+d2: 2881 start = 0 if k <= d2 else k-d2 # max(0, k-d2) 2882 sum = parent.twist_map(k-start)(x[start]) * y[k-start] 2883 for i from start < i <= d1: 2884 sum += parent.twist_map(k-i)(x[i]) * y[k-i] 2885 x.append(sum) 2886 for k from d1 >= k >= 0: 2887 start = 0 if k <= d2 else k-d2 # max(0, k-d2) 2888 end = k if k <= d1 else d1 # min(k, d1) 2889 sum = parent.twist_map(k-start)(x[start]) * y[k-start] 2890 for i from start < i <= end: 2891 sum += parent.twist_map(k-i)(x[i]) * y[k-i] 2892 x[k] = sum 2893 2894 2895 # Fast exponentiation 2896 # ------------------- 2897 2898 cdef void _inplace_pow(self, Py_ssize_t n): 2899 """ 2900 Replace self by self**n. 2901 """ 2902 while n & 1 == 0: 2903 self._inplace_rmul(self) 2904 n = n >> 1 2905 cdef SkewPolynomial_generic_dense selfpow = <SkewPolynomial_generic_dense>self._new_c(list(self.__coeffs),self._parent) 2906 n = n >> 1 2907 while n != 0: 2908 selfpow._inplace_rmul(selfpow) 2909 if n&1 == 1: 2910 self._inplace_rmul(selfpow) 2911 n = n >> 1 2912 2913 2914 cpdef _pow_(self,exp,modulus=None,side=Right): 2915 """ 2916 INPUT: 2917 2918 - ``exp`` -- an Integer 2919 2920 - ``modulus`` -- a skew polynomial over the same ring (default: None) 2921 2922 - ``side`` -- ``Left`` or ``Right`` (default: Right) 2923 2924 OUTPUT: 2925 2926 If ``modulus`` is None, return ``self**exp``. 2927 2928 Otherwise, return the remainder of self**exp in the ``side`` 2929 euclidean division by ``modulus``. 2930 2931 REMARK: 2932 2933 The quotient of the underlying skew polynomial ring by the 2934 principal ideal generated by ``modulus`` is in general *not* 2935 a ring. 2936 2937 As a consequence, Sage first computes exactly ``self**exp`` 2938 and then reduce it modulo ``modulus``. 2939 2940 However, if the base ring is a finite field, Sage uses the 2941 following optimized algorithm: 2942 2943 #. One first compute a central skew polynomial `N` which is 2944 divisible by ``modulus``. (Since `N` lies in center, the 2945 quotient `K[X,\sigma]/N` inherits a ring structure.) 2946 2947 #. One compute ``self**exp`` in the quotient ring `K[X,\sigma]/N` 2948 2949 #. One reduce modulo ``modulus`` the result computed in the 2950 previous step 2951 2952 EXAMPLES:: 2953 2954 sage: k.<t> = GF(5^3) 2955 sage: Frob = k.frobenius_endomorphism() 2956 sage: S.<x> = k['x',Frob] 2957 sage: a = x + t 2958 sage: b = a^10 # short form for ``a._pow_(10)`` 2959 sage: b == a*a*a*a*a*a*a*a*a*a 2960 True 2961 2962 sage: modulus = x^3 + t*x^2 + (t+3)*x - 2 2963 sage: br = a._pow_(10,modulus); br 2964 (t^2 + t)*x^2 + (3*t^2 + 1)*x + t^2 + t 2965 sage: br == b.rem(modulus) 2966 True 2967 sage: bl = a._pow_(10,modulus,side=Left); bl 2968 (4*t^2 + 2*t + 3)*x^2 + (3*t^2 + 1)*x + 2*t + 3 2969 sage: bl == b.rem(modulus,side=Left) 2970 True 2971 2972 sage: a._pow_(10^100,modulus) # quite fast 2973 (3*t^2 + 3)*x^2 + (t^2 + 2*t + 4)*x + 4*t^2 + 2*t + 1 2974 """ 2975 cdef SkewPolynomial_generic_dense r 2976 2977 if not PY_TYPE_CHECK_EXACT(exp, Integer) or \ 2978 PY_TYPE_CHECK_EXACT(exp, int): 2979 try: 2980 exp = Integer(exp) 2981 except TypeError: 2982 raise TypeError("non-integral exponents not supported") 2983 2984 if self.degree() <= 0: 2985 return self.parent()(self[0]**exp) 2986 if exp == 0: 2987 return self.parent()(1) 2988 if exp < 0: 2989 return (~self).pow(-exp,modulus,side=side) 2990 2991 if self == self.parent().gen(): # special case x**n should be faster! 2992 P = self.parent() 2993 R = P.base_ring() 2994 v = [R.zero_element()]*exp + [R.one_element()] 2995 r = <SkewPolynomial_generic_dense>self._parent(v) 2996 else: 2997 r = <SkewPolynomial_generic_dense>self._new_c(list(self.__coeffs),self._parent) 2998 sig_on() 2999 r._inplace_pow(exp) 3000 sig_off() 3001 3002 if modulus: 3003 sig_on() 3004 if side == Right: 3005 r._inplace_rrem(modulus) 3006 else: 3007 r._inplace_lrem(modulus) 3008 sig_off() 3009 3010 return r 3011 3012 3013 def __pow__(self,exp,modulus): 3014 """ 3015 INPUT: 3016 3017 - ``exp`` -- an Integer 3018 3019 - ``modulus`` -- a skew polynomial over the same ring 3020 (default: None) 3021 3022 OUTPUT: 3023 3024 If ``modulus`` is None, return ``self**exp``. 3025 3026 Otherwise, return the remainder of self**exp in the right 3027 euclidean division by ``modulus``. 3028 3029 .. SEEALSO:: 3030 3031 :meth:`~sage.rings.polynomial.skew_polynomial_element._pow_` 3032 3033 EXAMPLES:: 3034 3035 sage: k.<t> = GF(5^3) 3036 sage: Frob = k.frobenius_endomorphism() 3037 sage: S.<x> = k['x',Frob] 3038 sage: a = x + t 3039 sage: b = a^10 3040 sage: b == a*a*a*a*a*a*a*a*a*a 3041 True 3042 3043 sage: modulus = x^3 + t*x^2 + (t+3)*x - 2 3044 sage: bmod = a._pow_(10,modulus); bmod 3045 (t^2 + t)*x^2 + (3*t^2 + 1)*x + t^2 + t 3046 sage: bmod == b.rem(modulus) 3047 True 3048 """ 3049 return self._pow_(exp,modulus) 3050 3051 3052 def make_generic_skew_polynomial(parent, coeffs): 3053 return parent(coeffs) 3054 3055 3056 cdef class ConstantSkewPolynomialSection(Map): 3057 cpdef Element _call_(self, x): 3058 if x.degree() <= 0: 3059 try: 3060 return <Element>(x.constant_coefficient()) 3061 except AttributeError: 3062 return <Element>((<SkewPolynomial>x).constant_coefficient()) 3063 else: 3064 raise TypeError("not a constant polynomial") 3065 3066 3067 cdef class SkewPolynomialBaseringInjection(RingHomomorphism): 3068 def __init__(self, domain, codomain): 3069 assert codomain.base_ring() is domain, "domain must be basering" 3070 RingHomomorphism.__init__(self, Hom(domain,codomain)) 3071 self._an_element = codomain.gen() 3072 self._repr_type_str = "Polynomial base injection" 3073 self._new_constant_poly_ = self._an_element._new_constant_poly 3074 3075 cpdef Element _call_(self, x): 3076 return self._new_constant_poly_(x, self._codomain) 3077 3078 cpdef Element _call_with_args(self, x, args=(), kwds={}): 3079 try: 3080 return self._codomain._element_constructor_(x, *args, **kwds) 3081 except AttributeError: 3082 # if there is no element constructor, there is a custom call method. 3083 return self._codomain(x, *args, **kwds) 3084 3085 def section(self): 3086 return ConstantSkewPolynomialSection(self._codomain, self._domain) -
new file sage/rings/polynomial/skew_polynomial_finite_field.pxd
diff --git a/sage/rings/polynomial/skew_polynomial_finite_field.pxd b/sage/rings/polynomial/skew_polynomial_finite_field.pxd new file mode 100644
- + 1 from sage.rings.integer cimport Integer 2 3 from skew_polynomial_element cimport SkewPolynomial_generic_dense 4 from sage.matrix.matrix_dense cimport Matrix_dense 5 from skew_polynomial_element cimport CenterSkewPolynomial_generic_dense 6 7 from polynomial_element cimport Polynomial 8 from sage.structure.element cimport RingElement 9 10 11 cdef class SkewPolynomial_finite_field_dense (SkewPolynomial_generic_dense): 12 # cache 13 cdef list _conjugates 14 cdef Polynomial _norm 15 cdef _norm_factor 16 cdef _optbound 17 cdef dict _rdivisors 18 cdef dict _types 19 cdef _factorization 20 21 cdef inline void _init_cache(self) 22 23 # Karatsuba 24 #cpdef RingElement _mul_karatsuba(self, RingElement other, cutoff=*) 25 cpdef SkewPolynomial_finite_field_dense _mul_central(self, SkewPolynomial_finite_field_dense right) 26 cpdef RingElement _mul_(self, RingElement right) 27 cpdef rquo_rem_karatsuba(self, RingElement other, cutoff=*) 28 29 cdef SkewPolynomial_finite_field_dense _rgcd(self,SkewPolynomial_finite_field_dense other) 30 cdef SkewPolynomial_finite_field_dense _coeff_llcm(self, SkewPolynomial_finite_field_dense other) 31 32 # Inplace functions 33 cdef void _inplace_lrem(self, SkewPolynomial_finite_field_dense other) 34 cdef void _inplace_rrem(self, SkewPolynomial_finite_field_dense other) 35 cdef void _inplace_lfloordiv(self, SkewPolynomial_finite_field_dense other) 36 cdef void _inplace_rfloordiv(self, SkewPolynomial_finite_field_dense other) 37 cdef void _inplace_pow_mod(self, Integer n, SkewPolynomial_finite_field_dense mod) 38 cdef void _inplace_lmonic(self) 39 cdef void _inplace_rmonic(self) 40 cdef void _inplace_rgcd(self,SkewPolynomial_finite_field_dense other) 41 cdef Py_ssize_t _val_inplace_unit(self) 42 cdef SkewPolynomial_finite_field_dense _rquo_inplace_rem(self, SkewPolynomial_finite_field_dense other) 43 44 # Specific functions 45 cdef Matrix_dense _matphir_c(self) 46 cdef Matrix_dense _matmul_c(self) 47 48 # Finding divisors 49 cdef SkewPolynomial_finite_field_dense _rdivisor_c(P, CenterSkewPolynomial_generic_dense N) 50 51 # Finding factorizations 52 cdef _factor_c(self) 53 cdef _factor_uniform_c(self) 54 55 56 cdef class SkewPolynomial_finite_field_karatsuba: 57 cdef _parent 58 cdef Py_ssize_t _order 59 cdef Py_ssize_t _cutoff 60 cdef RingElement _zero 61 cdef _twist 62 cdef char _algo_matrix 63 cdef RingElement _t 64 cdef Matrix_dense _T, _Tinv 65 66 cdef list mul_step (self, list x, list y) 67 cdef list mul_step_matrix(self, list x, list y) 68 cdef list mul_iter(self, list x, list y, char flag) 69 cdef list _twinv 70 cdef list div_step(self, list a, Py_ssize_t ia, Py_ssize_t da, list b, Py_ssize_t ib, Py_ssize_t db) 71 cdef list div_iter(self, list a, Py_ssize_t ia, Py_ssize_t da, list b, Py_ssize_t ib, Py_ssize_t db) -
new file sage/rings/polynomial/skew_polynomial_finite_field.pyx
diff --git a/sage/rings/polynomial/skew_polynomial_finite_field.pyx b/sage/rings/polynomial/skew_polynomial_finite_field.pyx new file mode 100644
- + 1 r""" 2 This module implements skew polynomials over finite fields. 3 4 Let `k` be a finite field and `\sigma` be a ring automorphism 5 of `k` (i.e. a power of the Frobenius endomorphism). Let 6 Put `S = k[X,\sigma]`: as an addtive group, it is the usual ring 7 of polynomials with coefficients in `k` and the multiplication 8 on `S` is defined by the rule `X * a = \sigma(a) * X`. 9 10 We recall that: 11 12 #. `S` is a left (resp. right) euclidean noncommutative ring 13 14 #. in particular, every left (resp. right) ideal is principal 15 16 Since `k` is a finite field, we have the additional following 17 properties: 18 19 #. the center of `S`, denoted by `Z`, is the univariate 20 polynomial ring over `k^\sigma` (`k` fixed by `\sigma`) 21 in the variable `X^r` where `r` is the order of `\sigma` 22 23 #. `S[1/X]` is an Azumaya algebra over `Z[1/X^r]` (i.e. etale 24 locally `S[1/X]` is a matrix algebra over `Z[1/X^r]`) 25 26 #. in particular, we have a reduced norm map `N` from `S[1/X]` 27 to `Z[1/X^r]` (etale locally, it is the determinant); one 28 can prove that it maps `S` to `Z` 29 30 #. `N` has very good properties regarding to factorizations; in 31 particular: 32 33 #. if `a` is a skew polynomial, `a` always divides `N(a)` 34 35 #. if `a` is a skew polynomial, any factorization of `N(a)` 36 (in any order) lifts to a factorization of `a` (and we 37 have a precise control on the number of such lifts); as 38 a consequence there is an explicit (but complicated) 39 formula counting the number of factorizations of a skew 40 polynomial. 41 42 43 .. TODO:: 44 45 Try to replace as possible ``finite field`` by ``field 46 endowed with a finite order twist morphism``. It may cause 47 new phenomena due to the non trivality of the Brauer group. 48 49 EXAMPLES:: 50 51 We illustrate some properties listed above:: 52 53 sage: k.<t> = GF(5^3) 54 sage: Frob = k.frobenius_endomorphism() 55 sage: S.<x> = k['x',Frob]; S 56 Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5 57 sage: Z = S.center(); Z 58 Center of Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5: 59 Univariate Polynomial Ring in (x^3) over Finite Field of size 5 60 sage: a = x^5 + (2*t^2 + t + 1)*x^4 + (3*t^2 + 3*t + 2)*x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 61 62 We compute the reduced norm of `a`:: 63 64 sage: N = a.reduced_norm(); N 65 (x^3)^5 + 3*(x^3)^4 + 2*(x^3)^3 + 3*(x^3) + 4 66 67 Note that the parent of `N` is the center `Z` (and not `S`):: 68 69 sage: N.parent() 70 Center of Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5: 71 Univariate Polynomial Ring in (x^3) over Finite Field of size 5 72 sage: N.parent() == Z 73 True 74 75 sage: S(N) # coercion of N into S 76 x^15 + 3*x^12 + 2*x^9 + 3*x^3 + 4 77 78 `N` is a mutiple of `a`:: 79 80 sage: S(N).is_divisible_by(a) 81 True 82 sage: S(N).is_divisible_by(a,side=Left) 83 True 84 85 .. NOTE:: 86 87 We really need to coerce first `N` into `S`. Otherwise an 88 error occurs:: 89 90 sage: N.is_divisible_by(a) 91 Traceback (most recent call last): 92 ... 93 AttributeError: 'sage.rings.polynomial.skew_polynomial_element.CenterSkewPolynomial_generic_dense' object has no attribute 'is_divisible_by' 94 95 As a polynomial in `x^r` (here `r = 3`), ``N`` factors as a produit 96 of two irreducible polynomials:: 97 98 sage: N.factor() 99 ((x^3)^2 + (x^3) + 1) * ((x^3)^3 + 2*(x^3)^2 + 4*(x^3) + 4) 100 101 And so does `a`:: 102 103 sage: F = a.factor(); F 104 (x^3 + (t^2 + 2*t + 1)*x^2 + (4*t^2 + t + 2)*x + 3*t^2 + t + 4) * (x^2 + (t^2 + 4*t)*x + 2*t) 105 106 We can check that each of these two factors of `a` corresponds 107 to a factor of `N`:: 108 109 sage: F[0][0].reduced_norm() 110 (x^3)^3 + 2*(x^3)^2 + 4*(x^3) + 4 111 sage: F[1][0].reduced_norm() 112 (x^3)^2 + (x^3) + 1 113 114 Actually, `a` has exactly two factorizations corresponding to the 115 two possible orderings of the irreducible factors of `N`:: 116 117 sage: a.count_factorizations() 118 2 119 sage: for F in a.factorizations(): print F 120 (x^3 + (t^2 + 2*t + 1)*x^2 + (4*t^2 + t + 2)*x + 3*t^2 + t + 4) * (x^2 + (t^2 + 4*t)*x + 2*t) 121 (x^2 + (2*t + 2)*x + t + 1) * (x^3 + (3*t^2 + 1)*x^2 + (2*t + 4)*x + 4*t^2 + 3*t + 4) 122 123 AUTHOR:: 124 125 - Xavier Caruso (2012-06-29) 126 """ 127 128 ############################################################################# 129 # Copyright (C) 2012 Xavier Caruso <xavier.caruso@normalesup.org> 130 # 131 # Distributed under the terms of the GNU General Public License (GPL) 132 # 133 # http://www.gnu.org/licenses/ 134 #**************************************************************************** 135 136 137 include "../../ext/stdsage.pxi" 138 include "../../ext/interrupt.pxi" 139 140 import operator, copy, re 141 142 from sage.rings.integer cimport Integer 143 from sage.rings.integer_ring import ZZ 144 from sage.functions.other import ceil 145 146 from sage.structure.parent cimport Parent 147 148 from skew_polynomial_element cimport SkewPolynomial_generic_dense 149 from sage.matrix.constructor import matrix 150 from sage.matrix.constructor import zero_matrix 151 from sage.matrix.matrix_dense cimport Matrix_dense 152 from sage.matrix.matrix_space import MatrixSpace 153 from sage.structure.factorization import Factorization 154 155 from skew_polynomial_element cimport SkewPolynomial 156 from polynomial_ring import PolynomialRing_general 157 from polynomial_ring_constructor import PolynomialRing 158 159 from sage.rings.ring cimport Ring 160 from sage.structure.element cimport RingElement 161 162 from sage.structure.side import Left, Right 163 from sage.rings.infinity import Infinity 164 165 from sage.combinat.q_analogues import q_jordan 166 from sage.functions.other import factorial 167 from sage.combinat.permutation import Permutations 168 from sage.combinat.partition import Partition 169 from sage.misc.mrange import xmrange_iter 170 171 172 cdef class SkewPolynomial_finite_field_dense (SkewPolynomial_generic_dense): 173 def __init__(self, parent, x=None, int check=1, is_gen=False, int construct=0, **kwds): 174 SkewPolynomial_generic_dense.__init__ (self, parent, x, check, is_gen, construct, **kwds) 175 self._init_cache() 176 177 178 cdef inline void _init_cache(self): 179 """ 180 Initialize cached variables (set them to None). 181 """ 182 self._conjugates = [ self.__coeffs ] 183 self._norm = None 184 self._norm_factor = None 185 self._optbound = None 186 self._rdivisors = None 187 self._types = None 188 self._factorization = None 189 190 191 cdef SkewPolynomial _new_c(self, list coeffs, Parent P, char check=0): 192 """ 193 Fast creation of a new skew polynomial 194 """ 195 cdef SkewPolynomial_finite_field_dense f = <SkewPolynomial_finite_field_dense>PY_NEW_SAME_TYPE(self) 196 f._parent = P 197 f.__coeffs = coeffs 198 f._init_cache() 199 if check: 200 f.__normalize() 201 return f 202 203 204 # Skew multiplication 205 # ------------------- 206 207 def _mul_karatsuba(self,right,cutoff=None): 208 """ 209 Karatsuba multiplication 210 211 INPUT: 212 213 - ``right`` -- an other skew polynomial in the same ring 214 215 - ``cutoff`` -- ``None``, an integer or Infinity (default: None) 216 217 .. WARNING:: 218 219 ``cutoff`` need to be greater than or equal to the order of the 220 twist map acting on the base ring of the underlying skew polynomial 221 ring. 222 223 OUTPUT: 224 225 The result of the product self*right (computed by a variant of 226 Karatsuba`s algorithm) 227 228 .. NOTE:: 229 230 if ``cutoff`` is None, use the default cutoff which is the 231 maximum between 150 and the order of the twist map. 232 233 EXAMPLES:: 234 235 sage: k.<t> = GF(5^3) 236 sage: Frob = k.frobenius_endomorphism() 237 sage: S.<x> = k['x',Frob] 238 sage: a = S.random_element(5000) 239 sage: b = S.random_element(5000) 240 sage: timeit("c = a._mul_karatsuba(b)") # random, long time 241 5 loops, best of 3: 659 ms per loop 242 sage: timeit("c = a._mul_classical(b)") # random, long time 243 5 loops, best of 3: 1.9 s per loop 244 sage: a._mul_karatsuba(b) == a._mul_classical(b) 245 True 246 247 The operator ``*`` performs Karatsuba multiplication:: 248 249 sage: timeit("c = a*b") # random, long time 250 5 loops, best of 3: 653 ms per loop 251 """ 252 karatsuba_class = self._parent._karatsuba_class 253 if cutoff != None: 254 save_cutoff = karatsuba_class.get_cutoff() 255 karatsuba_class.set_cutoff(cutoff) 256 res = karatsuba_class.mul(self,right) 257 if cutoff != None: 258 karatsuba_class.set_cutoff(save_cutoff) 259 return res 260 261 262 def _mul_karatsuba_matrix(self,right): 263 """ 264 Karatsuba multiplication with multiplication step based 265 on an isomorphism between a quotient of the underlying 266 skew polynomial ring and a ring of matrices. 267 268 INPUT: 269 270 - ``right`` -- an other skew polynomial in the same ring 271 272 OUTPUT: 273 274 The result of the product self*right 275 276 EXAMPLES:: 277 278 sage: k.<t> = GF(5^3) 279 sage: Frob = k.frobenius_endomorphism() 280 sage: S.<x> = k['x',Frob] 281 sage: a = S.random_element(degree=20) 282 sage: b = S.random_element(degree=20) 283 sage: a._mul_karatsuba_matrix(b) == a*b 284 True 285 286 This routine is only efficient when the twisting map (here 287 ``Frob``) has a large order `r` and the degrees of ``self`` 288 and ``other`` have a very special shape (just below a power 289 of `2` times `r * floor(r/2)`):: 290 291 sage: k.<t> = GF(5^40) 292 sage: Frob = k.frobenius_endomorphism() 293 sage: S.<x> = k['x',Frob] 294 sage: a = S.random_element(degree=799) 295 sage: b = S.random_element(degree=799) 296 sage: timeit("c = a*b") # random, long time 297 5 loops, best of 3: 23.1 s per loop 298 sage: timeit("c = a._mul_karatsuba_matrix(b)") # random, long time 299 5 loops, best of 3: 12.2 s per loop 300 """ 301 karatsuba_class = self._parent._karatsuba_class 302 return karatsuba_class.mul_matrix(self,right) 303 304 305 cpdef SkewPolynomial_finite_field_dense _mul_central(self, SkewPolynomial_finite_field_dense right): 306 r""" 307 Return self * right 308 309 .. WARNING:: 310 311 Do you use this function! It is very slow due to a quite 312 slow interface with ``polynomial_zz_pex``. 313 314 ALGORITHM:: 315 316 Notations:: 317 318 - `S` is the underlyling skew polynomial ring 319 320 - `x` is the variable on `S` 321 322 - `k` is the base ring of `S` (it is a finite field) 323 324 - `\sigma` is the twisting automorphism acting on `k` 325 326 - `r` is the order of `\sigma` 327 328 - `t` is a generator of `k` over `k^\sigma` 329 330 #. We decompose the polynomial ``right`` as follows:: 331 332 .. MATH:: 333 334 right = \sum_{i=0}^{r-1} \sum_{j=0}^{r-1} y_{i,j} t^j x^i 335 336 where `y_{i,j}` are polynomials in the center `k^\sigma[x^r]`. 337 338 #. We compute all products `z_{i,j} = left * y_{i,j}`; since 339 all `y_{i,j}` lie in the center, we can compute all these 340 products as if `left` was a commutative polynomial (and we 341 can therefore use fast algorithms like FFT and/or fast 342 implementations) 343 344 #. We compute and return the sum 345 346 .. MATH:: 347 348 \sum_{i=0}^{r-1} \sum_{j=0}^{r-1} z_{i,j} t^j x^i 349 350 EXAMPLES:: 351 352 sage: k.<t> = GF(5^3) 353 sage: Frob = k.frobenius_endomorphism() 354 sage: S.<x> = k['x',Frob] 355 sage: a = S.random_element(degree=10) 356 sage: b = S.random_element(degree=10) 357 sage: a._mul_central(b) == a*b 358 True 359 360 TESTS:: 361 362 Here is an example where `k^\sigma` is not a prime field:: 363 364 sage: k.<t> = GF(5^6) 365 sage: Frob = k.frobenius_endomorphism(2) 366 sage: S.<x> = k['x',Frob] 367 sage: a = S.random_element(degree=10) 368 sage: b = S.random_element(degree=10) 369 sage: a._mul_central(b) == a*b 370 True 371 """ 372 skew_ring = self._parent 373 base_ring = skew_ring.base_ring() 374 commutative_ring = PolynomialRing(skew_ring.base_ring(),name='x') 375 cdef RingElement c 376 cdef RingElement zero = base_ring(0) 377 cdef Py_ssize_t i, j, k 378 cdef Py_ssize_t order = skew_ring._order 379 cdef Py_ssize_t degree = base_ring.degree() 380 381 left = commutative_ring(self.__coeffs) 382 cdef list y = [ c.polynomial() for c in right.__coeffs ] 383 cdef Py_ssize_t leny = len(y) 384 cdef list yc = leny * [zero] 385 cdef list res = (leny + len(self.__coeffs) - 1) * [zero] 386 cdef list term 387 cdef list twist = [ base_ring.gen() ] 388 for i from 0 <= i < order-1: 389 twist.append(skew_ring.twist_map(1)(twist[i])) 390 for i from 0 <= i < order: 391 for j from 0 <= j < degree: 392 for k from i <= k < leny by order: 393 yc[k] = y[k][j] 394 term = (left * commutative_ring(yc)).list() 395 for k from i <= k < len(term): 396 res[k] += term[k] * twist[(k-i)%order]**j 397 for k from i <= k < leny by order: 398 yc[k] = zero 399 return self._new_c(res,skew_ring,1) 400 401 402 cpdef RingElement _mul_(self, RingElement right): 403 """ 404 Compute self * right (in this order) 405 406 .. NOTE:: 407 408 Use skew Karatsuba's algorithm for skew 409 polynomials of large degrees. 410 411 INPUT: 412 413 - right -- a skew polynomial in the same ring 414 415 TESTS:: 416 417 sage: k.<t> = GF(5^3) 418 sage: Frob = k.frobenius_endomorphism() 419 sage: S.<x> = k['x',Frob] 420 sage: a = x^2 + t; a 421 x^2 + t 422 sage: b = x^2 + (t + 1)*x; b 423 x^2 + (t + 1)*x 424 sage: a * b 425 x^4 + (3*t^2 + 2)*x^3 + t*x^2 + (t^2 + t)*x 426 sage: a * b == b * a 427 False 428 """ 429 return self._parent._karatsuba_class.mul(self,right) 430 431 432 def _mul_classical(self,right): 433 """ 434 Compute self * right (in this order) using the 435 skew SchoolBook algorithm. 436 437 INPUT: 438 439 - ``right`` -- a skew polynomial in the same ring 440 441 TESTS:: 442 443 sage: k.<t> = GF(5^3) 444 sage: Frob = k.frobenius_endomorphism() 445 sage: S.<x> = k['x',Frob] 446 sage: a = x^2 + t; a 447 x^2 + t 448 sage: b = x^2 + (t + 1)*x; b 449 x^2 + (t + 1)*x 450 sage: a._mul_classical(b) 451 x^4 + (3*t^2 + 2)*x^3 + t*x^2 + (t^2 + t)*x 452 sage: a * b == b * a 453 False 454 """ 455 karatsuba_class = self._parent._karatsuba_class 456 save_cutoff = karatsuba_class.get_cutoff() 457 karatsuba_class.set_cutoff(Infinity) 458 res = karatsuba_class.mul(self,right) 459 karatsuba_class.set_cutoff(save_cutoff) 460 return res 461 462 463 cpdef rquo_rem_karatsuba(self, RingElement other, cutoff=None): 464 """ 465 Right euclidean division based on Karatsuba's algorithm. 466 467 DO NOT USE THIS! It is not efficient for usual degrees! 468 469 .. TODO:: 470 471 Try to understand why... 472 473 TESTS:: 474 475 sage: k.<t> = GF(5^3) 476 sage: Frob = k.frobenius_endomorphism() 477 sage: S.<x> = k['x',Frob] 478 479 sage: a = S.random_element(2000) 480 sage: b = S.random_element(1000) 481 sage: timeit("q,r = a.rquo_rem_karatsuba(b)") # random, long time 482 5 loops, best of 3: 104 ms per loop 483 sage: timeit("q,r = a.rquo_rem(b)") # random, long time 484 5 loops, best of 3: 79.6 ms per loop 485 sage: a.rquo_rem(b) == a.rquo_rem_karatsuba(b) 486 True 487 488 sage: a = S.random_element(10000) 489 sage: b = S.random_element(5000) 490 sage: timeit("q,r = a.rquo_rem_karatsuba(b)") # random, long time 491 5 loops, best of 3: 1.79 s per loop 492 sage: timeit("q,r = a.rquo_rem(b)") # random, long time 493 5 loops, best of 3: 1.93 s per loop 494 sage: a.rquo_rem(b) == a.rquo_rem_karatsuba(b) 495 True 496 """ 497 karatsuba_class = self._parent._karatsuba_class 498 if cutoff != None: 499 save_cutoff = karatsuba_class.get_cutoff() 500 karatsuba_class.set_cutoff(cutoff) 501 res = karatsuba_class.div(self,other) 502 if cutoff != None: 503 karatsuba_class.set_cutoff(save_cutoff) 504 return res 505 506 507 # We improve some functions 508 # ------------------------- 509 510 def rquo_rem(self,other): 511 """ 512 DEFINITION: 513 514 Let `a` and `b` be two skew polynomials over the same 515 ring. The *right euclidean division* of `a` by `b` is 516 a couple `(q,r)` such that 517 518 - `a = q*b + r` 519 520 - the degree of `r` is less than the degree of `b` 521 522 `q` (resp. `r`) is called the *quotient* (resp. the 523 remainder) of this euclidean division. 524 525 If the leading coefficient of `b` is a unit (e.g. if 526 `b` is monic) then `q` and `r` exist and are unique. 527 528 INPUT: 529 530 - ``other`` -- a skew polynomial ring over the same 531 base ring 532 533 OUTPUT: 534 535 - the quotient and the remainder of the left euclidean 536 division of this skew polynomial by ``other`` 537 538 .. NOTE:: 539 540 Doesn't work if the leading coefficient of the divisor 541 is not a unit. 542 543 EXAMPLES:: 544 545 sage: R.<t> = ZZ[] 546 sage: sigma = R.hom([t+1]) 547 sage: S.<x> = R['x',sigma] 548 sage: a = S.random_element(degree=4); a 549 t^2*x^4 + (-12*t^2 - 2*t - 1)*x^3 + (-95*t^2 + t + 2)*x^2 + (-t^2 + t)*x + 2*t - 8 550 sage: b = S.random_element(monic=True); b 551 x^2 + (4*t^2 - t - 2)*x - t^2 + t - 1 552 sage: q,r = a.rquo_rem(b) 553 sage: a == q*b + r 554 True 555 556 The leading coefficient of the divisor need to be invertible:: 557 558 sage: c = S.random_element(); c 559 (-4*t^2 + t)*x^2 - 2*t^2*x + 5*t^2 - 6*t - 4 560 sage: a.rquo_rem(c) 561 Traceback (most recent call last): 562 ... 563 NotImplementedError: the leading coefficient of the divisor is not invertible 564 """ 565 cdef list a = self.list() 566 cdef list b = other.list() 567 cdef Py_ssize_t i, j 568 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 569 parent = self._parent 570 if db < 0: 571 raise ZeroDivisionError 572 if da < db: 573 return parent(0), self 574 cdef RingElement inv = ~b[db] 575 cdef list q = [ ] 576 cdef Py_ssize_t order = parent._order 577 cdef list twinv = [ inv ], twb = [ b ] 578 cdef RingElement c, x 579 for i from 0 <= i < min(da-db,order-1): 580 twinv.append(parent.twist_map()(twinv[i])) 581 twb.append([ parent.twist_map()(x) for x in twb[i] ]) 582 for i from da-db >= i >= 0: 583 c = twinv[i%order] * a[i+db] 584 for j from 0 <= j < db: 585 a[i+j] -= c * twb[i%order][j] 586 q.append(c) 587 q.reverse() 588 return parent(q), parent(a[:db]) 589 590 591 cdef SkewPolynomial_finite_field_dense _rgcd(self,SkewPolynomial_finite_field_dense other): 592 """ 593 Fast right gcd. 594 """ 595 cdef SkewPolynomial_finite_field_dense A = self 596 cdef SkewPolynomial_finite_field_dense B = other 597 cdef SkewPolynomial_finite_field_dense swap 598 if len(B.__coeffs): 599 A = <SkewPolynomial_finite_field_dense>self._new_c(A.__coeffs[:],A._parent) 600 B = <SkewPolynomial_finite_field_dense>B._new_c(B.__coeffs[:],B._parent) 601 while len(B.__coeffs): 602 A._inplace_rrem(B) 603 swap = A; A = B; B = swap 604 return A 605 else: 606 return self 607 608 609 cdef void _inplace_rmul(self, SkewPolynomial_generic_dense right): 610 """ 611 Replace self by self*right 612 """ 613 self.__coeffs = self._parent._karatsuba_class.mul_list(self.__coeffs,right.__coeffs) 614 self._init_cache() 615 616 617 cdef void _inplace_lmul(self, SkewPolynomial_generic_dense left): 618 """ 619 Replace self by left*self 620 """ 621 self.__coeffs = self._parent._karatsuba_class.mul_list(left.__coeffs,self.__coeffs) 622 self._init_cache() 623 624 625 cpdef _pow_(self,right,modulus=None,side=Right): 626 """ 627 INPUT: 628 629 - ``exp`` -- an Integer 630 631 - ``modulus`` -- a skew polynomial over the same ring (default: None) 632 633 - ``side`` -- ``Left`` or ``Right`` (default: Right) 634 635 OUTPUT: 636 637 If ``modulus`` is None, return ``self**exp``. 638 639 Otherwise, return the remainder of self**exp in the ``side`` 640 euclidean division by ``modulus``. 641 642 REMARK: 643 644 The quotient of the underlying skew polynomial ring by the 645 principal ideal generated by ``modulus`` is in general *not* 646 a ring. 647 648 As a consequence, Sage first computes exactly ``self**exp`` 649 and then reduce it modulo ``modulus``. 650 651 However, if the base ring is a finite field, Sage uses the 652 following optimized algorithm: 653 654 #. One first compute a central skew polynomial `N` which is 655 divisible by ``modulus``. (Since `N` lies in center, the 656 quotient `K[X,\sigma]/N` inherits a ring structure.) 657 658 #. One compute ``self**exp`` in the quotient ring `K[X,\sigma]/N` 659 660 #. One reduce modulo ``modulus`` the result computed in the 661 previous step 662 663 EXAMPLES:: 664 665 sage: k.<t> = GF(5^3) 666 sage: Frob = k.frobenius_endomorphism() 667 sage: S.<x> = k['x',Frob] 668 sage: a = x + t 669 sage: b = a^10 # short form for ``a._pow_(10)`` 670 sage: b == a*a*a*a*a*a*a*a*a*a 671 True 672 673 sage: modulus = x^3 + t*x^2 + (t+3)*x - 2 674 sage: br = a._pow_(10,modulus); br 675 (t^2 + t)*x^2 + (3*t^2 + 1)*x + t^2 + t 676 sage: br == b.rem(modulus) 677 True 678 sage: bl = a._pow_(10,modulus,side=Left); bl 679 (4*t^2 + 2*t + 3)*x^2 + (3*t^2 + 1)*x + 2*t + 3 680 sage: bl == b.rem(modulus,side=Left) 681 True 682 683 sage: a._pow_(10^100,modulus) # rather fast 684 (3*t^2 + 3)*x^2 + (t^2 + 2*t + 4)*x + 4*t^2 + 2*t + 1 685 """ 686 cdef SkewPolynomial_finite_field_dense r 687 688 if not isinstance(right, Integer) or isinstance(right, int): 689 try: 690 right = Integer(right) 691 except TypeError: 692 raise TypeError("non-integral exponents not supported") 693 694 if self.degree() <= 0: 695 return self.parent()(self[0]**right) 696 if right == 0: 697 return self.parent()(1) 698 if right < 0: 699 return (~self).pow(-right,modulus,side=side) 700 701 if self == self.parent().gen(): # special case x**n should be faster! 702 P = self.parent() 703 R = P.base_ring() 704 v = [R.zero_element()]*right + [R.one_element()] 705 r = <SkewPolynomial_generic_dense>self._parent(v) 706 sig_on() 707 if modulus: 708 if side is Right: 709 r._inplace_rrem(modulus) 710 else: 711 r._inplace_lrem(modulus) 712 r._init_cache() 713 sig_off() 714 return r 715 716 mod = modulus 717 if not modulus is None: 718 try: 719 mod = self.parent()(mod.bound()) 720 except NotImplementedError: 721 mod = None 722 r = <SkewPolynomial_generic_dense>self._new_c(self.__coeffs,self._parent) 723 sig_on() 724 if mod: 725 r._inplace_pow_mod(right,mod) 726 else: 727 r._inplace_pow(right) 728 if (not modulus is None) and modulus != mod: 729 if side is Right: 730 r._inplace_rrem(modulus) 731 else: 732 r._inplace_lrem(modulus) 733 r._init_cache() 734 sig_off() 735 return r 736 737 738 # Inplace functions 739 # ----------------- 740 741 #cdef void _inplace_conjugate(self,n): 742 # cdef Py_ssize_t i 743 # cdef Morphism twist = <Morphism>self._parent.twist_map(n) 744 # cdef RingElement x 745 # for i from 0 <= i < len(self.__coeffs): 746 # x = twist(self.__coeffs[i]) 747 # self.__coeffs[i] = x 748 749 750 cdef void _inplace_lrem(self, SkewPolynomial_finite_field_dense other): 751 """ 752 Replace self by the remainder in the left euclidean division 753 of self by other (only for internal use). 754 """ 755 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 756 cdef list b = (<SkewPolynomial_finite_field_dense>other).__coeffs 757 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 758 cdef Py_ssize_t i, j 759 cdef RingElement c, inv 760 parent = self._parent 761 if db < 0: 762 raise ZeroDivisionError 763 if da >= db: 764 inv = ~b[db] 765 for i from da-db >= i >= 0: 766 c = parent.twist_map(-db)(inv*a[i+db]) 767 for j from 0 <= j < db: 768 a[i+j] -= b[j] * parent.twist_map(j)(c) 769 del a[db:] 770 self.__normalize() 771 self._init_cache() 772 773 774 cdef void _inplace_rrem(self, SkewPolynomial_finite_field_dense other): 775 """ 776 Replace self by the remainder in the right euclidean division 777 of self by other (only for internal use). 778 """ 779 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 780 cdef list b = (<SkewPolynomial_finite_field_dense>other).__coeffs 781 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 782 cdef Py_ssize_t i, j, order 783 cdef RingElement c, x, inv 784 cdef list twinv, twb 785 parent = self._parent 786 if db < 0: 787 raise ZeroDivisionError 788 if da >= db: 789 order = parent._order 790 inv = ~b[db] 791 twinv = [ inv ] 792 for i from 0 <= i < min(da-db,order-1): 793 twinv.append(parent.twist_map()(twinv[i])) 794 twb = (<SkewPolynomial_finite_field_dense>other)._conjugates 795 for i from len(twb)-1 <= i < min(da-db,order-1): 796 twb.append([ parent.twist_map()(x) for x in twb[i] ]) 797 for i from da-db >= i >= 0: 798 c = twinv[i%order] * a[i+db] 799 for j from 0 <= j < db: 800 a[i+j] -= c * twb[i%order][j] 801 del a[db:] 802 self.__normalize() 803 self._init_cache() 804 805 806 cdef void _inplace_lfloordiv(self, SkewPolynomial_finite_field_dense other): 807 """ 808 Replace self by the quotient in the left euclidean division 809 of self by other (only for internal use). 810 """ 811 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 812 cdef list b = (<SkewPolynomial_finite_field_dense>other).__coeffs 813 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 814 cdef Py_ssize_t i, j, deb 815 cdef RingElement c, inv 816 parent = self._parent 817 if db < 0: 818 raise ZeroDivisionError 819 if da < db: 820 (<SkewPolynomial_finite_field_dense>self).__coeffs = [ ] 821 else: 822 inv = ~b[db] 823 for i from da-db >= i >= 0: 824 c = a[i+db] = parent.twist_map(-db)(inv*a[i+db]) 825 if i < db: deb = db 826 else: deb = i 827 for j from deb <= j < db+i: 828 a[j] -= b[j-i] * parent.twist_map(j-i)(c) 829 del a[:db] 830 self.__normalize() 831 self._init_cache() 832 833 834 cdef void _inplace_rfloordiv(self, SkewPolynomial_finite_field_dense other): 835 """ 836 Replace self by the quotient in the right euclidean division 837 of self by other (only for internal use). 838 """ 839 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 840 cdef list b = (<SkewPolynomial_finite_field_dense>other).__coeffs 841 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 842 cdef Py_ssize_t i, j, deb, order 843 cdef RingElement c, x, inv 844 parent = self._parent 845 if db < 0: 846 raise ZeroDivisionError 847 if da < db: 848 (<SkewPolynomial_finite_field_dense>self).__coeffs = [ ] 849 else: 850 order = parent._order 851 inv = ~b[db] 852 twinv = [ inv ] 853 for i from 0 <= i < min(da-db,order-1): 854 twinv.append(parent.twist_map()(twinv[i])) 855 twb = (<SkewPolynomial_finite_field_dense>other)._conjugates 856 for i from len(twb)-1 <= i < min(da-db,order-1): 857 twb.append([ parent.twist_map()(x) for x in twb[i] ]) 858 for i from da-db >= i >= 0: 859 c = a[i+db] = twinv[i%order] * a[i+db] 860 if i < db: deb = db 861 else: deb = i 862 for j from deb <= j < db+i: 863 a[j] -= c * twb[i%order][j-i] 864 del a[:db] 865 self.__normalize() 866 self._init_cache() 867 868 869 cdef void _inplace_pow_mod(self, Integer n, SkewPolynomial_finite_field_dense mod): 870 """ 871 Replace self by the remainder in the euclidean 872 division of self**n by mod (only for internal use). 873 874 .. WARNING:: 875 876 Assume that mod is central 877 """ 878 while n&1 == 0: 879 self._inplace_rmul(self) 880 self._inplace_rrem(mod) 881 n = n >> 1 882 cdef SkewPolynomial_finite_field_dense selfpow = <SkewPolynomial_finite_field_dense>self._new_c(self.__coeffs[:], self._parent) 883 n = n >> 1 884 while n != 0: 885 selfpow._inplace_rmul(selfpow) 886 selfpow._inplace_rrem(mod) 887 if n&1 == 1: 888 self._inplace_rmul(selfpow) 889 self._inplace_rrem(mod) 890 n = n >> 1 891 892 893 cdef void _inplace_lmonic(self): 894 """ 895 Replace self by ``self.lmonic()`` (only for internal use). 896 """ 897 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 898 cdef Py_ssize_t da = len(a)-1, i 899 cdef RingElement inv = ~a[da] 900 parent = self._parent 901 a[da] = parent.base_ring()(1) 902 for i from 0 <= i < da: 903 a[i] *= parent.twist_map(i-da)(inv) 904 self._init_cache() 905 906 907 cdef void _inplace_rmonic(self): 908 """ 909 Replace self by ``self.rmonic()`` (only for internal use). 910 """ 911 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 912 cdef Py_ssize_t da = len(a)-1, i 913 cdef RingElement inv = ~a[da] 914 a[da] = self._parent.base_ring()(1) 915 for i from 0 <= i < da: 916 a[i] *= inv 917 self._init_cache() 918 919 920 cdef void _inplace_rgcd(self,SkewPolynomial_finite_field_dense other): 921 """ 922 Replace self by its right gcd with other (only for internal use). 923 """ 924 cdef SkewPolynomial_finite_field_dense B 925 cdef list swap 926 if len(other.__coeffs): 927 B = <SkewPolynomial_finite_field_dense>self._new_c(other.__coeffs[:],other._parent) 928 while len(B.__coeffs): 929 B._conjugates = [ B.__coeffs ] 930 self._inplace_rrem(B) 931 swap = self.__coeffs 932 self.__coeffs = B.__coeffs 933 B.__coeffs = swap 934 self._init_cache() 935 936 937 cdef SkewPolynomial_finite_field_dense _rquo_inplace_rem(self, SkewPolynomial_finite_field_dense other): 938 """ 939 Replace self by the remainder in the right euclidean division 940 of self by other and return the quotient (only for internal use). 941 """ 942 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 943 cdef list b = (<SkewPolynomial_finite_field_dense>other).__coeffs 944 cdef Py_ssize_t da = len(a)-1, db = len(b)-1 945 cdef Py_ssize_t i, j 946 cdef RingElement c, inv 947 cdef list q 948 parent = self._parent 949 if db < 0: 950 raise ZeroDivisionError 951 if da < db: 952 return self._new_c([],self._parent) 953 inv = ~b[db] 954 q = [ ] 955 for i from da-db >= i >= 0: 956 c = parent.twist_map(i)(inv) * a[i+db] 957 q.append(c) 958 for j from 0 <= j < db: 959 a[i+j] -= c * parent.twist_map(i)(b[j]) 960 del a[db:] 961 self.__normalize() 962 self._init_cache() 963 q.reverse() 964 return self._new_c(q,self._parent) 965 966 967 cdef Py_ssize_t _val_inplace_unit(self): 968 """ 969 Return `v` the valuation of self and replace self by 970 self >> v (only for internal use). 971 """ 972 cdef list a = (<SkewPolynomial_finite_field_dense>self).__coeffs 973 cdef Py_ssize_t val = 0 974 if len(a) < 0: 975 return -1 976 while a[0].is_zero(): 977 del a[0] 978 val += 1 979 self._init_cache() 980 return val 981 982 983 # lowest common divisor: other functions 984 # -------------------------------------- 985 986 cdef SkewPolynomial_finite_field_dense _coeff_llcm(self, SkewPolynomial_finite_field_dense other): 987 """ 988 (only for internal use) 989 """ 990 cdef SkewPolynomial_finite_field_dense A = <SkewPolynomial_finite_field_dense>self 991 cdef SkewPolynomial_finite_field_dense B = <SkewPolynomial_finite_field_dense>other 992 if not len(A.__coeffs) and not len(B.__coeffs): 993 raise ZeroDivisionError 994 A = <SkewPolynomial_finite_field_dense>A._new_c(A.__coeffs[:],A._parent) 995 B = <SkewPolynomial_finite_field_dense>A._new_c(B.__coeffs[:],B._parent) 996 cdef parent = self._parent 997 cdef RingElement one = self.base_ring().one_element() 998 cdef SkewPolynomial_finite_field_dense U = self._new_c([one],parent) 999 cdef SkewPolynomial_finite_field_dense V = self._new_c([],parent) 1000 cdef SkewPolynomial_finite_field_dense Q, R 1001 cdef list swap 1002 while len(B.__coeffs): 1003 Q = A._rquo_inplace_rem(B) 1004 Q._inplace_rmul(V) 1005 U._inplace_sub(Q) 1006 swap = A.__coeffs 1007 A.__coeffs = B.__coeffs 1008 B.__coeffs = swap 1009 swap = U.__coeffs 1010 U.__coeffs = V.__coeffs 1011 V.__coeffs = swap 1012 return V 1013 1014 1015 # Specific functions 1016 # ------------------ 1017 1018 cdef Matrix_dense _matphir_c(self): 1019 r""" 1020 Return the matrix of the multiplication by `X^r` on 1021 the quotient `K[X,\sigma] / K[X,\sigma]*self`. 1022 """ 1023 cdef Py_ssize_t i, j, col, exp, n 1024 cdef Py_ssize_t d = self.degree() 1025 cdef Py_ssize_t r = self.parent()._order 1026 parent = self._parent 1027 cdef k = parent.base_ring() 1028 cdef Matrix_dense phir = <Matrix_dense?>zero_matrix(k,d,d) 1029 cdef RingElement zero = k(0) 1030 cdef RingElement one = k(1) 1031 if r < d: 1032 for i from 0 <= i < d-r: 1033 phir.set_unsafe(r+i,i,one) 1034 col = d-r 1035 exp = d 1036 else: 1037 col = 0 1038 exp = r 1039 cdef SkewPolynomial_finite_field_dense powx = <SkewPolynomial_finite_field_dense>self._new_c([zero,one], parent) 1040 cdef SkewPolynomial_finite_field_dense v 1041 if (exp % 2 == 1): 1042 v = <SkewPolynomial_finite_field_dense>self._new_c([zero,one], parent) 1043 if self.degree() == 1: 1044 v._inplace_rrem(self) 1045 else: 1046 v = <SkewPolynomial_finite_field_dense>self._new_c([one], parent) 1047 exp = exp >> 1 1048 n = 1 1049 while exp != 0: 1050 powx._inplace_lmul(powx.conjugate(n)) 1051 powx._inplace_rrem(self) 1052 n = n << 1 1053 if (exp % 2 == 1): 1054 v = v.conjugate(n) 1055 v._inplace_rmul(powx) 1056 v._inplace_rrem(self) 1057 exp = exp >> 1 1058 l = v.list() 1059 for i from 0 <= i < len(l): 1060 phir.set_unsafe(i,col,l[i]) 1061 for j from col+1 <= j < d: 1062 v <<= 1 1063 v = v.conjugate(1) 1064 v._inplace_rrem(self) 1065 for i from 0 <= i <= v.degree(): 1066 phir.set_unsafe(i,j,v.__coeffs[i]) 1067 return phir 1068 1069 def smurf(self): 1070 return self._matphir_c() 1071 1072 1073 cdef Matrix_dense _matmul_c(self): 1074 r""" 1075 Return the matrix of the multiplication by self on 1076 `K[X,\sigma]` considered as a free module over `K[X^r]` 1077 (here `r` is the order of `\sigma`). 1078 1079 .. WARNING:: 1080 1081 Does not work if self is not monic. 1082 """ 1083 cdef Py_ssize_t i, j, deb, k, r = self.parent()._order 1084 cdef Py_ssize_t d = self.degree () 1085 cdef Ring base_ring = <Ring?>self.parent().base_ring() 1086 cdef RingElement minusone = <RingElement?>base_ring(-1) 1087 cdef RingElement zero = <RingElement?>base_ring(0) 1088 cdef Polk = PolynomialRing (base_ring, 'xr') 1089 cdef Matrix_dense M = <Matrix_dense?>zero_matrix(Polk,r,r) 1090 cdef list l = self.list() 1091 for j from 0 <= j < r: 1092 for i from 0 <= i < r: 1093 if i < j: 1094 pol = [ zero ] 1095 deb = i-j+r 1096 else: 1097 pol = [ ] 1098 deb = i-j 1099 for k from deb <= k <= d by r: 1100 pol.append(l[k]) 1101 M.set_unsafe(i,j,Polk(pol)) 1102 for i from 0 <= i <= d: 1103 l[i] = self._parent.twist_map()(l[i]) 1104 return M 1105 1106 1107 def reduced_norm(self): 1108 r""" 1109 Return the reduced norm of this skew polynomial. 1110 1111 .. NOTE:: 1112 1113 The result is cached. 1114 1115 ALGORITHM: 1116 1117 If `r` (= the order of the twist map) is small compared 1118 to `d` (= the degree of this skew polynomial), the reduced 1119 norm is computed as the determinant of the multiplication 1120 by `P` (= this skew polynomial) acting on `K[X,\sigma]` 1121 (= the underlying skew ring) viewed as a free module of 1122 rank `r` over `K[X^r]`. 1123 1124 Otherwise, the reduced norm is computed as the characteristic 1125 polynomial (considered as a polynomial of the variable `X^r`) 1126 of the left multiplication by `X` on the quotient 1127 `K[X,\sigma] / K[X,\sigma]*P` (which is a `K`-vector space 1128 of dimension `d`). 1129 1130 EXAMPLES:: 1131 1132 sage: k.<t> = GF(5^3) 1133 sage: Frob = k.frobenius_endomorphism() 1134 sage: S.<x> = k['x',Frob] 1135 sage: a = S.random_element(degree=3,monic=True); a 1136 x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 2 1137 sage: N = a.reduced_norm(); N 1138 (x^3)^3 + 4*(x^3)^2 + 4 1139 1140 Note that the parent of `N` is the center of the `S` 1141 (and not `S` itself):: 1142 1143 sage: N.parent() 1144 Center of Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5: 1145 Univariate Polynomial Ring in (x^3) over Finite Field of size 5 1146 sage: N.parent() == S.center() 1147 True 1148 1149 In any case, coercion works fine:: 1150 1151 sage: S(N) 1152 x^9 + 4*x^6 + 4 1153 sage: N + a 1154 x^9 + 4*x^6 + x^3 + (2*t^2 + 3)*x^2 + (4*t^2 + t + 4)*x + 2*t^2 + 1 1155 1156 We check that `N` is a multiple of `a`:: 1157 1158 sage: S(N).is_divisible_by(a) 1159 True 1160 sage: S(N).is_divisible_by(a,side=Left) 1161 True 1162 1163 .. NOTE:: 1164 1165 We really need to coerce first `N` into `S`. Otherwise an 1166 error occurs:: 1167 1168 sage: N.is_divisible_by(a) 1169 Traceback (most recent call last): 1170 ... 1171 AttributeError: 'sage.rings.polynomial.skew_polynomial_element.CenterSkewPolynomial_generic_dense' object has no attribute 'is_divisible_by' 1172 1173 We check that the reduced norm is a multiplicative map:: 1174 1175 sage: a = S.random_element(degree=5) 1176 sage: b = S.random_element(degree=7) 1177 sage: a.reduced_norm() * b.reduced_norm() == (a*b).reduced_norm() 1178 True 1179 """ 1180 if self._norm is None: 1181 center = self.parent().center() 1182 if self.is_zero(): 1183 self._norm = center(0) 1184 else: 1185 section = center._embed_basering.section() 1186 exp = (self.parent().base_ring().cardinality() - 1) / (center.base_ring().cardinality() - 1) 1187 order = self.parent()._order 1188 lc = section(self.leading_coefficient()**exp) 1189 if order < self.degree(): 1190 M = self._matmul_c() 1191 self._norm = center([ lc*section(x) for x in M.determinant().monic().list() ]) 1192 else: 1193 charpoly = self._matphir_c().characteristic_polynomial() 1194 self._norm = center([ lc*section(x) for x in charpoly.list() ]) 1195 return self._norm 1196 1197 1198 def reduced_norm_factor(self): 1199 """ 1200 Return the reduced norm of this polynomial 1201 factorized in the centre. 1202 1203 EXAMPLES: 1204 1205 sage: k.<t> = GF(5^3) 1206 sage: Frob = k.frobenius_endomorphism() 1207 sage: S.<x> = k['x',Frob] 1208 1209 sage: a = (x^2 + 1) * (x+3) 1210 sage: a.reduced_norm_factor() 1211 ((x^3) + 3) * ((x^3) + 2)^2 1212 """ 1213 if self._norm_factor is None: 1214 self._norm_factor = self.reduced_norm().factor() 1215 return self._norm_factor 1216 1217 1218 def is_central(self): 1219 """ 1220 Return True if this skew polynomial lies in the center. 1221 1222 EXAMPLES:: 1223 1224 sage: k.<t> = GF(5^3) 1225 sage: Frob = k.frobenius_endomorphism() 1226 sage: S.<x> = k['x',Frob] 1227 1228 sage: x.is_central() 1229 False 1230 sage: (t*x^3).is_central() 1231 False 1232 sage: (x^6 + x^3).is_central() 1233 True 1234 """ 1235 center = self.parent().center() 1236 try: 1237 center(self) 1238 return True 1239 except ValueError: 1240 return False 1241 1242 1243 def bound(self): 1244 """ 1245 Return a bound of this skew polynomial (i.e. a multiple 1246 of this skew polynomial lying in the center). 1247 1248 .. NOTE:: 1249 1250 Since `b` is central, it divides a skew polynomial 1251 on the left iff it divides it on the right 1252 1253 ALGORITHM: 1254 1255 #. Sage first checks whether ``self`` is itself in the 1256 center. It if is, it returns ``self`` 1257 1258 #. If an optimal bound was previously computed and 1259 cached, Sage returns it 1260 1261 #. Otherwise, Sage returns the reduced norm of ``self`` 1262 1263 As a consequence, the output of this function may depend 1264 on previous computations (an example is given below). 1265 1266 EXAMPLES:: 1267 1268 sage: k.<t> = GF(5^3) 1269 sage: Frob = k.frobenius_endomorphism() 1270 sage: S.<x> = k['x',Frob] 1271 sage: Z = S.center() 1272 1273 sage: a = x^2 + (4*t + 2)*x + 4*t^2 + 3 1274 sage: b = a.bound(); b 1275 (x^3)^2 + (x^3) + 4 1276 1277 Note that the parent of `b` is the center of the `S` 1278 (and not `S` itself):: 1279 1280 sage: b.parent() 1281 Center of Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5: 1282 Univariate Polynomial Ring in (x^3) over Finite Field of size 5 1283 sage: b.parent() == Z 1284 True 1285 1286 We check that `b` is divisible by `a`:: 1287 1288 sage: S(b).is_divisible_by(a) 1289 True 1290 sage: S(b).is_divisible_by(a,side=Left) 1291 True 1292 1293 Actually, `b` is the reduced norm of `a`:: 1294 1295 sage: b == a.reduced_norm() 1296 True 1297 1298 Now, we compute the optimal bound of `a` and see that 1299 it affects the behaviour of ``bound()``:: 1300 1301 sage: a.optimal_bound() 1302 (x^3) + 3 1303 sage: a.bound() 1304 (x^3) + 3 1305 1306 We finally check that if `a` is a central skew polynomial, 1307 then ``a.bound()`` returns simply `a`:: 1308 1309 sage: a = S(Z.random_element(degree=4)); a 1310 2*x^12 + x^9 + 2*x^3 1311 sage: b = a.bound(); b 1312 2*(x^3)^4 + (x^3)^3 + 2*(x^3) 1313 sage: a == b 1314 True 1315 """ 1316 center = self.parent().center() 1317 try: 1318 return center(self) 1319 except ValueError: 1320 pass 1321 if not self._optbound is None: 1322 return center(self._optbound) 1323 return self.reduced_norm() 1324 1325 1326 def optimal_bound(self): 1327 """ 1328 Return the optimal bound of this skew polynomial (i.e. 1329 the monic multiple of this skew polynomial of minimal 1330 degree lying in the center). 1331 1332 .. NOTE:: 1333 1334 The result is cached. 1335 1336 EXAMPLES:: 1337 1338 sage: k.<t> = GF(5^3) 1339 sage: Frob = k.frobenius_endomorphism() 1340 sage: S.<x> = k['x',Frob] 1341 sage: Z = S.center() 1342 1343 sage: a = x^2 + (4*t + 2)*x + 4*t^2 + 3 1344 sage: b = a.optimal_bound(); b 1345 (x^3) + 3 1346 1347 Note that the parent of `b` is the center of the `S` 1348 (and not `S` itself):: 1349 1350 sage: b.parent() 1351 Center of Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5: 1352 Univariate Polynomial Ring in (x^3) over Finite Field of size 5 1353 sage: b.parent() == Z 1354 True 1355 1356 We check that `b` is divisible by `a`:: 1357 1358 sage: S(b).is_divisible_by(a) 1359 True 1360 sage: S(b).is_divisible_by(a,side=Left) 1361 True 1362 """ 1363 center = self.parent().center() 1364 if self._optbound is None: 1365 try: 1366 self._optbound = center(self).monic() 1367 except ValueError: 1368 bound = self._matphir_c().minimal_polynomial() 1369 section = center._embed_basering.section() 1370 self._optbound = [ section(x) for x in bound.list() ] 1371 return center(self._optbound) 1372 1373 1374 def is_irreducible(self): 1375 """ 1376 Return true if this skew polynomial is irreducible. 1377 1378 EXAMPLES:: 1379 1380 sage: k.<t> = GF(5^3) 1381 sage: Frob = k.frobenius_endomorphism() 1382 sage: S.<x> = k['x',Frob] 1383 1384 sage: a = x^2 + t*x + 1 1385 sage: a.is_irreducible() 1386 False 1387 sage: a.factor() 1388 (x + 4*t^2 + 4*t + 1) * (x + 3*t + 2) 1389 1390 sage: a = x^2 + t*x + t + 1 1391 sage: a.is_irreducible() 1392 True 1393 sage: a.factor() 1394 x^2 + t*x + t + 1 1395 1396 Skew polynomials of degree `1` are of course irreducible:: 1397 1398 sage: a = x + t 1399 sage: a.is_irreducible() 1400 True 1401 1402 A random irreducible skew polynomial is irreducible:: 1403 1404 sage: a = S.random_irreducible(degree=4,monic=True); a # random 1405 x^4 + (t + 1)*x^3 + (3*t^2 + 2*t + 3)*x^2 + 3*t*x + 3*t 1406 sage: a.is_irreducible() 1407 True 1408 1409 By convention, constant skew polynomials are not irreducible:: 1410 1411 sage: S(1).is_irreducible() 1412 False 1413 sage: S(0).is_irreducible() 1414 False 1415 """ 1416 return self.reduced_norm().is_irreducible() 1417 1418 1419 def type(self,N): 1420 """ 1421 INPUT: 1422 1423 - ``N`` -- an irreducible polynomial in the 1424 center of the underlying skew polynomial ring 1425 1426 OUTPUT: 1427 1428 The `N`-type of this skew polynomial 1429 1430 .. NOTE:: 1431 1432 The result is cached. 1433 1434 DEFINITION: 1435 1436 The `N`-type of a skew polynomial `a` is the Partition 1437 `(t_0, t_1, t_2, ...)` defined by 1438 1439 .. MATH:: 1440 1441 t_0 + \cdots + t_i = \frac{\deg gcd(a,N^i)}{\deg N} 1442 1443 where `\deg N` is the degree of `N` considered as an 1444 element in the center. 1445 1446 This notion has an important mathematic interest because 1447 it corresponds to the Jordan type of the `N`-typical part 1448 of the associated Galois representation. 1449 1450 EXAMPLES:: 1451 1452 sage: k.<t> = GF(5^3) 1453 sage: Frob = k.frobenius_endomorphism() 1454 sage: S.<x> = k['x',Frob] 1455 sage: Z = S.center(); x3 = Z.gen() 1456 1457 sage: a = x^4 + x^3 + (4*t^2 + 4)*x^2 + (t^2 + 2)*x + 4*t^2 1458 sage: N = x3^2 + x3 + 1 1459 sage: a.type(N) 1460 [1] 1461 sage: N = x3 + 1 1462 sage: a.type(N) 1463 [2] 1464 1465 sage: a = x^3 + (3*t^2 + 1)*x^2 + (3*t^2 + t + 1)*x + t + 1 1466 sage: N = x3 + 1 1467 sage: a.type(N) 1468 [2, 1] 1469 1470 If `N` does not divide the reduced map of `a`, the type 1471 is empty:: 1472 1473 sage: N = x3 + 2 1474 sage: a.type(N) 1475 [] 1476 1477 If `a = N`, the type is just `[r]` where `r` is the order 1478 of the twist map ``Frob``:: 1479 1480 sage: N = x3^2 + x3 + 1 1481 sage: S(N).type(N) 1482 [3] 1483 1484 `N` must be irreducible:: 1485 1486 sage: N = (x3 + 1) * (x3 + 2) 1487 sage: a.type(N) 1488 Traceback (most recent call last): 1489 ... 1490 ValueError: N is not irreducible 1491 """ 1492 try: 1493 return self._types[N] 1494 except (KeyError, TypeError): 1495 if not N.is_irreducible(): 1496 raise ValueError("N is not irreducible") 1497 skew_ring = self._parent 1498 if self._norm_factor is None: 1499 m = -1 1500 else: 1501 i = [ n for n,_ in self._norm_factor ].index(N) 1502 m = self._norm_factor[i][1] 1503 NS = skew_ring(N) 1504 type = [ ] 1505 degN = N.degree() 1506 while True: 1507 d = self.gcd(NS) 1508 deg = d.degree()/degN 1509 if deg == 0: 1510 break 1511 if m >= 0: 1512 if deg == 1: 1513 type += m * [1] 1514 break 1515 m -= deg 1516 self = self // d 1517 type.append(deg) 1518 type = Partition(type) 1519 if self._types is None: 1520 self._types = { N: type } 1521 else: 1522 self._types[N] = type 1523 return type 1524 1525 1526 # Finding divisors 1527 # ---------------- 1528 1529 cdef SkewPolynomial_finite_field_dense _rdivisor_c(P, CenterSkewPolynomial_generic_dense N): 1530 """ 1531 cython procedure computing an irreducible monic right divisor 1532 of `P` whose reduced norm is `N` 1533 1534 .. WARNING:: 1535 1536 `N` needs to be an irreducible factor of the 1537 reduced norm of `P`. This function does not check 1538 this (and his behaviour is not defined if the 1539 require property doesn't hold). 1540 """ 1541 cdef skew_ring = P._parent 1542 cdef Py_ssize_t d = N.degree() 1543 cdef Py_ssize_t e = P.degree()/d 1544 cdef SkewPolynomial_finite_field_dense D 1545 if e == 1: 1546 D = <SkewPolynomial_finite_field_dense>P._new_c(list(P.__coeffs),skew_ring) 1547 D._inplace_rmonic() 1548 return D 1549 1550 E = N.parent().base_ring().extension(N,name='xr') 1551 PE = PolynomialRing(E,name='T') 1552 cdef Integer exp 1553 if skew_ring.characteristic() != 2: 1554 exp = Integer((E.cardinality()-1)/2) 1555 cdef SkewPolynomial_finite_field_dense NS = <SkewPolynomial_finite_field_dense>skew_ring(N) 1556 cdef SkewPolynomial_finite_field_dense Q = <SkewPolynomial_finite_field_dense>(NS // P) 1557 cdef SkewPolynomial_finite_field_dense R, X 1558 cdef Matrix_dense M = <Matrix_dense?>MatrixSpace(E,e,e)(0) 1559 cdef Matrix_dense V = <Matrix_dense?>MatrixSpace(E,e,1)(0) 1560 cdef Matrix_dense W 1561 cdef Py_ssize_t i, j, t, r = skew_ring._order 1562 cdef Polynomial dd, xx, yy, zz 1563 1564 while 1: 1565 R = <SkewPolynomial_finite_field_dense>skew_ring.random_element((e*r-1,e*r-1)) 1566 R._inplace_lmul(Q) 1567 X = <SkewPolynomial_finite_field_dense>Q._new_c(Q.__coeffs[:],Q._parent) 1568 for j from 0 <= j < e: 1569 for i from 0 <= i < e: 1570 M.set_unsafe(i, j, E([skew_ring._retraction(X[t*r+i]) for t in range(d)])) 1571 X._inplace_lmul(R) 1572 X._inplace_rrem(NS) 1573 for i from 0 <= i < e: 1574 V.set_unsafe(i, 0, E([skew_ring._retraction(X[t*r+i]) for t in range(d)])) 1575 W = M._solve_right_nonsingular_square(V) 1576 if M*W != V: 1577 skew_ring._new_retraction_map() 1578 continue 1579 xx = PE(W.list()+[E(-1)]) 1580 if skew_ring.characteristic() == 2: 1581 yy = PE.gen() 1582 zz = PE.gen() 1583 for i from 1 <= i < d: 1584 zz = (zz*zz) % xx 1585 yy += zz 1586 dd = xx.gcd(yy) 1587 if dd.degree() != 1: continue 1588 else: 1589 yy = PE.gen().__pow__(exp,xx) - 1 1590 dd = xx.gcd(yy) 1591 if dd.degree() != 1: 1592 yy += 2 1593 dd = xx.gcd(yy) 1594 if dd.degree() != 1: continue 1595 D = P._rgcd(R + skew_ring.center()((dd[0]/dd[1]).list())) 1596 if D.degree() == 0: 1597 continue 1598 D._inplace_rmonic() 1599 D._init_cache() 1600 return D 1601 1602 1603 def irreducible_divisor(self,side=Right,distribution=None): 1604 """ 1605 INPUT: 1606 1607 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1608 1609 - ``distribution`` -- None (default) or ``uniform`` 1610 1611 - None: no particular specification 1612 1613 - ``uniform``: the returned irreducible divisor is 1614 uniformly distributed 1615 1616 .. NOTE:: 1617 1618 ``uniform`` is a little bit slower. 1619 1620 OUTPUT: 1621 1622 - an irreducible monic ``side`` divisor of ``self`` 1623 1624 EXAMPLES:: 1625 1626 sage: k.<t> = GF(5^3) 1627 sage: Frob = k.frobenius_endomorphism() 1628 sage: S.<x> = k['x',Frob] 1629 sage: a = x^6 + 3*t*x^5 + (3*t + 1)*x^3 + (4*t^2 + 3*t + 4)*x^2 + (t^2 + 2)*x + 4*t^2 + 3*t + 3 1630 1631 sage: dr = a.irreducible_divisor(); dr # random 1632 x^3 + (2*t^2 + t + 4)*x^2 + (4*t + 1)*x + 4*t^2 + t + 1 1633 sage: a.is_divisible_by(dr) 1634 True 1635 1636 sage: dl = a.irreducible_divisor(side=Left); dl # random 1637 x^3 + (2*t^2 + t + 1)*x^2 + (4*t^2 + 3*t + 3)*x + 4*t^2 + 2*t + 1 1638 sage: a.is_divisible_by(dl,side=Left) 1639 True 1640 1641 Right divisors are cached. Hence, if we ask again for a 1642 right divisor, we will get the same answer:: 1643 1644 sage: a.irreducible_divisor() # random 1645 x^3 + (2*t^2 + t + 4)*x^2 + (4*t + 1)*x + 4*t^2 + t + 1 1646 1647 However the algorithm is probabilistic. Hence, if we first 1648 reinitialiaze `a`, we may get a different answer:: 1649 1650 sage: a = x^6 + 3*t*x^5 + (3*t + 1)*x^3 + (4*t^2 + 3*t + 4)*x^2 + (t^2 + 2)*x + 4*t^2 + 3*t + 3 1651 sage: a.irreducible_divisor() # random 1652 x^3 + (t^2 + 3*t + 4)*x^2 + (t + 2)*x + 4*t^2 + t + 1 1653 1654 We can also generate uniformly distributed irreducible monic 1655 divisors as follows:: 1656 1657 sage: a.irreducible_divisor(distribution="uniform") # random 1658 x^3 + (4*t + 2)*x^2 + (2*t^2 + 2*t + 2)*x + 2*t^2 + 2 1659 sage: a.irreducible_divisor(distribution="uniform") # random 1660 x^3 + (t^2 + 2)*x^2 + (3*t^2 + 1)*x + 4*t^2 + 2*t 1661 sage: a.irreducible_divisor(distribution="uniform") # random 1662 x^3 + x^2 + (4*t^2 + 2*t + 4)*x + t^2 + 3 1663 1664 By convention, the zero skew polynomial has no irreducible 1665 divisor: 1666 1667 sage: S(0).irreducible_divisor() 1668 Traceback (most recent call last): 1669 ... 1670 ValueError: 0 has no irreducible divisor 1671 """ 1672 if self.is_zero(): 1673 raise ValueError("0 has no irreducible divisor") 1674 if not (distribution is None or distribution == "uniform"): 1675 raise ValueError("distribution must be None or 'uniform'") 1676 if distribution == "uniform": 1677 skew_ring = self._parent 1678 center = skew_ring.center() 1679 cardcenter = center.base_ring().cardinality() 1680 gencenter = center.gen() 1681 count = [ ] 1682 total = 0 1683 F = self.reduced_norm_factor() 1684 for n,_ in F: 1685 if n == gencenter: 1686 total += 1 1687 else: 1688 degn = n.degree() 1689 P = self.gcd(skew_ring(n)) 1690 m = P.degree()/degn 1691 cardL = cardcenter**degn 1692 total += (cardL**m - 1)/(cardL - 1) 1693 count.append(total) 1694 if total == 0: 1695 raise ValueError("No irreducible divisor having given reduced norm") 1696 random = ZZ.random_element(total) 1697 for i in range(len(F)): 1698 if random < count[i]: 1699 N = F[i][0] 1700 break 1701 else: 1702 N = self.reduced_norm_factor()[0][0] 1703 return self.irreducible_divisor_with_norm(N,side=side,distribution=distribution) 1704 1705 1706 def irreducible_divisor_with_norm(self,N,side=Right,distribution=None): # Ajouter side 1707 """ 1708 INPUT: 1709 1710 - ``N`` -- an irreducible polynomial in the center 1711 of the underlying skew polynomial ring 1712 1713 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1714 1715 - ``distribution`` -- None (default) or ``uniform`` 1716 1717 - None: no particular specification 1718 1719 - ``uniform``: the returned irreducible divisor is 1720 uniformly distributed 1721 1722 .. NOTE:: 1723 1724 ``uniform`` is a little bit slower. 1725 1726 OUTPUT: 1727 1728 - an irreducible monic ``side`` divisor of ``self`` 1729 whose reduced norm is similar to `N` (i.e. `N` times 1730 a unit). 1731 1732 EXAMPLES:: 1733 1734 sage: k.<t> = GF(5^3) 1735 sage: Frob = k.frobenius_endomorphism() 1736 sage: S.<x> = k['x',Frob] 1737 sage: Z = S.center(); x3 = Z.gen() 1738 sage: a = x^6 + 3*x^3 + 2 1739 1740 sage: d1 = a.irreducible_divisor_with_norm(x3+1); d1 # random 1741 x + t^2 + 3*t 1742 sage: a.is_divisible_by(d1) 1743 True 1744 sage: d1.reduced_norm() 1745 (x^3) + 1 1746 1747 sage: d2 = a.irreducible_divisor_with_norm(x3+2); d2 # random 1748 x + 2*t^2 + 3*t + 2 1749 sage: a.is_divisible_by(d2) 1750 True 1751 sage: d2.reduced_norm() 1752 (x^3) + 2 1753 1754 sage: d3 = a.irreducible_divisor_with_norm(x3+3) 1755 Traceback (most recent call last): 1756 ... 1757 ValueError: No irreducible divisor having given reduced norm 1758 1759 We can also generate uniformly distributed irreducible monic 1760 divisors as follows:: 1761 1762 sage: a.irreducible_divisor_with_norm(x3+1,distribution="uniform") # random 1763 x + 3*t^2 + 3*t + 1 1764 sage: a.irreducible_divisor_with_norm(x3+1,distribution="uniform") # random 1765 x + 1 1766 sage: a.irreducible_divisor_with_norm(x3+1,distribution="uniform") # random 1767 x + 2*t^2 + 4*t 1768 """ 1769 cdef SkewPolynomial_finite_field_dense cP1 1770 cdef CenterSkewPolynomial_generic_dense cN 1771 if self.is_zero(): 1772 raise "No irreducible divisor having given reduced norm" 1773 skew_ring = self._parent 1774 center = skew_ring.center() 1775 try: 1776 N = center(N) 1777 except TypeError: 1778 raise TypeError("N must be a polynomial in the center") 1779 cardcenter = center.base_ring().cardinality() 1780 gencenter = center.gen() 1781 1782 if N == gencenter: 1783 if self[0] == 0: 1784 return skew_ring.gen() 1785 else: 1786 raise ValueError("No irreducible divisor having given reduced norm") 1787 1788 D = None 1789 try: 1790 D = self._rdivisors[N] 1791 except (KeyError, TypeError): 1792 if N.is_irreducible(): 1793 cP1 = <SkewPolynomial_finite_field_dense>self._rgcd(self._parent(N)) 1794 cN = <CenterSkewPolynomial_generic_dense>N 1795 if cP1.degree() > 0: 1796 D = cP1._rdivisor_c(cN) 1797 if self._rdivisors is None: 1798 self._rdivisors = { N: D } 1799 else: 1800 self._rdivisors[N] = D 1801 distribution = "" 1802 if D is None: 1803 raise ValueError("No irreducible divisor having given reduced norm") 1804 1805 NS = self._parent(N) 1806 degN = N.degree() 1807 if side is Right: 1808 if distribution == "uniform": 1809 P1 = self._rgcd(NS) 1810 if P1.degree() != degN: 1811 Q1 = NS // P1 1812 deg = P1.degree()-1 1813 while True: 1814 R = Q1*skew_ring.random_element((deg,deg)) 1815 if P1.gcd(R) == 1: 1816 break 1817 D = P1.gcd(D*R) 1818 return D 1819 else: 1820 deg = NS.degree()-1 1821 P1 = self.lgcd(NS) 1822 while True: 1823 if distribution == "uniform": 1824 while True: 1825 R = skew_ring.random_element((deg,deg)) 1826 if NS.gcd(R) == 1: 1827 break 1828 D = NS.gcd(D*R) 1829 Dp = NS // D 1830 LDp = P1.gcd(Dp) 1831 LD = P1 // LDp 1832 if LD.degree() == degN: 1833 return LD 1834 distribution = "uniform" 1835 1836 1837 def irreducible_divisors(self,side=Right): 1838 """ 1839 INPUT: 1840 1841 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1842 1843 OUTPUT: 1844 1845 An iterator over all irreducible monic ``side`` divisors 1846 of this skew polynomial 1847 1848 EXAMPLES: 1849 1850 sage: k.<t> = GF(5^3) 1851 sage: Frob = k.frobenius_endomorphism() 1852 sage: S.<x> = k['x',Frob] 1853 sage: a = x^4 + 2*t*x^3 + 3*t^2*x^2 + (t^2 + t + 1)*x + 4*t + 3 1854 sage: iter = a.irreducible_divisors(); iter 1855 <generator object at 0x...> 1856 sage: iter.next() # random 1857 x + 2*t^2 + 4*t + 4 1858 sage: iter.next() # random 1859 x + 3*t^2 + 4*t + 1 1860 1861 We can use this function to build the list of all monic 1862 irreducible divisors of `a`:: 1863 1864 sage: rightdiv = [ d for d in a.irreducible_divisors() ] 1865 sage: leftdiv = [ d for d in a.irreducible_divisors(side=Left) ] 1866 1867 We do some checks:: 1868 1869 sage: len(rightdiv) == a.count_irreducible_divisors() 1870 True 1871 sage: len(rightdiv) == len(Set(rightdiv)) # check no duplicates 1872 True 1873 sage: for d in rightdiv: 1874 ... if not a.is_divisible_by(d): 1875 ... print "Found %s which is not a right divisor" % d 1876 ... elif not d.is_irreducible(): 1877 ... print "Found %s which is not irreducible" % d 1878 1879 sage: len(leftdiv) == a.count_irreducible_divisors(side=Left) 1880 True 1881 sage: len(leftdiv) == len(Set(leftdiv)) # check no duplicates 1882 True 1883 sage: for d in leftdiv: 1884 ... if not a.is_divisible_by(d,side=Left): 1885 ... print "Found %s which is not a left divisor" % d 1886 ... elif not d.is_irreducible(): 1887 ... print "Found %s which is not irreducible" % d 1888 1889 Note that left divisors and right divisors differ:: 1890 1891 sage: Set(rightdiv) == Set(leftdiv) 1892 False 1893 1894 Note that the algorithm is probabilistic. As a consequence, if we 1895 build again the list of right monic irreducible divisors of `a`, we 1896 may get a different ordering:: 1897 1898 sage: rightdiv2 = [ d for d in a.irreducible_divisors() ] 1899 sage: rightdiv == rightdiv2 1900 False 1901 sage: Set(rightdiv) == Set(rightdiv2) 1902 True 1903 """ 1904 return self._irreducible_divisors(side) 1905 1906 1907 def _irreducible_divisors(self,side): # prendre side en compte 1908 """ 1909 Return an iterator over all irreducible monic 1910 divisors of this skew polynomial. 1911 1912 Do not use this function. Use instead 1913 ``self.irreducible_divisors()``. 1914 """ 1915 if self.is_zero(): 1916 return 1917 skew_ring = self._parent 1918 center = skew_ring.center() 1919 kfixed = center.base_ring() 1920 F = self.reduced_norm_factor() 1921 oppside = side.opposite() 1922 for N,_ in F: 1923 if N == center.gen(): 1924 yield skew_ring.gen() 1925 continue 1926 degN = N.degree() 1927 NS = skew_ring(N) 1928 P = self.gcd(NS,side=side) 1929 m = P.degree()/degN 1930 if m == 1: 1931 yield P 1932 continue 1933 degrandom = P.degree() - 1 1934 Q,_ = NS.quo_rem(P,side=side) 1935 P1 = self.irreducible_divisor_with_norm(N,side=side) 1936 Q1,_ = P.quo_rem(P1,side=side) 1937 while True: 1938 R = skew_ring.random_element((degrandom,degrandom)) 1939 if side is Right: 1940 g = (R*Q).rem(P,side=Left) 1941 else: 1942 g = (Q*R).rem(P) 1943 if g.gcd(P,side=oppside) != 1: continue 1944 L = Q1 1945 V = L 1946 for i in range(1,m): 1947 if side is Right: 1948 L = (g*L).gcd(P,side=Left) 1949 else: 1950 L = (L*g).gcd(P) 1951 V = V.gcd(L,side=oppside) 1952 if V == 1: break 1953 rng = xmrange_iter([kfixed]*degN,center) 1954 for i in range(m): 1955 for pol in xmrange_iter([rng]*i): 1956 f = skew_ring(1) 1957 for j in range(i): 1958 coeff = pol.pop() 1959 f = (g*f+coeff).rem(P,side=oppside) 1960 if side is Right: 1961 d = (f*Q1).gcd(P,side=Left) 1962 else: 1963 d = (Q1*f).gcd(P) 1964 d,_ = P.quo_rem(d,side=oppside) 1965 yield d 1966 1967 1968 def count_irreducible_divisors(self,side=Right): 1969 """ 1970 INPUT: 1971 1972 - ``side`` -- ``Left`` or ``Right`` (default: Right) 1973 1974 OUTPUT: 1975 1976 The number of irreducible monic ``side`` divisors of 1977 this skew polynomial. 1978 1979 .. NOTE:: 1980 1981 Actually, one can prove that there are always as 1982 many left irreducible monic divisors as right 1983 irreducible monic divisors. 1984 1985 EXAMPLES:: 1986 1987 sage: k.<t> = GF(5^3) 1988 sage: Frob = k.frobenius_endomorphism() 1989 sage: S.<x> = k['x',Frob] 1990 1991 We illustrate that a skew polynomial may have a number of irreducible 1992 divisors greater than its degree. 1993 1994 sage: a = x^4 + (4*t + 3)*x^3 + t^2*x^2 + (4*t^2 + 3*t)*x + 3*t 1995 sage: a.count_irreducible_divisors() 1996 12 1997 sage: a.count_irreducible_divisors(side=Left) 1998 12 1999 2000 We illustrate that an irreducible polynomial in the center have 2001 in general a lot of irreducible divisors in the skew polynomial 2002 ring:: 2003 2004 sage: Z = S.center(); x3 = Z.gen() 2005 sage: N = x3^5 + 4*x3^4 + 4*x3^2 + 4*x3 + 3; N 2006 (x^3)^5 + 4*(x^3)^4 + 4*(x^3)^2 + 4*(x^3) + 3 2007 sage: N.is_irreducible() 2008 True 2009 sage: S(N).count_irreducible_divisors() 2010 9768751 2011 """ 2012 if self.is_zero(): 2013 return 0 2014 skew_ring = self.parent() 2015 cardcenter = skew_ring.center().base_ring().cardinality() 2016 gencenter = skew_ring.center().gen() 2017 F = self.reduced_norm_factor() 2018 val = self.valuation() 2019 self >>= val 2020 count = 0 2021 if val > 0: 2022 count = 1 2023 for N,_ in F: 2024 if N == gencenter: 2025 continue 2026 degN = N.degree() 2027 P = self.gcd(skew_ring(N), side=side) 2028 m = P.degree()/degN 2029 cardL = cardcenter**degN 2030 count += (cardL**m - 1)/(cardL - 1) 2031 return count 2032 2033 2034 # Finding factorizations 2035 # ---------------------- 2036 2037 cdef _factor_c(self): 2038 """ 2039 Compute a factorization of ``self`` 2040 """ 2041 cdef skew_ring = self._parent 2042 cdef Py_ssize_t degQ, degrandom, m, mP, i 2043 cdef CenterSkewPolynomial_generic_dense N 2044 cdef SkewPolynomial_finite_field_dense poly = <SkewPolynomial_finite_field_dense>self.rmonic() 2045 cdef val = poly._val_inplace_unit() 2046 if val == -1: 2047 return Factorization([], sort=False, unit=skew_ring.zero_element()) 2048 cdef list factors = [ (skew_ring.gen(), val) ] 2049 cdef SkewPolynomial_finite_field_dense P, Q, P1, NS, g, right, Pn 2050 cdef SkewPolynomial_finite_field_dense right2 = skew_ring(1) << val 2051 cdef RingElement unit = <RingElement>self.leading_coefficient() 2052 cdef Polynomial gencenter = skew_ring.center().gen() 2053 cdef Py_ssize_t p = skew_ring.characteristic() 2054 cdef F = self.reduced_norm_factor() 2055 2056 for N,m in F: 2057 if N == gencenter: 2058 continue 2059 degN = N.degree() 2060 if poly.degree() == degN: 2061 factors.append((poly,1)) 2062 break 2063 NS = <SkewPolynomial_finite_field_dense>skew_ring(N) 2064 P1 = None 2065 while 1: 2066 P = <SkewPolynomial_finite_field_dense>poly._rgcd(NS) 2067 P._inplace_rmonic() 2068 mP = P.degree() / degN 2069 if mP == 0: break 2070 if mP == 1: 2071 factors.append((P,1)) 2072 poly._inplace_rfloordiv(P) 2073 for i from 1 <= i < m: 2074 if poly.degree() == degN: 2075 factors.append((poly,1)) 2076 break 2077 P = poly._rgcd(NS) 2078 P._inplace_rmonic() 2079 factors.append((P,1)) 2080 poly._inplace_rfloordiv(P) 2081 break 2082 if P1 is None: 2083 P1 = P._rdivisor_c(N) 2084 Q = <SkewPolynomial_finite_field_dense>NS._new_c(NS.__coeffs[:], NS._parent) 2085 Q._inplace_rfloordiv(P) 2086 Q._inplace_lmul(P1) 2087 factors.append((P1,1)) 2088 right = <SkewPolynomial_finite_field_dense>P1._new_c(P1.__coeffs[:], P1._parent) 2089 m -= (mP-1) 2090 degrandom = P.degree() 2091 while mP > 2: 2092 while 1: 2093 g = <SkewPolynomial_finite_field_dense>skew_ring.random_element((degrandom,degrandom)) 2094 g._inplace_lmul(Q) 2095 g._inplace_rgcd(P) 2096 Pn = right._coeff_llcm(g) 2097 if len(Pn.__coeffs)-1 == degN: break 2098 Pn._inplace_rmonic() 2099 factors.append((Pn,1)) 2100 right._inplace_lmul(Pn) 2101 degrandom -= degN 2102 mP -= 1 2103 poly._inplace_rfloordiv(right) 2104 P1,_ = P.rquo_rem(right) 2105 factors.reverse() 2106 return Factorization(factors, sort=False, unit=unit) 2107 2108 2109 cdef _factor_uniform_c(self): 2110 """ 2111 Compute a uniformly distrbuted factorization of ``self`` 2112 """ 2113 skew_ring = self._parent 2114 cdef Integer cardE, cardcenter = skew_ring.center().base_ring().cardinality() 2115 cdef CenterSkewPolynomial_generic_dense gencenter = <CenterSkewPolynomial_generic_dense>skew_ring.center().gen() 2116 cdef SkewPolynomial_finite_field_dense gen = <SkewPolynomial_finite_field_dense>skew_ring.gen() 2117 2118 cdef list factorsN = [ ] 2119 cdef dict dict_divisor = { } 2120 cdef dict dict_type = { } 2121 cdef dict dict_right = { } 2122 cdef CenterSkewPolynomial_generic_dense N 2123 cdef Py_ssize_t m 2124 cdef list type 2125 2126 for N,m in self.reduced_norm_factor(): 2127 factorsN += m * [N] 2128 if N == gencenter: continue 2129 type = list(self.type(N)) 2130 dict_type[N] = type 2131 if type[0] > 1: 2132 dict_divisor[N] = self.irreducible_divisor_with_norm(N) 2133 dict_right[N] = skew_ring(1) 2134 cdef list indices = list(Permutations(len(factorsN)).random_element()) 2135 2136 cdef RingElement unit = self.leading_coefficient() 2137 cdef SkewPolynomial_finite_field_dense left = self._new_c(self.__coeffs[:],skew_ring) 2138 left._inplace_rmonic() 2139 cdef SkewPolynomial_finite_field_dense right = <SkewPolynomial_finite_field_dense>skew_ring(1) 2140 cdef SkewPolynomial_finite_field_dense L, R 2141 cdef SkewPolynomial_finite_field_dense NS, P, Q, D, D1, D2, d 2142 cdef list factors = [ ] 2143 cdef list maxtype 2144 cdef Py_ssize_t i, j, degN, deg 2145 cdef count, maxcount 2146 2147 for i in indices: 2148 N = factorsN[i-1] 2149 if N == gencenter: 2150 D1 = gen 2151 else: 2152 type = dict_type[N] 2153 NS = skew_ring(N) 2154 P = left.gcd(NS) 2155 if type[0] == 1: 2156 D1 = P 2157 else: 2158 R = right._new_c(right.__coeffs[:],skew_ring) 2159 R._inplace_rfloordiv(dict_right[N]) 2160 D = R._coeff_llcm(dict_divisor[N]) 2161 maxtype = list(type) 2162 maxtype[-1] -= 1 2163 degN = N.degree() 2164 cardE = cardcenter ** degN 2165 maxcount = q_jordan(Partition(maxtype),cardE) 2166 Q = NS // P 2167 deg = P.degree()-1 2168 while 1: 2169 while 1: 2170 R = <SkewPolynomial_finite_field_dense>skew_ring.random_element((deg,deg)) 2171 R._inplace_lmul(Q) 2172 if P._rgcd(R).degree() == 0: 2173 break 2174 D1 = P._rgcd(D*R) 2175 D1._inplace_rmonic() 2176 2177 L = left._new_c(list(left.__coeffs),skew_ring) 2178 L._inplace_rfloordiv(D1) 2179 degN = N.degree() 2180 for j in range(len(type)): 2181 if type[j] == 1: 2182 newtype = type[:-1] 2183 break 2184 d = L._rgcd(NS) 2185 d._inplace_rmonic() 2186 deg = d.degree() / degN 2187 if deg < type[j]: 2188 newtype = type[:] 2189 newtype[j] = deg 2190 break 2191 L._inplace_rfloordiv(d) 2192 count = q_jordan(Partition(newtype),cardE) 2193 if ZZ.random_element(maxcount) < count: 2194 break 2195 dict_type[N] = newtype 2196 2197 D2 = D._new_c(list(D.__coeffs),skew_ring) 2198 D2._inplace_rmonic() 2199 while D2 == D1: 2200 while 1: 2201 R = <SkewPolynomial_finite_field_dense>skew_ring.random_element((deg,deg)) 2202 R._inplace_lmul(Q) 2203 if P._rgcd(R).degree() == 0: 2204 break 2205 D2 = P._rgcd(D*R) 2206 D2._inplace_rmonic() 2207 dict_divisor[N] = D1._coeff_llcm(D2) 2208 factors.append((D1,1)) 2209 left._inplace_rfloordiv(D1) 2210 right._inplace_lmul(D1) 2211 dict_right[N] = right._new_c(list(right.__coeffs),skew_ring) 2212 2213 factors.reverse() 2214 return Factorization(factors,sort=False,unit=unit) 2215 2216 2217 def factor(self,distribution=None): 2218 """ 2219 Return a factorization of this skew polynomial. 2220 2221 INPUT: 2222 2223 - ``distribution`` -- None (default) or ``uniform`` 2224 2225 - None: no particular specification 2226 2227 - ``uniform``: the returned factorization is uniformly 2228 distributed among all possible factorizations 2229 2230 .. NOTE:: 2231 2232 ``uniform`` is a little bit slower. 2233 2234 OUTPUT: 2235 2236 - ``Factorization`` -- a factorization of self as a 2237 product of a unit and a product of irreducible monic 2238 factors 2239 2240 EXAMPLES:: 2241 2242 sage: k.<t> = GF(5^3) 2243 sage: Frob = k.frobenius_endomorphism() 2244 sage: S.<x> = k['x',Frob] 2245 sage: a = x^3 + (t^2 + 4*t + 2)*x^2 + (3*t + 3)*x + t^2 + 1 2246 sage: F = a.factor(); F # random 2247 (x + t^2 + 4) * (x + t + 3) * (x + t) 2248 sage: F.value() == a 2249 True 2250 2251 The result of the factorization is cached. Hence, if we try 2252 again to factor `a`, we will get the same answer:: 2253 2254 sage: a.factor() # random 2255 (x + t^2 + 4) * (x + t + 3) * (x + t) 2256 2257 However, the algorithm is probabilistic. Hence if we first 2258 reinitialiaze `a`, we may get a different answer:: 2259 2260 sage: a = x^3 + (t^2 + 4*t + 2)*x^2 + (3*t + 3)*x + t^2 + 1 2261 sage: F = a.factor(); F # random 2262 (x + t^2 + t + 2) * (x + 2*t^2 + t + 4) * (x + t) 2263 sage: F.value() == a 2264 True 2265 2266 There is no guarantee on the distribution of the factorizations 2267 we get that way. (For this particular `a` for example, we get the 2268 uniform distribution on the subset of all factorizations ending 2269 by the factor `x + t`.) 2270 2271 If we rather want uniform distribution among all factorizations, 2272 we need to specify it as follows:: 2273 2274 sage: a.factor(distribution="uniform") # random 2275 (x + t^2 + 4) * (x + t) * (x + t + 3) 2276 sage: a.factor(distribution="uniform") # random 2277 (x + 2*t^2) * (x + t^2 + t + 1) * (x + t^2 + t + 2) 2278 sage: a.factor(distribution="uniform") # random 2279 (x + 2*t^2 + 3*t) * (x + 4*t + 2) * (x + 2*t + 2) 2280 2281 By convention, the zero skew polynomial has no factorization: 2282 2283 sage: S(0).factor() 2284 Traceback (most recent call last): 2285 ... 2286 ValueError: factorization of 0 not defined 2287 """ 2288 if not (distribution is None or distribution == "uniform"): 2289 raise ValueError("distribution must be None or 'uniform'") 2290 if self.is_zero(): 2291 raise ValueError("factorization of 0 not defined") 2292 sig_on() 2293 if distribution is None: 2294 if self._factorization is None: 2295 self._factorization = self._factor_c() 2296 F = self._factorization 2297 else: 2298 F = self._factor_uniform_c() 2299 if self._factorization is None: 2300 self._factorization = F 2301 sig_off() 2302 return F 2303 2304 2305 def count_factorizations(self): 2306 """ 2307 Return the number of factorizations (as a product of a 2308 unit and a product of irreducible monic factors) of this 2309 skew polynomial. 2310 2311 EXAMPLES:: 2312 2313 sage: k.<t> = GF(5^3) 2314 sage: Frob = k.frobenius_endomorphism() 2315 sage: S.<x> = k['x',Frob] 2316 sage: a = x^4 + (4*t + 3)*x^3 + t^2*x^2 + (4*t^2 + 3*t)*x + 3*t 2317 sage: a.count_factorizations() 2318 216 2319 2320 We illustrate that an irreducible polynomial in the center have 2321 in general a lot of distinct factorizations in the skew polynomial 2322 ring:: 2323 2324 sage: Z = S.center(); x3 = Z.gen() 2325 sage: N = x3^5 + 4*x3^4 + 4*x3^2 + 4*x3 + 3; N 2326 (x^3)^5 + 4*(x^3)^4 + 4*(x^3)^2 + 4*(x^3) + 3 2327 sage: N.is_irreducible() 2328 True 2329 sage: S(N).count_factorizations() 2330 30537115626 2331 """ 2332 if self.is_zero(): 2333 raise ValueError("factorization of 0 not defined") 2334 cardcenter = self._parent.center().base_ring().cardinality() 2335 gencenter = self._parent.center().gen() 2336 F = self.reduced_norm_factor() 2337 summ = 0 2338 count = 1 2339 for N,m in F: 2340 summ += m 2341 if m == 1: continue 2342 if N != gencenter: 2343 count *= q_jordan(self.type(N),cardcenter**N.degree()) 2344 count /= factorial(m) 2345 return count * factorial(summ) 2346 2347 def count_factorisations(self): 2348 """ 2349 Return the number of factorisations (as a product of a 2350 unit and a product of irreducible monic factors) of this 2351 skew polynomial. 2352 2353 EXAMPLES:: 2354 2355 sage: k.<t> = GF(5^3) 2356 sage: Frob = k.frobenius_endomorphism() 2357 sage: S.<x> = k['x',Frob] 2358 sage: a = x^4 + (4*t + 3)*x^3 + t^2*x^2 + (4*t^2 + 3*t)*x + 3*t 2359 sage: a.count_factorisations() 2360 216 2361 2362 We illustrate that an irreducible polynomial in the center have 2363 in general a lot of distinct factorisations in the skew polynomial 2364 ring:: 2365 2366 sage: Z = S.center(); x3 = Z.gen() 2367 sage: N = x3^5 + 4*x3^4 + 4*x3^2 + 4*x3 + 3; N 2368 (x^3)^5 + 4*(x^3)^4 + 4*(x^3)^2 + 4*(x^3) + 3 2369 sage: N.is_irreducible() 2370 True 2371 sage: S(N).count_factorisations() 2372 30537115626 2373 """ 2374 return self.count_factorizations() 2375 2376 2377 # Not optimized (many calls to reduced_norm, reduced_norm_factor,_rdivisor_c, which are slow) 2378 def factorizations(self): 2379 """ 2380 Return an iterator over all factorizations (as a product 2381 of a unit and a product of irreducible monic factors) of 2382 this skew polynomial. 2383 2384 .. NOTE:: 2385 2386 The algorithm is probabilistic. As a consequence, if 2387 we execute two times with the same input we can get 2388 the list of all factorizations in two differents orders. 2389 2390 EXAMPLES:: 2391 2392 sage: k.<t> = GF(5^3) 2393 sage: Frob = k.frobenius_endomorphism() 2394 sage: S.<x> = k['x',Frob] 2395 sage: a = x^3 + (t^2 + 1)*x^2 + (2*t + 3)*x + t^2 + t + 2 2396 sage: iter = a.factorizations(); iter 2397 <generator object at 0x...> 2398 sage: iter.next() # random 2399 (x + 3*t^2 + 4*t) * (x + 2*t^2) * (x + 4*t^2 + 4*t + 2) 2400 sage: iter.next() # random 2401 (x + 3*t^2 + 4*t) * (x + 3*t^2 + 2*t + 2) * (x + 4*t^2 + t + 2) 2402 2403 We can use this function to build the list of factorizations 2404 of `a`:: 2405 2406 sage: factorizations = [ F for F in a.factorizations() ] 2407 2408 We do some checks:: 2409 2410 sage: len(factorizations) == a.count_factorizations() 2411 True 2412 sage: len(factorizations) == len(Set(factorizations)) # check no duplicates 2413 True 2414 sage: for F in factorizations: 2415 ... if F.value() != a: 2416 ... print "Found %s which is not a correct factorization" % d 2417 ... continue 2418 ... for d,_ in F: 2419 ... if not d.is_irreducible(): 2420 ... print "Found %s which is not a correct factorization" % d 2421 """ 2422 if self.is_zero(): 2423 raise ValueError("factorization of 0 not defined") 2424 unit = self.leading_coefficient() 2425 poly = self.rmonic() 2426 for factors in self._factorizations_rec(): 2427 yield Factorization(factors,sort=False,unit=unit) 2428 def _factorizations_rec(self): 2429 if self.is_irreducible(): 2430 yield [ (self,1) ] 2431 else: 2432 for div in self._irreducible_divisors(Right): 2433 poly = self // div 2434 # Here, we should update poly._norm, poly._norm_factor, poly._rdivisors 2435 for factors in poly._factorizations_rec(): 2436 factors.append((div,1)) 2437 yield factors 2438 2439 2440 def factorisations(self): 2441 """ 2442 Return an iterator over all factorisations (as a product 2443 of a unit and a product of irreducible monic factors) of 2444 this skew polynomial. 2445 2446 .. NOTE:: 2447 2448 The algorithm is probabilistic. As a consequence, if 2449 we execute two times with the same input we can get 2450 the list of all factorizations in two differents orders. 2451 2452 EXAMPLES:: 2453 2454 sage: k.<t> = GF(5^3) 2455 sage: Frob = k.frobenius_endomorphism() 2456 sage: S.<x> = k['x',Frob] 2457 sage: a = x^3 + (t^2 + 1)*x^2 + (2*t + 3)*x + t^2 + t + 2 2458 sage: iter = a.factorisations(); iter 2459 <generator object at 0x...> 2460 sage: iter.next() # random 2461 (x + 3*t^2 + 4*t) * (x + 2*t^2) * (x + 4*t^2 + 4*t + 2) 2462 sage: iter.next() # random 2463 (x + 3*t^2 + 4*t) * (x + 3*t^2 + 2*t + 2) * (x + 4*t^2 + t + 2) 2464 2465 We can use this function to build the list of factorizations 2466 of `a`:: 2467 2468 sage: factorisations = [ F for F in a.factorisations() ] 2469 2470 We do some checks:: 2471 2472 sage: len(factorisations) == a.count_factorisations() 2473 True 2474 sage: len(factorisations) == len(Set(factorisations)) # check no duplicates 2475 True 2476 sage: for F in factorisations: 2477 ... if F.value() != a: 2478 ... print "Found %s which is not a correct factorization" % d 2479 ... continue 2480 ... for d,_ in F: 2481 ... if not d.is_irreducible(): 2482 ... print "Found %s which is not a correct factorization" % d 2483 """ 2484 return self.factorizations() 2485 2486 2487 # Karatsuba class 2488 ################# 2489 2490 cdef class SkewPolynomial_finite_field_karatsuba: 2491 """ 2492 A special class implementing Karatsuba multiplication and 2493 euclidean division on skew polynomial rings over finite fields. 2494 2495 Only for internal use! 2496 2497 TESTS:: 2498 2499 sage: k.<t> = GF(5^3) 2500 sage: Frob = k.frobenius_endomorphism() 2501 sage: S.<x> = k['x',Frob] 2502 2503 An instance of this class is automatically created when the skew 2504 ring is created. It is accessible via the key work ``_karatsuba_class``:: 2505 2506 sage: KarClass = S._karatsuba_class; KarClass 2507 Special class for Karatsuba multiplication and euclidean division on 2508 Skew Polynomial Ring in x over Finite Field in t of size 5^3 twisted by t |--> t^5 2509 2510 Each instance of this class is equipped with a ``cutoff`` variable. The 2511 default value is the maximum between 150 and the order of the twist map 2512 acting on the finite field:: 2513 2514 sage: KarClass.get_cutoff() 2515 150 2516 2517 We can set a new value to this variable as follows:: 2518 2519 sage: KarClass.set_cutoff(27) 2520 sage: KarClass.get_cutoff() 2521 27 2522 """ 2523 def __init__(self,parent,cutoff=0): 2524 """ 2525 Initialize a new instance of this class. 2526 """ 2527 self._parent = parent 2528 self._order = parent._order 2529 self._zero = parent.base_ring()(0) 2530 self.set_cutoff(cutoff) 2531 self._algo_matrix = 0 2532 self._t = None 2533 self._T = None 2534 self._Tinv = None 2535 2536 2537 def __repr__(self): 2538 """ 2539 Return a representation of ``self`` 2540 """ 2541 return "Special class for Karatsuba multiplication and euclidean division on\n%s" % self._parent 2542 2543 2544 def set_cutoff(self,cutoff=0): 2545 """ 2546 Set a new cutoff for all Karatsuba multiplications 2547 and euclidean divisions performed by this class. 2548 2549 INPUT: 2550 2551 - ``cutoff`` -- a nonnegative integer or +Infinity 2552 (default: 0) 2553 2554 .. NOTE:: 2555 2556 If ``cutoff`` is `0`, the new cutoff is set to the 2557 default value which is the maximum between 150 and 2558 the order of the twist map acting on the underlying 2559 finite field. 2560 2561 EXAMPLES:: 2562 2563 sage: k.<t> = GF(5^3) 2564 sage: Frob = k.frobenius_endomorphism() 2565 sage: S.<x> = k['x',Frob] 2566 sage: KarClass = S._karatsuba_class 2567 2568 sage: KarClass.set_cutoff(27) 2569 sage: KarClass.get_cutoff() 2570 27 2571 2572 sage: KarClass.set_cutoff(Infinity) 2573 sage: KarClass.get_cutoff() 2574 +Infinity 2575 2576 sage: KarClass.set_cutoff(0) # set the default value 2577 sage: KarClass.get_cutoff() 2578 150 2579 2580 The cutoff can't be less than the order of the twist map:: 2581 2582 sage: KarClass.set_cutoff(2) 2583 Traceback (most recent call last): 2584 ... 2585 ValueError: cutoff must be 0 or >= 3 2586 """ 2587 if cutoff == 0: 2588 self._cutoff = max(150,self._order) 2589 else: 2590 if cutoff < self._order: 2591 raise ValueError("cutoff must be 0 or >= %s" % self._order) 2592 if cutoff == Infinity: 2593 self._cutoff = -1 2594 else: 2595 self._cutoff = cutoff 2596 2597 2598 def get_cutoff(self): 2599 """ 2600 Return the current cutoff 2601 2602 EXAMPLES:: 2603 2604 sage: k.<t> = GF(5^3) 2605 sage: Frob = k.frobenius_endomorphism() 2606 sage: S.<x> = k['x',Frob] 2607 sage: KarClass = S._karatsuba_class 2608 sage: KarClass.get_cutoff() 2609 150 2610 sage: KarClass.set_cutoff(27) 2611 sage: KarClass.get_cutoff() 2612 27 2613 """ 2614 if self._cutoff < 0: 2615 return Infinity 2616 else: 2617 return self._cutoff 2618 2619 2620 def mul(self,left,right): 2621 """ 2622 INPUT: 2623 2624 - ``left`` -- a skew polynomial in the skew polynomial 2625 right attached to the class 2626 2627 - ``right`` -- a skew polynomial in the skew polynomial 2628 right attached to the class 2629 2630 OUTPUT: 2631 2632 The product left * right (computed by Karatsuba's algorithm) 2633 2634 EXAMPLES:: 2635 2636 sage: k.<t> = GF(5^3) 2637 sage: Frob = k.frobenius_endomorphism() 2638 sage: S.<x> = k['x',Frob] 2639 sage: KarClass = S._karatsuba_class 2640 2641 sage: a = S.random_element(degree=500) 2642 sage: b = S.random_element(degree=500) 2643 sage: c = KarClass.mul(a,b) 2644 sage: c == a*b 2645 True 2646 2647 .. NOTE:: 2648 2649 Behind the scene, the operator ``*`` calls actually 2650 the function ``KarClass.mul()``. 2651 """ 2652 cdef Py_ssize_t dx = left.degree(), dy = right.degree() 2653 cdef list x = left.list() 2654 cdef list y = right.list() 2655 cdef list res 2656 if self._cutoff < 0: 2657 res = self.mul_step(x,y) 2658 else: 2659 if dx < dy: 2660 res = self.mul_iter(y,x,1) 2661 else: 2662 res = self.mul_iter(x,y,0) 2663 return self._parent(res) 2664 2665 2666 def mul_matrix(self,left,right): 2667 """ 2668 INPUT: 2669 2670 - ``left`` -- an element in the skew polynomial ring 2671 attached to the class 2672 2673 - ``right`` -- an element in the skew polynomial ring 2674 attached to the class 2675 2676 OUTPUT: 2677 2678 The product left * right (computed by Karatsuba-matrix 2679 algorithm) 2680 2681 EXAMPLES:: 2682 2683 sage: k.<t> = GF(5^3) 2684 sage: Frob = k.frobenius_endomorphism() 2685 sage: S.<x> = k['x',Frob] 2686 sage: KarClass = S._karatsuba_class 2687 2688 sage: a = S.random_element(degree=500) 2689 sage: b = S.random_element(degree=500) 2690 sage: c = KarClass.mul_matrix(a,b) 2691 sage: c == a*b 2692 True 2693 """ 2694 cdef Py_ssize_t dx = left.degree(), dy = right.degree() 2695 cdef list x = left.list() 2696 cdef list y = right.list() 2697 cdef list res 2698 cdef Py_ssize_t save_cutoff = self._cutoff 2699 cdef Py_ssize_t i, j 2700 self._cutoff = int(self._order * self._order / 2) 2701 self._algo_matrix = 1 2702 if self._t is None: 2703 self._t = self._parent.base_ring().gen() 2704 if self._T is None: 2705 self._T = <Matrix_dense>zero_matrix(self._parent.base_ring(),self._order) 2706 for i from 0 <= i < self._order: 2707 map = self._parent.twist_map(-i) 2708 for j from 0 <= j < self._order: 2709 self._T.set_unsafe(i,j,map(self._t**j)) 2710 if self._Tinv is None: 2711 self._Tinv = <Matrix_dense>self._T.inverse() 2712 if dx < dy: 2713 res = self.mul_iter(y,x,1) 2714 else: 2715 res = self.mul_iter(x,y,0) 2716 self._cutoff = save_cutoff 2717 self._algo_matrix = 0 2718 return self._parent(res) 2719 2720 2721 def mul_list(self,x,y): 2722 """ 2723 INPUT: 2724 2725 - ``x`` -- a list of coefficients 2726 2727 - ``y`` -- a list of coefficients 2728 2729 OUTPUT: 2730 2731 The list of coefficients of the product `a * b` 2732 where `a` (resp. `b`) is the skew polynomial whose 2733 coefficients are given by `x` (resp. `y`). 2734 2735 EXAMPLES:: 2736 2737 sage: k.<t> = GF(5^3) 2738 sage: Frob = k.frobenius_endomorphism() 2739 sage: S.<u> = k['u',Frob] 2740 sage: KarClass = S._karatsuba_class 2741 2742 sage: x = [ k.random_element() for _ in range(500) ] 2743 sage: y = [ k.random_element() for _ in range(500) ] 2744 sage: z = KarClass.mul_list(x,y) 2745 sage: S(z) == S(x)*S(y) 2746 True 2747 """ 2748 if self._cutoff < 0: 2749 return self.mul_step(x,y) 2750 cdef Py_ssize_t dx = len(x)-1, dy = len(y)-1 2751 if dx < dy: 2752 return self.mul_iter(y,x,1) 2753 else: 2754 return self.mul_iter(x,y,0) 2755 2756 2757 cdef list mul_step(self, list x, list y): 2758 """ 2759 Multiplication step in Karatsuba algorithm 2760 """ 2761 cdef Py_ssize_t dx = len(x)-1, dy = len(y) - 1 2762 if dx < 0 or dy < 0: return [ ] 2763 cdef Py_ssize_t i, j, start, end = dx if dx < self._order else self._order-1 2764 cdef list twists = [ y ] 2765 for i from 0 <= i < end: 2766 twists.append([ self._parent.twist_map()(a) for a in twists[i] ]) 2767 cdef list res = [ ] 2768 for j from 0 <= j <= dx+dy: 2769 start = 0 if j <= dy else j-dy 2770 end = j if j <= dx else dx 2771 sum = x[start] * twists[start % self._order][j-start] 2772 for i from start < i <= end: 2773 sum += x[i] * twists[i % self._order][j-i] 2774 res.append(sum) 2775 return res 2776 2777 2778 cdef list mul_step_matrix(self, list x, list y): 2779 r""" 2780 Multiplication step in Karatsuba-matrix algorithm. It 2781 is based on the ring isomorphism: 2782 2783 .. MATH:: 2784 2785 S / NS \simeq M_r(k^\sigma) 2786 2787 with the standard notations:: 2788 2789 - `S` is the underlying skew polynomial ring 2790 2791 - `k` is the base ring of `S` (it's a finite field) 2792 2793 - `\sigma` is the twisting automorphism on `k` 2794 2795 - `r` is the order of `\sigma` 2796 2797 - `N` is a polynomial of definition of the extension 2798 `k / k^\sigma` 2799 2800 .. WARNING:: 2801 2802 The polynomials `x` and `y` must have degree 2803 at most `r^2/2` 2804 """ 2805 cdef Py_ssize_t dx = len(x)-1, dy = len(y) - 1 2806 cdef Py_ssize_t i, j 2807 cdef Py_ssize_t r = self._order 2808 cdef Mx, My, M 2809 cdef list row 2810 cdef RingElement c 2811 k = self._parent.base_ring() 2812 Mx = self._T * matrix(k, r, r, x + [self._zero] * (r*r-len(x))) 2813 My = self._T * matrix(k, r, r, y + [self._zero] * (r*r-len(y))) 2814 for i from 0 <= i < r: 2815 frob = self._parent.twist_map(i) 2816 row = Mx.row(i).list() 2817 row = map(frob,row) 2818 row = [ self._t*c for c in row[r-i:] ] + row[:r-i] 2819 Mx.set_row(i,row) 2820 row = My.row(i).list() 2821 row = map(frob,row) 2822 row = [ self._t*c for c in row[r-i:] ] + row[:r-i] 2823 My.set_row(i,row) 2824 M = Mx * My 2825 for i from 0 <= i < r: 2826 frob = self._parent.twist_map(-i) 2827 row = M.row(i).list() 2828 row = row[i:] + [ c/self._t for c in row[:i] ] 2829 row = map(frob,row) 2830 M.set_row(i,row) 2831 M = self._Tinv * M 2832 return M.list()[:dx+dy+1] 2833 2834 2835 cdef list mul_iter(self, list x, list y, char flag): # Assume dx >= dy 2836 """ 2837 Karatsuba's recursive iteration 2838 2839 .. WARNING:: 2840 2841 This function assumes that len(x) >= len(y). 2842 """ 2843 cdef Py_ssize_t i, j, k 2844 cdef Py_ssize_t dx = len(x)-1, dy = len(y)-1 2845 if self._algo_matrix: 2846 if dx < self._cutoff: 2847 if flag: 2848 return self.mul_step_matrix(y,x) 2849 else: 2850 return self.mul_step_matrix(x,y) 2851 else: 2852 if dy < self._cutoff: 2853 if flag: 2854 return self.mul_step(y,x) 2855 else: 2856 return self.mul_step(x,y) 2857 cdef Py_ssize_t dp = self._order * (1 + ceil(dx/self._order/2)) 2858 cdef list x1 = x[:dp], x2 = x[dp:] 2859 cdef list y1 = y[:dp], y2 = y[dp:] 2860 cdef list p1, p2, res 2861 res = self.mul_iter(x1,y1,flag) 2862 if dy >= dp: 2863 res.append(self._zero) 2864 p1 = self.mul_iter(x2,y2,flag) 2865 res.extend(p1) 2866 for i from 0 <= i <= dx-dp: x1[i] += x2[i] 2867 for i from 0 <= i <= dy-dp: y1[i] += y2[i] 2868 p2 = self.mul_iter(x1,y1,flag) 2869 j = dx + dy - 2*dp 2870 k = min(2*dp-2, len(res)-dp-1) 2871 for i from k >= i > j: 2872 res[i+dp] += p2[i] - res[i] 2873 for i from j >= i >= 0: 2874 res[i+dp] += p2[i] - p1[i] - res[i] 2875 else: 2876 if dx-dp < dy: 2877 p2 = self.mul_iter(y1,x2,not flag) 2878 else: 2879 p2 = self.mul_iter(x2,y1,flag) 2880 for i from 0 <= i < dy: 2881 res[i+dp] += p2[i] 2882 res.extend(p2[dy:]) 2883 return res 2884 2885 2886 def div(self,left,right): 2887 """ 2888 INPUT: 2889 2890 - ``left`` -- a skew polynomial in the skew polynomial 2891 right attached to the class 2892 2893 - ``right`` -- a skew polynomial in the skew polynomial 2894 right attached to the class 2895 2896 OUTPUT: 2897 2898 The quotient and the remainder of the right euclidien division 2899 of left by right (computed by Karatsuba's algorithm). 2900 2901 .. WARNING:: 2902 2903 This algorithm is only efficient for very very large 2904 degrees. Do not use it! 2905 2906 .. TODO:: 2907 2908 Try to understand why... 2909 2910 EXAMPLES:: 2911 2912 sage: k.<t> = GF(5^3) 2913 sage: Frob = k.frobenius_endomorphism() 2914 sage: S.<x> = k['x',Frob] 2915 sage: KarClass = S._karatsuba_class 2916 2917 sage: a = S.random_element(degree=1000) 2918 sage: b = S.random_element(degree=500) 2919 sage: q,r = KarClass.div(a,b) 2920 sage: q2,r2 = a.quo_rem(b) 2921 sage: q == q2 2922 True 2923 sage: r == r2 2924 True 2925 """ 2926 cdef list a = left.list() 2927 cdef list b = right.list() 2928 cdef Py_ssize_t da = len(a)-1, db = len(b)-1, i 2929 cdef list q 2930 self._twinv = [ ~b[db] ] 2931 for i from 0 <= i < min(da-db,self._order-1): 2932 self._twinv.append(self._parent.twist_map()(self._twinv[i])) 2933 if self._cutoff < 0: 2934 q = self.div_step(a,0,da,b,0,db) 2935 else: 2936 q = self.div_iter(a,0,da,b,0,db) 2937 return self._parent(q), self._parent(a[:db]) 2938 2939 2940 cdef list div_step(self, list a, Py_ssize_t ia, Py_ssize_t da, list b, Py_ssize_t ib, Py_ssize_t db): 2941 """ 2942 Division step in Karatsuba's algorithm 2943 """ 2944 cdef Py_ssize_t i, j 2945 cdef list q = [ ] 2946 cdef list twb = [ b[ib:ib+db] ] 2947 for i from 0 <= i < min(da-db,self._order-1): 2948 twb.append([ self._parent.twist_map()(x) for x in twb[i] ]) 2949 for i from da-db >= i >= 0: 2950 c = self._twinv[i%self._order] * a[ia+i+db] 2951 for j from 0 <= j < db: 2952 a[ia+i+j] -= c * twb[i%self._order][j] 2953 q.append(c) 2954 q.reverse() 2955 return q 2956 2957 2958 cdef list div_iter(self, list a, Py_ssize_t ia, Py_ssize_t da, list b, Py_ssize_t ib, Py_ssize_t db): 2959 """ 2960 Karatsuba recursive iteration 2961 """ 2962 cdef Py_ssize_t delta = da - db 2963 if delta < self._cutoff: 2964 return self.div_step(a,ia,da,b,ib,db) 2965 cdef Py_ssize_t i 2966 cdef Py_ssize_t dp = self._order * (1 + int(delta/self._order/2)) 2967 cdef list q = self.div_iter(a,ia+db,delta,b,ib+db-dp,dp), pr 2968 if db < delta: 2969 pr = self.mul_iter(q,b[ib:ib+db-dp],0) 2970 else: 2971 pr = self.mul_iter(b[ib:ib+db-dp],q,1) 2972 for i from 0 <= i < len(pr): 2973 a[i+ia+dp] -= pr[i] 2974 cdef list qq = self.div_iter(a,ia,db+dp-1,b,ib,db) 2975 qq.extend(q) 2976 return qq -
new file sage/rings/polynomial/skew_polynomial_ring.py
diff --git a/sage/rings/polynomial/skew_polynomial_ring.py b/sage/rings/polynomial/skew_polynomial_ring.py new file mode 100644
- + 1 r""" 2 Skew Univariate Polynomial Rings 3 4 Sage implements dense skew univariate polynomials over commutative rings. 5 6 DEFINITION: 7 8 Given a ring `R` and a ring endomorphism `\sigma` of `R`, the ring of 9 skew polynomials `R[x,\sigma]` is the usual abelian group polynomial 10 `R[x]` equipped with the modification multiplication deduced from the 11 rule `X a = \sigma(a) X`. 12 13 .. TODO:: 14 15 Add derivations. 16 17 EXAMPLES:: 18 19 sage: R.<t> = ZZ[] 20 sage: sigma = R.hom([t+1]) 21 sage: S.<x> = SkewPolynomialRing(R,sigma); S 22 Skew Polynomial Ring in x over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 23 24 One can also use a shorter syntax:: 25 26 sage: S.<x> = R['x',sigma]; S 27 Skew Polynomial Ring in x over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 28 29 Be careful, with the latter syntax, one cannot omit the name of the 30 variable neither in LHS nor in RHS. If we omit it in LHS, the variable 31 is not created:: 32 33 sage: Sy = R['y',sigma]; Sy 34 Skew Polynomial Ring in y over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 35 sage: y.parent() 36 Traceback (most recent call last): 37 ... 38 NameError: name 'y' is not defined 39 40 If we omit it in RHS, sage tries to create a polynomial ring and fails:: 41 42 sage: Sz.<z> = R[sigma] 43 Traceback (most recent call last): 44 ... 45 ValueError: variable names must be alphanumeric, but one is 'Ring endomorphism of Univariate Polynomial Ring in t over Integer Ring 46 Defn: t |--> t + 1' which is not. 47 48 As for polynomials, skew polynomial rings with different variable names 49 are not equal:: 50 51 sage: R['x',sigma] == R['y',sigma] 52 False 53 54 Of course, skew polynomial rings with different twist maps are not 55 equal as well 56 57 sage: R['x',sigma] == R['x',sigma^2] 58 False 59 60 Saving and loading of polynomial rings works:: 61 62 sage: loads(dumps(R['x',sigma])) == R['x',sigma] 63 True 64 65 There is a coercion map from the base ring of the skew polynomial rings:: 66 67 sage: S.has_coerce_map_from(R) 68 True 69 sage: x.parent() 70 Skew Polynomial Ring in x over Univariate Polynomial Ring in t over Integer Ring twisted by t |--> t + 1 71 sage: t.parent() 72 Univariate Polynomial Ring in t over Integer Ring 73 sage: y = x+t; y 74 x + t 75 sage: y.parent() is S 76 True 77 78 AUTHOR: 79 80 - Xavier Caruso (2012-06-29) 81 """ 82 83 ############################################################################# 84 # Copyright (C) 2012 Xavier Caruso <xavier.caruso@normalesup.org> 85 # 86 # Distributed under the terms of the GNU General Public License (GPL) 87 # 88 # http://www.gnu.org/licenses/ 89 #**************************************************************************** 90 91 from sage.structure.unique_representation import UniqueRepresentation 92 93 from sage.rings.infinity import Infinity 94 95 from sage.structure.element import Element 96 import sage.algebras.algebra 97 import sage.categories.basic as categories 98 import sage.rings.ring as ring 99 import sage.rings.ring_element as ring_element 100 import sage.rings.integral_domain as integral_domain 101 import sage.rings.principal_ideal_domain as principal_ideal_domain 102 import sage.rings.polynomial.polynomial_element_generic as polynomial_element_generic 103 import sage.rings.rational_field as rational_field 104 from sage.rings.integer_ring import is_IntegerRing, IntegerRing 105 from sage.rings.integer import Integer 106 from sage.libs.pari.all import pari_gen 107 108 from sage.structure.parent_gens import normalize_names 109 from sage.misc.prandom import randint 110 111 from sage.categories.morphism import Morphism 112 from sage.categories.morphism import IdentityMorphism 113 from sage.categories.homset import Hom 114 from sage.categories.map import Section 115 from sage.rings.morphism import RingHomomorphism 116 117 from sage.matrix.matrix_space import MatrixSpace 118 119 import sage.misc.latex as latex 120 121 from sage.rings.polynomial.polynomial_ring import PolynomialRing_general 122 from sage.rings.polynomial.polynomial_element import Polynomial_generic_dense 123 from sage.rings.polynomial.polynomial_element import PolynomialBaseringInjection 124 import copy, re 125 126 127 def is_SkewPolynomialRing(S): 128 """ 129 Return True if S is a skew polynomial ring. 130 131 EXAMPLES:: 132 133 sage: from sage.rings.polynomial.skew_polynomial_ring import is_SkewPolynomialRing 134 sage: is_SkewPolynomialRing(QQ['x']) 135 False 136 137 sage: k.<t> = GF(5^3) 138 sage: Frob = k.frobenius_endomorphism() 139 sage: is_SkewPolynomialRing(k['x',Frob]) 140 True 141 """ 142 return isinstance(S, SkewPolynomialRing_general) 143 144 145 ######################################################################################### 146 147 class SectionSkewPolynomialCenterInjection(Section): 148 def _call_ (self, x): 149 order = self.inverse()._order 150 section = self.inverse()._embed.section() 151 lx = x.list() 152 l = [ ] 153 mod = 0 154 for c in lx: 155 if mod == 0: 156 l.append(section(c)) 157 else: 158 if not c.is_zero(): 159 raise ValueError("%s is not in the center" % x) 160 mod += 1 161 if mod == order: 162 mod = 0 163 return self.codomain()(l) 164 165 166 class SkewPolynomialCenterInjection(RingHomomorphism): 167 def __init__(self,domain,codomain,embed,order): 168 RingHomomorphism.__init__(self,Hom(domain,codomain)) 169 self._embed = embed 170 self._order = order 171 self._codomain = codomain 172 self._section = SectionSkewPolynomialCenterInjection(self) 173 174 def _repr_(self): 175 return "Embedding of the center of %s into this ring" % self._codomain 176 177 def _call_(self,x): 178 k = self._codomain.base_ring () 179 l = [ ] 180 lz = [ k(0) ] * (self._order-1) 181 for c in x.list(): 182 l += [ self._embed(c) ] + lz 183 return self._codomain (l) 184 185 def section(self): 186 return self._section 187 188 189 class CenterSkewPolynomialRing(PolynomialRing_general): 190 """ 191 A specific class for the center of a skew polynomial ring. 192 """ 193 194 def __init__ (self, skew_ring, names=None, sparse=False, element_class=None): 195 if not isinstance (skew_ring, SkewPolynomialRing_general): 196 raise TypeError("%s is not a Skew Polynomial Ring" % skew_ring) 197 self._skew_ring = skew_ring 198 base_ring = skew_ring.base_ring() 199 kfixed, embed = skew_ring._map.fixed_points() 200 self._embed_basering = embed 201 order = skew_ring._map.order() 202 if order == Infinity: 203 raise NotImplementedError 204 self.__is_sparse = sparse 205 self._PolynomialRing_general__is_sparse = sparse 206 if element_class: 207 self._polynomial_class = element_class 208 else: 209 if sparse: 210 raise NotImplementedError("sparse skew polynomials are not implemented") 211 else: 212 self._polynomial_class = sage.rings.polynomial.skew_polynomial_element.CenterSkewPolynomial_generic_dense 213 214 # Algebra.__init__ also calls __init_extra__ of Algebras(...).parent_class, which 215 # tries to provide a conversion from the base ring, if it does not exist. 216 # This is for algebras that only do the generic stuff in their initialisation. 217 # But here, we want to use PolynomialBaseringInjection. Hence, we need to 218 # wipe the memory and construct the conversion from scratch. 219 sage.algebras.algebra.Algebra.__init__(self, kfixed, names=names, normalize=True, category=None) 220 221 if names is None: 222 if order == 1: 223 self._variable_name = skew_ring.variable_name() 224 self._latex_variable_name = skew_ring.latex_variable_names()[0] 225 self._parenthesis = False 226 else: 227 self._variable_name = skew_ring.variable_name () + "^" + str(order) 228 self._latex_variable_name = skew_ring.latex_variable_names()[0] + "^{" + str (order) + "}" 229 self._parenthesis = True 230 else: 231 self._variable_name = sage.algebras.algebra.Algebra.variable_name(self) 232 self._latex_variable_name = sage.algebras.algebra.Algebra.latex_variable_names(self)[0] 233 self._parenthesis = False 234 self._names = [ self._variable_name ] 235 self.__generator = self._polynomial_class (self, [0,1], is_gen=True) 236 base_inject = PolynomialBaseringInjection(kfixed,self) 237 center_inject = SkewPolynomialCenterInjection (self, skew_ring, embed, order) 238 self._unset_coercions_used() 239 self._populate_coercion_lists_( 240 coerce_list = [base_inject], 241 convert_list = [list, base_inject], 242 embedding = center_inject) 243 244 def _repr_ (self): 245 """ 246 Return a string representation of this ring. 247 """ 248 s = "Center of %s:\n" % self._skew_ring 249 s += PolynomialRing_general._repr_(self) 250 return s 251 252 def _latex_ (self): 253 """ 254 Return a latex representation of this ring. 255 """ 256 return "%s[%s]"%(latex.latex(self.base_ring()), self._latex_variable_name) 257 258 def gen (self,n=0): 259 """ 260 Return the generator of this ring. 261 262 EXAMPLES:: 263 264 sage: k.<t> = GF(2^10) 265 &n