Revision 3321:ad80c5ad0ca0, 65.3 KB checked in by David Roe <roed@…>, 6 years ago (diff)

Working on extended rationals.

Line
1r"""
2Elements
3
4AUTHORS:
5   -- David Harvey (2006-10-16): changed CommutativeAlgebraElement to derive
6   from CommutativeRingElement instead of AlgebraElement
7   -- David Harvey (2006-10-29): implementation and documentation of new
8   arithmetic architecture
9   -- William Stein (2006-11): arithmetic architecture -- pushing it through to completion.
10
11
12\subsection{The Abstract Element Class Heierarchy}
13This is the abstract class heierchary, i.e., these are all
14abstract base classes.
15\begin{verbatim}
16SageObject
17    Element
18        ModuleElement
20
21        MonoidElement
22            MultiplicativeGroupElement
23
24        RingElement
25            CommutativeRingElement
26                IntegralDomainElement
27                    DedekindDomainElement
28                        PrincipalIdealDomainElement
29                            EuclideanDomainElement
30            FieldElement
31                FiniteFieldElement
32            AlgebraElement   (note -- can't derive from module, since no multiple inheritence)
33                CommutativeAlgebraElement
34            InfinityElement
35                PlusInfinityElement
36                MinusInfinityElement
37\end{verbatim}
38
39\subsection{How to Define a New Element Class}
40Elements typically define a method \code{_new_c}, e.g.,
41\begin{verbatim}
42    cdef _new_c(self, defining data):
43        cdef FreeModuleElement_generic_dense x
44        x = PY_NEW(FreeModuleElement_generic_dense)
45        x._parent = self._parent
46        x._entries = v
47\end{verbatim}
48that creates a new sibling very quickly from defining data
49with assumed properties.
50
51SAGE has a special system in place for handling arithmetic operations
52for all Element subclasses. There are various rules that must be
53followed by both arithmetic implementors and callers.
54
55A quick summary for the impatient:
56\begin{itemize}
57 \item DO NOT OVERRIDE _add_c. EVER. THANKS.
58 \item DO NOT CALL _add_c_impl DIRECTLY.
59 \item To implement addition for a python class, override def _add_().
60 \item To implement addition for a pyrex class, override cdef _add_c_impl().
61 \item If you want to add x and y, whose parents you know are IDENTICAL,
62   you may call _add_(x, y) (from python or pyrex) or _add_c(x, y) (from
63   pyrex -- this will be faster). This will be the fastest way to guarantee
64   that the correct implementation gets called. Of course you can still
65   always use "x + y".
66\end{itemize}
67
68Now in more detail. The aims of this system are to provide (1) an efficient
69calling protocol from both python and pyrex, (2) uniform coercion semantics
70across SAGE, (3) ease of use, (4) readability of code.
71
72We will take addition of RingElements as an example; all other operators
73and classes are similar. There are four relevant functions.
74
76
77   This function is called by python or pyrex when the binary "+" operator
78   is encountered. It ASSUMES that at least one of its arguments is a
79   RingElement; only a really twisted programmer would violate this
80   condition. It has a fast pathway to deal with the most common case
81   where the arguments have the same parent. Otherwise, it uses the coercion
82   module to work out how to make them have the same parent. After any
83   necessary coercions have been performed, it calls _add_c to dispatch to
84   the correct underlying addition implementation.
85
86   Note that although this function is declared as def, it doesn't have the
87   usual overheads associated with python functions (either for the caller
88   or for __add__ itself). This is because python has optimised calling
89   protocols for such special functions.
90
92
93   DO ***NOT*** OVERRIDE THIS FUNCTION.
94
95   The two arguments to this function MUST have the SAME PARENT.
96   Its return value MUST have the SAME PARENT as its arguments.
97
98   If you want to add two objects from pyrex, and you know that their
99   parents are the same object, you are encouraged to call this function
100   directly, instead of using "x + y".
101
103   It takes necessary steps to decide whether a pyrex implementation of
105   The code is optimised in favour of reaching _add_c_impl as soon as
106   possible.
107
109
110   This is the function you should override to implement addition in a
111   python subclass of RingElement.
112
113   WARNING: if you override this in a *pyrex* class, it won't get called.
115   keep this in mind whenever you move a class down from python to pyrex.
116
117   The two arguments to this function are guaranteed to have the
118   SAME PARENT. Its return value MUST have the SAME PARENT as its
119   arguments.
120
121   If you want to add two objects from python, and you know that their
122   parents are the same object, you are encouraged to call this function
123   directly, instead of using "x + y".
124
125   The default implementation of this function is to call _add_c_impl,
126   so if no-one has defined a python implementation, the correct pyrex
127   implementation will get called.
128
130
131   This is the function you should override to implement addition in a
132   pyrex subclass of RingElement.
133
134   The two arguments to this function are guaranteed to have the
135   SAME PARENT. Its return value MUST have the SAME PARENT as its
136   arguments.
137
138   DO ***NOT*** CALL THIS FUNCTION DIRECTLY.
139
140   (Exception: you know EXACTLY what you are doing, and you know exactly
141   which implementation you are trying to call; i.e. you're not trying to
142   write generic code. In particular, if you call this directly, your code
143   will not work correctly if you run it on a python class derived from
144   a pyrex class where someone has redefined _add_ in python.)
145
146   The default implementation of this function is to raise a
147   NotImplementedError, which will happen if no-one has supplied
149
150"""
151
152
153##################################################################
154# Generic element, so all this functionality must be defined
155# by any element.  Derived class must call __init__
156##################################################################
157
158include "../ext/cdefs.pxi"
159include "../ext/stdsage.pxi"
160include "../ext/python.pxi"
161
162import operator
163
164from sage.structure.parent      cimport Parent
165from sage.structure.parent_base cimport ParentWithBase
166
167# This classes uses element.pxd.  To add data members, you
168# must change that file.
169
170def make_element(_class, _dict, parent):
171    """
172    Used for unpickling Element objects (and subclasses).
173
174    This should work for any Python class deriving from Element, as long
175    as it doesn't implement some screwy __new__() method.
176
178    """
179    new_object = _class.__new__(_class)
180    new_object._set_parent(parent)
181    new_object.__dict__ = _dict
182    return new_object
183
184
185
186def is_Element(x):
187    """
188    Return True if x is of type Element.
189
190    EXAMPLES:
191        sage: is_Element(2/3)
192        True
193        sage: is_Element(QQ^3)
194        False
195    """
196    return IS_INSTANCE(x, Element)
197
198cdef class Element(sage_object.SageObject):
199    """
200    Generic element of a structure. All other types of elements
201    (RingElement, ModuleElement, etc) derive from this type.
202
203    Subtypes must either call __init__() to set _parent, or may
204    set _parent themselves if that would be more efficient.
205    """
206
207    def __init__(self, parent):
208        r"""
209        INPUT:
210            parent -- a SageObject
211        """
212        #if parent is None:
213        #    raise RuntimeError, "bug -- can't set parent to None"
214        self._parent = parent
215
216    def _set_parent(self, parent):
217        r"""
218        INPUT:
219            parent -- a SageObject
220        """
221        self._parent = parent
222
223    cdef _set_parent_c(self, ParentWithBase parent):
224        self._parent = parent
225
226    def _repr_(self):
227        return "Generic element of a structure"
228
229    def __reduce__(self):
230        return (make_element, (self.__class__, self.__dict__, self._parent))
231
232    def __hash__(self):
233        return hash(str(self))
234
235    def _im_gens_(self, codomain, im_gens):
236        """
237        Return the image of self in codomain under the map that sends
238        the images of the generators of the parent of self to the
239        tuple of elements of im_gens.
240        """
241        raise NotImplementedError
242
243    cdef base_extend_c(self, ParentWithBase R):
244        if HAS_DICTIONARY(self):
245            return self.base_extend(R)
246        else:
247            return self.base_extend_c_impl(R)
248
249    cdef base_extend_c_impl(self, ParentWithBase R):
250        cdef ParentWithBase V
251        V = self._parent.base_extend(R)
252        return (<Parent>V)._coerce_c(self)
253
254    def base_extend(self, R):
255        return self.base_extend_c_impl(R)
256
257    def base_ring(self):
258        """
259        Returns the base ring of this element's parent (if that makes sense).
260        """
261        return self._parent.base_ring()
262
263    def category(self):
264        from sage.categories.category import Elements
265        return Elements(self._parent)
266
267    def parent(self, x=None):
268        """
269        Returns parent of this element; or, if the optional argument x is
270        supplied, the result of coercing x into the parent of this element.
271        """
272        if x is None:
273            return self._parent
274        else:
275            return self._parent(x)
276
277    def __xor__(self, right):
278        raise RuntimeError, "Use ** for exponentiation, not '^', which means xor\n"+\
279              "in Python, and has the wrong precedence."
280
281    def _coeff_repr(self, no_space=True):
282        if self._is_atomic():
283            s = str(self)
284        else:
285            s = "(%s)"%self
286        if no_space:
287            return s.replace(' ','')
288        return s
289
290    def _latex_coeff_repr(self):
291        try:
292            s = self._latex_()
293        except AttributeError:
294            s = str(self)
295        if self._is_atomic():
296            return s
297        else:
298            return "\\left(%s\\right)"%s
299
300    def _is_atomic(self):
301        """
302        Return True if and only if parenthesis are not required when
303        *printing* out any of $x - s$, $x + s$, $x^s$ and $x/s$.
304
305        EXAMPLES:
306            sage: n = 5; n._is_atomic()
307            True
308            sage: n = x+1; n._is_atomic()
309            False
310        """
311        if self._parent.is_atomic_repr():
312            return True
313        s = str(self)
314        return PyBool_FromLong(s.find("+") == -1 and s.find("-") == -1 and s.find(" ") == -1)
315
316    def is_zero(self):
317        return PyBool_FromLong(self == self._parent(0))
318
319    def _richcmp_(left, right, op):
320        return left._richcmp(right, op)
321
322    cdef _richcmp(left, right, int op):
323        """
324        Compare left and right.
325        """
326        cdef int r
327        if not have_same_parent(left, right):
328            try:
329                _left, _right = canonical_coercion_c(left, right)
330                r = cmp(_left, _right)
331            except TypeError:
332                r = cmp(type(left), type(right))
333                if r == 0:
334                    r = -1
335        else:
336            if HAS_DICTIONARY(left):   # fast check
337                r = left.__cmp__(right)
338            else:
339                r = left._cmp_c_impl(right)
340
341        return left._rich_to_bool(op, r)
342
343    cdef _rich_to_bool(self, int op, int r):
344        if op == 0:  #<
345            return PyBool_FromLong(r  < 0)
346        elif op == 2: #==
347            return PyBool_FromLong(r == 0)
348        elif op == 4: #>
349            return PyBool_FromLong(r  > 0)
350        elif op == 1: #<=
351            return PyBool_FromLong(r <= 0)
352        elif op == 3: #!=
353            return PyBool_FromLong(r != 0)
354        elif op == 5: #>=
355            return PyBool_FromLong(r >= 0)
356
357    ####################################################################
358    # For a derived Pyrex class, you **must** put the following in
359    # your subclasses, in order for it to take advantage of the
360    # above generic comparison code.  You must also define
361    # _cmp_c_impl.
362    # This is simply how Python works.
363    #
364    # For a *Python* class just define __cmp__ as always.
365    # But note that when this get called you can assume that
366    # both inputs have identical parents.
367    ####################################################################
368    def __richcmp__(left, right, int op):
369        return (<Element>left)._richcmp(right, op)
370
371    cdef int _cmp_c_impl(left, Element right) except -2:
372        ### For derived SAGEX code, you *MUST* ALSO COPY the __richcmp__ above
373        ### into your class!!!  For Python code just use __cmp__.
374        raise NotImplementedError, "BUG: sort algorithm for elements of '%s' not implemented"%right.parent()
375
376    def __cmp__(left, right):
377        return left._cmp_c_impl(right)   # default
378
379def is_ModuleElement(x):
380    """
381    Return True if x is of type ModuleElement.
382
383    This is even faster than using isinstance inline.
384
385    EXAMPLES:
386        sage: is_ModuleElement(2/3)
387        True
388        sage: is_ModuleElement((QQ^3).0)
389        True
390        sage: is_ModuleElement('a')
391        False
392    """
393    return IS_INSTANCE(x, ModuleElement)
394
395
396cdef class ModuleElement(Element):
397    """
398    Generic element of a module.
399    """
400    ##################################################
401    def is_zero(self):
402        return PyBool_FromLong(self == self._parent(0))
403
404    ##################################################
406    ##################################################
408        """
409        Top-level addition operator for ModuleElements.
410
411        See extensive documentation at the top of element.pyx.
412        """
413
414        # Try fast pathway if they are both ModuleElements and the parents
415        # match.
416
417        # (We know at least one of the arguments is a ModuleElement. So if
418        # their types are *equal* (fast to check) then they are both
419        # ModuleElements. Otherwise use the slower test via PY_TYPE_CHECK.)
420        if have_same_parent(left, right):
423
424    cdef ModuleElement _add_c(left, ModuleElement right):
425        """
427
428        DO NOT OVERRIDE THIS FUNCTION.
429
430        See extensive documentation at the top of element.pyx.
431        """
432        if HAS_DICTIONARY(left):   # fast check
433            # TODO: this bit will be unnecessarily slow if someone derives
434            # from the pyrex class *without* overriding _add_, since then
435            # we'll be making an unnecessary python call to _add_, which will
436            # end up in _add_c_impl anyway. There must be a simple way to
437            # distinguish this situation. It's complicated because someone
438            # can even override it at the instance level (without overriding
439            # it in the class.)
441        else:
442            # Must be a pure Pyrex class.
444
445
446    cdef ModuleElement _add_c_impl(left, ModuleElement right):
447        """
448        Pyrex classes should override this function to implement addition.
449        DO NOT CALL THIS FUNCTION DIRECTLY.
450        See extensive documentation at the top of element.pyx.
451        """
452        raise TypeError, arith_error_message(left, right, operator.add)
453
454
455    def _add_(ModuleElement left, ModuleElement right):
456        """
457        Python classes should override this function to implement addition.
458
459        See extensive documentation at the top of element.pyx.
460        """
462
463    ##################################################
464    # Subtraction
465    ##################################################
466
467    def __sub__(left, right):
468        """
469        Top-level subtraction operator for ModuleElements.
470        See extensive documentation at the top of element.pyx.
471        """
472        if have_same_parent(left, right):
473            return (<ModuleElement>left)._sub_c(<ModuleElement>right)
474        return bin_op_c(left, right, operator.sub)
475
476    cdef ModuleElement _sub_c(left, ModuleElement right):
477        """
478        Subtraction dispatcher for ModuleElements.
479
480        DO NOT OVERRIDE THIS FUNCTION.
481
482        See extensive documentation at the top of element.pyx.
483        """
484
485        if HAS_DICTIONARY(left):   # fast check
486            return left._sub_(right)
487        else:
488            # Must be a pure Pyrex class.
489            return left._sub_c_impl(right)
490
491
492    cdef ModuleElement _sub_c_impl(left, ModuleElement right):
493        """
494        Pyrex classes should override this function to implement subtraction.
495
496        DO NOT CALL THIS FUNCTION DIRECTLY.
497
498        See extensive documentation at the top of element.pyx.
499        """
500        # default implementation is to use the negation and addition
501        # dispatchers:
503
504
505    def _sub_(ModuleElement left, ModuleElement right):
506        """
507        Python classes should override this function to implement subtraction.
508
509        See extensive documentation at the top of element.pyx.
510        """
511        return left._sub_c_impl(right)
512
513    ##################################################
514    # Negation
515    ##################################################
516
517    def __neg__(self):
518        """
519        Top-level negation operator for ModuleElements.
520        See extensive documentation at the top of element.pyx.
521        """
522        # We ASSUME that self is a ModuleElement. No type checks.
523        return (<ModuleElement>self)._neg_c()
524
525
526    cdef ModuleElement _neg_c(self):
527        """
528        Negation dispatcher for ModuleElements.
529        DO NOT OVERRIDE THIS FUNCTION.
530        See extensive documentation at the top of element.pyx.
531        """
532
533        if HAS_DICTIONARY(self):   # fast check
534            return self._neg_()
535        else:
536            # Must be a pure Pyrex class.
537            return self._neg_c_impl()
538
539
540    cdef ModuleElement _neg_c_impl(self):
541        """
542        Pyrex classes should override this function to implement negation.
543        DO NOT CALL THIS FUNCTION DIRECTLY.
544        See extensive documentation at the top of element.pyx.
545        """
546        # default implementation is to try multiplying by -1.
547        return bin_op_c(self._parent._base(-1), self, operator.mul)
548
549
550    def _neg_(ModuleElement self):
551        """
552        Python classes should override this function to implement negation.
553
554        See extensive documentation at the top of element.pyx.
555        """
556        return self._neg_c_impl()
557
558    ##################################################
559    # Module element multiplication (scalars, etc.)
560    ##################################################
561    def __mul__(left, right):
562        return module_element_generic_multiply_c(left, right)
563
564    cdef ModuleElement _multiply_by_scalar(self, right):
565        # self * right,  where right need not be a ring element in the base ring
566        # This does type checking and canonical coercion then calls _lmul_c_impl.
567        if PY_TYPE_CHECK(right, Element):
568            if (<Element>right)._parent is self._parent._base:
569                # No coercion needed
570                return self._lmul_c(right)
571            else:
572                # Otherwise we do an explicit canonical coercion.
573                try:
574                    return self._lmul_c( self._parent._base._coerce_c(right) )
575                except TypeError:
576                    # that failed -- try to base extend right then do the multiply:
577                    self = self.base_extend((<RingElement>right)._parent)
578                    return (<ModuleElement>self)._lmul_c(right)
579        else:
580            # right is not an element at all
581            return (<ModuleElement>self)._lmul_c(self._parent._base._coerce_c(right))
582
583    cdef ModuleElement _rmultiply_by_scalar(self, left):
584        # left * self, where left need not be a ring element in the base ring
585        # This does type checking and canonical coercion then calls _rmul_c_impl.
586        #
587        # INPUT:
588        #    self -- a module element
589        #    left -- a scalar
590        # OUTPUT:
591        #    left * self
592        #
593        if PY_TYPE_CHECK(left, Element):
594            if (<Element>left)._parent is self._parent._base:
595                # No coercion needed
596                return self._rmul_c(left)
597            else:
598                # Otherwise we do an explicit canonical coercion.
599                try:
600                    return self._rmul_c(self._parent._base._coerce_c(left))
601                except TypeError:
602                    # that failed -- try to base extend self then do the multiply:
603                    self = self.base_extend((<RingElement>left)._parent)
604                    return (<ModuleElement>self)._rmul_c(left)
605        else:
606            # now left is not an element at all.
607            return (<ModuleElement>self)._rmul_c(self._parent._base._coerce_c(left))
608
609    cdef ModuleElement _lmul_nonscalar_c(left, right):
610        # Compute the product left * right, where right is assumed to be a nonscalar (so no coercion)
611        # This is a last resort.
612        if HAS_DICTIONARY(left):
613            return left._lmul_nonscalar(right)
614        else:
615            return left._lmul_nonscalar_c_impl(right)
616
617    cdef ModuleElement _lmul_nonscalar_c_impl(left, right):
618        raise TypeError
619
620    def _lmul_nonscalar(left, right):
621        return left._lmul_nonscalar_c_impl(right)
622
623    cdef ModuleElement _rmul_nonscalar_c(right, left):
624        if HAS_DICTIONARY(right):
625            return right._rmul_nonscalar(left)
626        else:
627            return right._rmul_nonscalar_c_impl(left)
628
629    cdef ModuleElement _rmul_nonscalar_c_impl(right, left):
630        raise TypeError
631
632    def _rmul_nonscalar(right, left):
633        return right._rmul_nonscalar_c_impl(left)
634
635
636    # rmul -- left * self
637    cdef ModuleElement _rmul_c(self, RingElement left):
638        """
639        DO NOT OVERRIDE THIS FUNCTION.  OK to call.
640        """
641        if HAS_DICTIONARY(self):
642            return self._rmul_(left)
643        else:
644            return self._rmul_c_impl(left)
645
646    cdef ModuleElement _rmul_c_impl(self, RingElement left):
647        """
648        Default module left scalar multiplication, which is to try to
649        canonically coerce the scalar to the integers and do that
650        multiplication, which is always defined.
651        """
652        n = int(left)
653        if n != left:
654            raise TypeError, "left (=%s) must be an integer."%left
655        a = self
656        if n < 0:
657            a = -a
658            n = -n
659        sum = self._parent(0)
660        asum = a
661        while True:
662            if n&1 > 0: sum = sum + asum
663            n = n >> 1
664            if n != 0:
665                asum = asum + asum
666            else:
667                break
668        return sum
669
670    def _rmul_(self, left):
671        return self._rmul_c_impl(left)
672
673
674    # lmul -- self * right
675
676    cdef ModuleElement _lmul_c(self, RingElement right):
677        """
678        DO NOT OVERRIDE THIS FUNCTION.
679        """
680        if HAS_DICTIONARY(self):
681            return self._lmul_(right)
682        else:
683            return self._lmul_c_impl(right)
684
685    cdef ModuleElement _lmul_c_impl(self, RingElement right):
686        """
687        Default module right scalar multiplication, which is to try to
688        canonically coerce the scalar to the integers and do that
689        multiplication, which is always defined.
690        """
691        return self._rmul_c(right)
692
693    def _lmul_(self, right):
694        return self._lmul_c_impl(right)
695
696    cdef RingElement coerce_to_base_ring(self, x):
697        if PY_TYPE_CHECK(x, Element) and (<Element>x)._parent is self._parent._base:
698            return x
699        try:
700            return self._parent._base._coerce_c(x)
701        except AttributeError:
702            return self._parent._base(x)
703
704    ##################################################
705    # Other properties
706    ##################################################
708        """
709        Return the additive order of self.
710        """
712
714        """
715        Return the additive order of self.
716        """
717        raise NotImplementedError
718
719def module_element_generic_multiply(left, right):
720    return module_element_generic_multiply_c(left, right)
721
722cdef module_element_generic_multiply_c(left, right):
723    cdef int is_element
724    if PY_TYPE_CHECK(right, ModuleElement) and not PY_TYPE_CHECK(right, RingElement):
725        # do left * (a module element right)
726        is_element = PY_TYPE_CHECK(left, Element)
727        if is_element  and  (<Element>left)._parent is (<ModuleElement>right)._parent._base:
728            # No coercion needed
729            return (<ModuleElement>right)._rmul_c(left)
730        else:
731            try:
732                return (<ModuleElement>right)._rmul_nonscalar_c(left)
733            except TypeError:
734                pass
735            # Otherwise we have to do an explicit canonical coercion.
736            try:
737                return (<ModuleElement>right)._rmul_c(
738                    (<Parent>(<ModuleElement>right)._parent._base)._coerce_c(left))
739            except TypeError:
740                if is_element:
741                    # that failed -- try to base extend right then do the multiply:
742                    right = right.base_extend((<RingElement>left)._parent)
743                    return (<ModuleElement>right)._rmul_c(left)
744    else:
745        # do (module element left) * right
746        # This is the symmetric case of above.
747        is_element = PY_TYPE_CHECK(right, Element)
748        if is_element  and  (<Element>right)._parent is (<ModuleElement>left)._parent._base:
749            # No coercion needed
750            return (<ModuleElement>left)._lmul_c(right)
751        else:
752            try:
753                return (<ModuleElement>left)._lmul_nonscalar_c(right)
754            except TypeError:
755                pass
756            # Otherwise we have to do an explicit canonical coercion.
757            try:
758                return (<ModuleElement>left)._lmul_c(
759                    (<Parent>(<ModuleElement>left)._parent._base)._coerce_c(right))
760            except TypeError:
761                if is_element:
762                    # that failed -- try to base extend right then do the multiply:
763                    left = left.base_extend((<RingElement>right)._parent)
764                    return (<ModuleElement>left)._rmul_c(right)
765    raise TypeError
766
767########################################################################
768# Monoid
769########################################################################
770
771def is_MonoidElement(x):
772    """
773    Return True if x is of type MonoidElement.
774    """
775    return IS_INSTANCE(x, MonoidElement)
776
777cdef class MonoidElement(Element):
778    """
779    Generic element of a monoid.
780    """
781
782    #############################################################
783    # Multiplication
784    #############################################################
785    def __mul__(left, right):
786        """
787        Top-level multiplication operator for monoid elements.
788        See extensive documentation at the top of element.pyx.
789        """
790        if have_same_parent(left, right):
791            return (<MonoidElement>left)._mul_c(<MonoidElement>right)
792        try:
793            return bin_op_c(left, right, operator.mul)
794        except TypeError, msg:
795            if isinstance(left, (int, long)) and left==1:
796                return right
797            elif isinstance(right, (int, long)) and right==1:
798                return left
799            raise TypeError, msg
800
801
802
803
804    cdef MonoidElement _mul_c(left, MonoidElement right):
805        """
806        Multiplication dispatcher for RingElements.
807        DO NOT OVERRIDE THIS FUNCTION.
808        See extensive documentation at the top of element.pyx.
809        """
810        if HAS_DICTIONARY(left):   # fast check
811            return left._mul_(right)
812        else:
813            return left._mul_c_impl(right)
814
815
816    cdef MonoidElement _mul_c_impl(left, MonoidElement right):
817        """
818        Pyrex classes should override this function to implement multiplication.
819        DO NOT CALL THIS FUNCTION DIRECTLY.
820        See extensive documentation at the top of element.pyx.
821        """
822        raise TypeError
823
824    def _mul_(left, right):
825        return left._mul_c_impl(right)
826
827    #############################################################
828    # Other generic functions that should be available to
829    # any monoid element.
830    #############################################################
831    def order(self):
832        """
833        Return the multiplicative order of self.
834        """
835        return self.multiplicative_order()
836
837    def multiplicative_order(self):
838        """
839        Return the multiplicative order of self.
840        """
841        raise NotImplementedError
842
843    def __pow__(self, nn, dummy):
844        """
845        Retern the (integral) power of self.
846        """
847        cdef int cn
848
849        n = int(nn)
850        if n != nn:
851            raise NotImplementedError, "non-integral exponents not supported"
852
853        a = self
854        if n < 0:
855            n = -n
856            a = ~self
857
858        if n < 4:
859            # These cases will probably be called often
860            # and don't benifit from the code below
861            cn = n
862            if cn == 0:
863                return (<Element>a)._parent(1)
864            elif cn == 1:
865                return a
866            elif cn == 2:
867                return a*a
868            elif cn == 3:
869                return a*a*a
870
871        # One multiplication can be saved by starting with
872        # the smallest power needed rather than with 1
873        apow = a
874        while n&1 == 0:
875            apow = apow*apow
876            n = n >> 1
877        power = apow
878        n = n >> 1
879
880        while n != 0:
881            apow = apow*apow
882            if n&1 != 0: power = power*apow
883            n = n >> 1
884
885        return power
886
887
889    """
890    Return True if x is of type AdditiveGroupElement.
891    """
893
895    """
896    Generic element of an additive group.
897    """
898    def order(self):
899        """
900        Return additive order of element
901        """
903
904    def __invert__(self):
905        raise NotImplementedError, "multiplicative inverse not defined for additive group elements"
906
907    cdef ModuleElement _rmul_c_impl(self, RingElement left):
908        return self._lmul_c_impl(left)
909
910    cdef ModuleElement _lmul_c_impl(self, RingElement right):
911        cdef int m
912        try:
913            m = int(right)  # a little worrisome.
914            if m<0:
915                return (-self)*(-m)
916            if m==1:
917                return self
918            P = self.scheme()(0)
919            if m==0:
920                return P
921            power = P
922            i = 0
923            apow2 = self
924            while ((m>>i) > 0):
925                if((m>>i) & 1):
926                    power = power + apow2
927                apow2 = apow2 + apow2
928                i = i + 1
929            return power
930        except OverflowError:
931            m0 = int(right)
932            if m0<0:
933                return (-self)*(-m0)
934            if m0==1:
935                return self
936            P = self.scheme()(0)
937            if m0==0:
938                return P
939            power = P
940            i = 0
941            apow2 = self
942            while ((m0>>i) > 0):
943                if((m0>>i) & 1):
944                    power = power + apow2
945                apow2 = apow2 + apow2
946                i = i + 1
947            return power
948
949def is_MultiplicativeGroupElement(x):
950    """
951    Return True if x is of type MultiplicativeGroupElement.
952    """
953    return IS_INSTANCE(x, MultiplicativeGroupElement)
954
955cdef class MultiplicativeGroupElement(MonoidElement):
956    """
957    Generic element of a multiplicative group.
958    """
959    def order(self):
960        """
961        Return the multiplicative order of self.
962        """
963        return self.multiplicative_order()
964
966        raise ArithmeticError, "addition not defined in a multiplicative group"
967
968    def __div__(left, right):
969        if have_same_parent(left, right):
970            return left._div_(right)
971        return bin_op_c(left, right, operator.div)
972
973    cdef MultiplicativeGroupElement _div_c(self, MultiplicativeGroupElement right):
974        """
975        Multiplication dispatcher for MultiplicativeGroupElements.
976        DO NOT OVERRIDE THIS FUNCTION.
977        See extensive documentation at the top of element.pyx.
978        """
979        if HAS_DICTIONARY(self):   # fast check
980            return self._div_(right)
981        else:
982            return self._div_c_impl(right)
983
984    cdef MultiplicativeGroupElement _div_c_impl(self, MultiplicativeGroupElement right):
985        """
986        Pyrex classes should override this function to implement division.
987        DO NOT CALL THIS FUNCTION DIRECTLY.
988        See extensive documentation at the top of element.pyx.
989        """
990        return self._parent.fraction_field()(self, right)
991
992    def _div_(MultiplicativeGroupElement self, MultiplicativeGroupElement right):
993        """
994        Python classes should override this function to implement division.
995        """
996        return self._div_c_impl(right)
997
998
999    def __invert__(self):
1000        return 1/self
1001
1002
1003def is_RingElement(x):
1004    """
1005    Return True if x is of type RingElement.
1006    """
1007    return IS_INSTANCE(x, RingElement)
1008
1009cdef class RingElement(ModuleElement):
1010    ##################################################
1011    def is_zero(self):
1012        return PyBool_FromLong(self == self.parent()(0))
1013
1014    def is_one(self):
1015        return PyBool_FromLong(self == self.parent()(1))
1016
1017    ##################################
1018    # Multiplication
1019    ##################################
1020
1021    # The default behavior for scalars is just to coerce into the parent ring.
1022    cdef ModuleElement _lmul_c_impl(self, RingElement right):
1023        return self._mul_c(<RingElement>(self._parent(right)))
1024
1025    cdef ModuleElement _rmul_c_impl(self, RingElement left):
1026        return (<RingElement>(self._parent)(left))._mul_c(self)
1027
1028    def __mul__(self, right):
1029        """
1030        Top-level multiplication operator for ring elements.
1031        See extensive documentation at the top of element.pyx.
1032        """
1033        # Try fast pathway if they are both RingElements and the parents match.
1034        # (We know at least one of the arguments is a RingElement. So if their
1035        # types are *equal* (fast to check) then they are both RingElements.
1036        # Otherwise use the slower test via PY_TYPE_CHECK.)
1037        if have_same_parent(self, right):
1038            return (<RingElement>self)._mul_c(<RingElement>right)
1039
1040        if not (PY_TYPE_CHECK(self, Element) and PY_TYPE_CHECK(right, Element)):
1041            # one of self or right is not even an Element.
1042            return bin_op_c(self, right, operator.mul)
1043
1044        # Now we can assume both self and right are of a class that derives
1045        # from Element (so they have a parent).  If one is a ModuleElement,
1046        # do some special code.
1047        if PY_TYPE_CHECK(self, ModuleElement) and PY_TYPE_CHECK(right, ModuleElement):
1048            # We may assume both are module elements.
1049            if (<Element>self)._parent is (<Element>right)._parent._base:
1050                return (<ModuleElement>right)._rmul_c(self)
1051            elif (<Element>self)._parent._base is (<Element>right)._parent:
1052                return (<ModuleElement>self)._lmul_c(right)
1053            if not PY_TYPE_CHECK(right, RingElement):
1054                # Now self must be a ring element:
1055                # If the parent is the same as the base ring, good
1056                if (<RingElement>self)._parent is (<ModuleElement>right)._parent._base:
1057                    return (<ModuleElement>right)._rmul_c(self)
1058                elif PY_TYPE_CHECK(right, Matrix):
1059                    return (<Matrix>right)._rmultiply_by_scalar(left)
1060                else:
1061                    # Otherwise we have to do an explicit canonical coercion.
1062                    try:
1063                        return (<ModuleElement>right)._rmul_c(
1064                            (<Parent>(<ModuleElement>right)._parent._base)._coerce_c(self))
1065                    except TypeError:
1066                        # that failed -- try to base extend right then do the multiply:
1067                        right = right.base_extend((<RingElement>self)._parent)
1068                        return (<ModuleElement>right)._rmul_c(self)
1069            elif PY_TYPE_CHECK(right, Matrix):  # matrix is a ring element
1070                return (<Matrix>right)._rmultiply_by_scalar(self)
1071
1072        # General case.
1073        return bin_op_c(self, right, operator.mul)
1074
1075    cdef RingElement _mul_c(self, RingElement right):
1076        """
1077        Multiplication dispatcher for RingElements.
1078        DO NOT OVERRIDE THIS FUNCTION.
1079        See extensive documentation at the top of element.pyx.
1080        """
1081        if HAS_DICTIONARY(self):   # fast check
1082            return self._mul_(right)
1083        else:
1084            return self._mul_c_impl(right)
1085
1086    cdef RingElement _mul_c_impl(self, RingElement right):
1087        """
1088        Pyrex classes should override this function to implement multiplication.
1089        DO NOT CALL THIS FUNCTION DIRECTLY.
1090        See extensive documentation at the top of element.pyx.
1091        """
1092        raise TypeError, arith_error_message(self, right, operator.mul)
1093
1094    def _mul_(RingElement self, RingElement right):
1095        """
1096        Python classes should override this function to implement multiplication.
1097        See extensive documentation at the top of element.pyx.
1098        """
1099        return self._mul_c_impl(right)
1100
1101    def __pow__(self, m, dummy):
1102        """
1103        Retern the (integral) power of self.
1104
1105        EXAMPLE:
1106            sage: a = Integers(389)['x']['y'](37)
1107            sage: a^2
1108            202
1109            sage: a^388
1110            1
1111            sage: a^(2^120)
1112            81
1113            sage: a^0
1114            1
1115            sage: a^1 == a
1116            True
1117            sage: a^2 * a^3 == a^5
1118            True
1119            sage: (a^3)^2 == a^6
1120            True
1121            sage: a^57 * a^43 == a^100
1122            True
1123            sage: a^(-1) == 1/a
1124            True
1125            sage: a^200 * a^(-64) == a^136
1126            True
1127        """
1128        cdef int cn
1129
1130        n = int(m)
1131        if n != m:
1132            raise ValueError, "n must be an integer"
1133
1134        if n < 0:
1135            n = -n
1136            a = ~self
1137        else:
1138            a = self
1139
1140        if n < 4:
1141            # These cases will probably be called often
1142            # and don't benifit from the code below
1143            cn = n
1144            if cn == 0:
1145                return (<Element>a)._parent(1)
1146            elif cn == 1:
1147                return a
1148            elif cn == 2:
1149                return a*a
1150            elif cn == 3:
1151                return a*a*a
1152
1153        # One multiplication can be saved by starting with
1154        # the smallest power needed rather than with 1
1155        apow = a
1156        while n&1 == 0:
1157            apow = apow*apow
1158            n = n >> 1
1159        power = apow
1160        n = n >> 1
1161
1162        while n != 0:
1163            apow = apow*apow
1164            if n&1 != 0: power = power*apow
1165            n = n >> 1
1166
1167        return power
1168
1169
1170    ##################################
1171    # Division
1172    ##################################
1173
1174    def __truediv__(self, right):
1175        # in sage all divs are true
1176        if not PY_TYPE_CHECK(self, Element):
1177            return bin_op_c(self, right, operator.div)
1178        return self.__div__(right)
1179
1180    def __div__(self, right):
1181        """
1182        Top-level multiplication operator for ring elements.
1183        See extensive documentation at the top of element.pyx.
1184        """
1185        if have_same_parent(self, right):
1186            return (<RingElement>self)._div_c(<RingElement>right)
1187        return bin_op_c(self, right, operator.div)
1188
1189
1190
1191    cdef RingElement _div_c(self, RingElement right):
1192        """
1193        Multiplication dispatcher for RingElements.
1194        DO NOT OVERRIDE THIS FUNCTION.
1195        See extensive documentation at the top of element.pyx.
1196        """
1197        if HAS_DICTIONARY(self):   # fast check
1198            return self._div_(right)
1199        else:
1200            return self._div_c_impl(right)
1201
1202    cdef RingElement _div_c_impl(self, RingElement right):
1203        """
1204        Pyrex classes should override this function to implement division.
1205        DO NOT CALL THIS FUNCTION DIRECTLY.
1206        See extensive documentation at the top of element.pyx.
1207        """
1208        try:
1209            return self._parent.fraction_field()(self, right)
1210        except AttributeError:
1211            raise TypeError, arith_error_message(self, right, operator.div)
1212
1213    def _div_(RingElement self, RingElement right):
1214        """
1215        Python classes should override this function to implement division.
1216        """
1217        return self._div_c_impl(right)
1218
1219    def __pos__(self):
1220        return self
1221
1222    def __invert__(self):
1223        return 1/self
1224
1225    ##################################################
1226
1227    def order(self):
1228        """
1229        Return the additive order of self.
1230        """
1232
1234        """
1235        Return the additive order of self.
1236        """
1237        raise NotImplementedError
1238
1239    def multiplicative_order(self):
1240        r"""
1241        Return the multiplicative order of self, if self is a unit, or raise
1242        \code{ArithmeticError} otherwise.
1243        """
1244        if not self.is_unit():
1245            raise ArithmeticError, "self (=%s) must be a unit to have a multiplicative order."
1246        raise NotImplementedError
1247
1248    def is_unit(self):
1249        if self == 1 or self == -1:
1250            return True
1251        raise NotImplementedError
1252
1253    def is_nilpotent(self):
1254        """
1255        Return True if self is nilpotent, i.e., some power of self
1256        is 0.
1257        """
1258        if self.is_unit():
1259            return False
1260        if self.is_zero():
1261            return True
1262        raise NotImplementedError
1263
1264
1265
1266
1267def is_CommutativeRingElement(x):
1268    """
1269    Return True if x is of type CommutativeRingElement.
1270    """
1271    return IS_INSTANCE(x, CommutativeRingElement)
1272
1273cdef class CommutativeRingElement(RingElement):
1274    def _im_gens_(self, codomain, im_gens):
1275        if len(im_gens) == 1 and self.parent().gen(0) == 1:
1276            return codomain(self)
1277        raise NotImplementedError
1278
1279    def inverse_mod(self, I):
1280        r"""
1281        Return an inverse of self modulo the ideal $I$, if defined,
1282        i.e., if $I$ and self together generate the unit ideal.
1283        """
1284        raise NotImplementedError
1285
1286    def mod(self, I):
1287        r"""
1288        Return a representative for self modulo the ideal I (or the ideal
1289        generated by the elements of I if I is not an ideal.)
1290
1291        EXAMPLE:  Integers
1292        Reduction of 5 modulo an ideal:
1293            sage: n = 5
1294            sage: n.mod(3*ZZ)
1295            2
1296
1297        Reduction of 5 modulo the ideal generated by 3.
1298            sage: n.mod(3)
1299            2
1300
1301        Reduction of 5 modulo the ideal generated by 15 and 6, which is $(3)$.
1302            sage: n.mod([15,6])
1303            2
1304
1305
1306        EXAMPLE: Univiate polynomials
1307            sage: R.<x> = PolynomialRing(QQ)
1308            sage: f = x^3 + x + 1
1309            sage: f.mod(x + 1)
1310            -1
1311
1312        When little is implemented about a given ring, then mod may
1313        return simply return $f$.  For example, reduction is not
1314        implemented for $\Z[x]$ yet. (TODO!)
1315
1316            sage: R.<x> = PolynomialRing(ZZ)
1317            sage: f = x^3 + x + 1
1318            sage: f.mod(x + 1)
1319            x^3 + x + 1
1320
1321
1322
1323        EXAMPLE: Multivariate polynomials
1324        We reduce a polynomial in two variables modulo a polynomial
1325        and an ideal:
1326            sage: R.<x,y,z> = PolynomialRing(QQ, 3)
1327            sage: (x^2 + y^2 + z^2).mod(x+y+z)
1328            2*z^2 + 2*y*z + 2*y^2
1329
1330        Notice above that $x$ is eliminated.  In the next example,
1331        both $y$ and $z$ are eliminated.
1332
1333            sage: (x^2 + y^2 + z^2).mod( (x - y, y - z) )
1334            3*z^2
1335            sage: f = (x^2 + y^2 + z^2)^2; f
1336            z^4 + 2*y^2*z^2 + y^4 + 2*x^2*z^2 + 2*x^2*y^2 + x^4
1337            sage: f.mod( (x - y, y - z) )
1338            9*z^4
1339
1340        In this example $y$ is eliminated.
1341            sage: (x^2 + y^2 + z^2).mod( (x^3, y - z) )
1342            2*z^2 + x^2
1343        """
1344        from sage.rings.all import is_Ideal
1345        if not is_Ideal(I) or not I.ring() is self.parent():
1346            I = self.parent().ideal(I)
1347            #raise TypeError, "I = %s must be an ideal in %s"%(I, self.parent())
1348        return I.reduce(self)
1349
1350cdef class Vector(ModuleElement):
1351    def __mul__(left, right):
1352        if PY_TYPE_CHECK(left, Vector):
1353            # left is the vector
1354            # Possibilities:
1355            #     left * matrix
1356            if PY_TYPE_CHECK(right, Matrix):
1357                return (<Matrix>right)._vector_times_matrix_c(<Vector>left)
1358            #     left * vector
1359            if PY_TYPE_CHECK(right, Vector):
1360                return (<Vector>left)._vector_times_vector_c(<Vector>right)
1361            #     left * scalar
1362            return (<ModuleElement>left)._multiply_by_scalar(right)
1363
1364        else:
1365            # right is the vector
1366            # Possibilities:
1367            #     matrix * right
1368            if PY_TYPE_CHECK(left, Matrix):
1369                return (<Matrix>left)._matrix_times_vector_c(<Vector>right)
1370            #     vector * right
1371            if PY_TYPE_CHECK(left, Vector):
1372                return (<Vector>left)._vector_times_vector_c(<Vector>right)
1373            #     scalar * right
1374            return (<ModuleElement>right)._rmultiply_by_scalar(left)
1375
1376    cdef Vector _vector_times_vector_c(Vector left, Vector right):
1377        if left._degree != right._degree:
1378            raise TypeError, "incompatible degrees"
1379        left, right = canonical_base_coercion_c(left, right)
1380        if HAS_DICTIONARY(left):
1381            return left._vector_times_vector(right)
1382        else:
1383            return left._vector_times_vector_c_impl(right)
1384    cdef Vector _vector_times_vector_c_impl(Vector left, Vector right):
1385        raise TypeError,arith_error_message(left, right, operator.mul)
1386
1387    def  _vector_times_vector(left, right):
1388        return self.vector_time_vector_c_impl(right)
1389
1390    def __div__(self, right):
1391        if PY_TYPE_CHECK(self, Vector):
1392            right = (<Vector>self)._parent._base._coerce_c(right)
1393            return (<Vector>self)._lmul_c(~right)
1394        raise TypeError, arith_error_message(self, right, operator.div)
1395
1396
1397cdef have_same_base(Element x, Element y):
1398    return x._parent._base is y._parent._base
1399
1400
1401def is_Vector(x):
1402    return IS_INSTANCE(x, Vector)
1403
1404cdef class Matrix(AlgebraElement):
1405    cdef int is_sparse_c(self):
1406        raise NotImplementedError
1407
1408    cdef int is_dense_c(self):
1409        raise NotImplementedError
1410
1411    def __mul__(left, right):
1412        if PY_TYPE_CHECK(left, Matrix):
1413            # left is the matrix
1414            # Possibilities:
1415            #     left * matrix
1416            if PY_TYPE_CHECK(right, Matrix):
1417                return (<Matrix>left)._matrix_times_matrix_c(<Vector>right)
1418            #     left * vector
1419            if PY_TYPE_CHECK(right, Vector):
1420                return (<Matrix>left)._matrix_times_vector_c(<Vector>right)
1421            #     left * scalar
1422            return (<Matrix>left)._multiply_by_scalar(right)
1423
1424        else:
1425            # right is the matrix
1426            # Possibilities:
1427            #     matrix * right
1428            if PY_TYPE_CHECK(left, Matrix):
1429                return (<Matrix>left)._matrix_times_matrix_c(<Matrix>right)
1430            #     vector * right
1431            if PY_TYPE_CHECK(left, Vector):
1432                return (<Matrix>right)._vector_times_matrix_c(<Vector>left)
1433            #     scalar * right
1434            return (<Matrix>right)._rmultiply_by_scalar(left)
1435
1436    cdef Vector _vector_times_matrix_c(matrix_right, Vector vector_left):
1437        if vector_left._degree != matrix_right._nrows:
1438            raise TypeError, "incompatible dimensions"
1439        matrix_right, vector_left = canonical_base_coercion_c(matrix_right, vector_left)
1440        if HAS_DICTIONARY(matrix_right):
1441            return matrix_right._vector_times_matrix(vector_left)
1442        else:
1443            return matrix_right._vector_times_matrix_c_impl(vector_left)
1444
1445    cdef Vector _vector_times_matrix_c_impl(matrix_right, Vector vector_left):
1446        raise TypeError
1447
1448    def _vector_times_matrix(matrix_right, vector_left):
1449        return matrix_right._vector_times_matrix_c_impl(vector_left)
1450
1451
1452    cdef Vector _matrix_times_vector_c(matrix_left, Vector vector_right):
1453        if matrix_left._ncols != vector_right._degree:
1454            raise TypeError, "incompatible dimensions"
1455        matrix_left, vector_right = canonical_base_coercion_c(matrix_left, vector_right)
1456        if HAS_DICTIONARY(matrix_left):
1457            return matrix_left._matrix_times_vector(vector_right)
1458        else:
1459            return matrix_left._matrix_times_vector_c_impl(vector_right)
1460
1461    cdef Vector _matrix_times_vector_c_impl(matrix_left, Vector vector_right):
1462        raise TypeError
1463    def _matrix_times_vector(matrix_left, vector_right):
1464        return matrix_left._matrix_times_vector_c_impl(vector_right)
1465
1466
1467    cdef Matrix _matrix_times_matrix_c(left, Matrix right):
1468        cdef int sl, sr
1469        if left._ncols != right._nrows:
1470            raise TypeError, "incompatible dimensions"
1471        left, right = canonical_base_coercion_c(left, right)
1472        sl = left.is_sparse_c(); sr = right.is_sparse_c()
1473        if sl != sr:  # is dense and one is sparse
1474            if sr:  # left is dense
1475                right = right.dense_matrix()
1476            else:
1477                left = left.dense_matrix()
1478        if HAS_DICTIONARY(left):
1479            return left._matrix_times_matrix(right)
1480        else:
1481            return left._matrix_times_matrix_c_impl(right)
1482
1483    cdef Matrix _matrix_times_matrix_c_impl(left, Matrix right):
1484        raise TypeError
1485    def _matrix_time_matrix(left, right):
1486        return left._matrix_times_matrix_c_impl(right)
1487
1488
1489def is_Matrix(x):
1490    return IS_INSTANCE(x, Matrix)
1491
1492def is_IntegralDomainElement(x):
1493    """
1494    Return True if x is of type IntegralDomainElement.
1495    """
1496    return IS_INSTANCE(x, IntegralDomainElement)
1497
1498cdef class IntegralDomainElement(CommutativeRingElement):
1499    def is_nilpotent(self):
1500        return self.is_zero()
1501
1502
1503def is_DedekindDomainElement(x):
1504    """
1505    Return True if x is of type DedekindDomainElement.
1506    """
1507    return IS_INSTANCE(x, DedekindDomainElement)
1508
1509cdef class DedekindDomainElement(IntegralDomainElement):
1510    pass
1511
1512def is_PrincipalIdealDomainElement(x):
1513    """
1514    Return True if x is of type PrincipalIdealDomainElement.
1515    """
1516    return IS_INSTANCE(x, PrincipalIdealDomainElement)
1517
1518cdef class PrincipalIdealDomainElement(DedekindDomainElement):
1519    def lcm(self, right):
1520        """
1521        Returns the least common multiple of self and right.
1522        """
1523        if not PY_TYPE_CHECK(right, Element) or not ((<Element>right)._parent is self._parent):
1524            return bin_op_c(self, right, lcm)
1525        return self._lcm(right)
1526
1527    def gcd(self, right):
1528        """
1529        Returns the gcd of self and right, or 0 if both are 0.
1530        """
1531        if not PY_TYPE_CHECK(right, Element) or not ((<Element>right)._parent is self._parent):
1532            return bin_op_c(self, right, gcd)
1533        return self._gcd(right)
1534
1535    def xgcd(self, right):
1536        r"""
1537        Return the extended gcd of self and other, i.e., elements $r, s, t$ such that
1538        $$1539 r = s \cdot self + t \cdot other. 1540$$
1541        """
1542        if not PY_TYPE_CHECK(right, Element) or not ((<Element>right)._parent is self._parent):
1543            return bin_op_c(self, right, xgcd)
1544        return self._xgcd(right)
1545
1546
1547# This is pretty nasty low level stuff. The idea is to speed up construction
1548# of EuclideanDomainElements (in particular Integers) by skipping some tp_new
1549# calls up the inheritance tree.
1550PY_SET_TP_NEW(EuclideanDomainElement, Element)
1551
1552def is_EuclideanDomainElement(x):
1553    """
1554    Return True if x is of type EuclideanDomainElement.
1555    """
1556    return IS_INSTANCE(x, EuclideanDomainElement)
1557
1558cdef class EuclideanDomainElement(PrincipalIdealDomainElement):
1559
1560    def degree(self):
1561        raise NotImplementedError
1562
1563    def _gcd(self, other):
1564        """
1565        Return the greatest common divisor of self and other.
1566
1567        Algorithm 3.2.1 in Cohen, GTM 138.
1568        """
1569        A = self
1570        B = other
1571        while not B.is_zero():
1572            Q, R = A.quo_rem(B)
1573            A = B
1574            B = R
1575        return A
1576
1578        raise NotImplementedError
1579
1580    def quo_rem(self, other):
1581        raise NotImplementedError
1582
1583    def __floordiv__(self,right):
1584        """
1585        Quotient of division of self by other.  This is denoted //.
1586        """
1587        Q, _ = self.quo_rem(right)
1588        return Q
1589
1590    def __mod__(self, other):
1591        """
1592        Remainder of division of self by other.
1593
1594        EXAMPLES:
1595            sage: R.<x> = ZZ[]
1596            sage: x % (x+1)
1597            -1
1598            sage: (x**3 + x - 1) % (x**2 - 1)
1599            2*x - 1
1600        """
1601        _, R = self.quo_rem(other)
1602        return R
1603
1604def is_FieldElement(x):
1605    """
1606    Return True if x is of type FieldElement.
1607    """
1608    return IS_INSTANCE(x, FieldElement)
1609
1610cdef class FieldElement(CommutativeRingElement):
1611
1612    def is_unit(self):
1613        """
1614        Return True if self is a unit in its parent ring.
1615
1616        EXAMPLES:
1617            sage: a = 2/3; a.is_unit()
1618            True
1619
1620        On the other hand, 2 is not a unit, since its parent is ZZ.
1621            sage: a = 2; a.is_unit()
1622            False
1623            sage: parent(a)
1624            Integer Ring
1625
1626        However, a is a unit when viewed as an element of QQ:
1627            sage: a = QQ(2); a.is_unit()
1628            True
1629        """
1630        return PyBool_FromLong(not self.is_zero())
1631
1632    def _gcd(self, FieldElement other):
1633        """
1634        Return the greatest common divisor of self and other.
1635        """
1636        if self.is_zero() and other.is_zero():
1637            return self
1638        else:
1639            return self.parent()(1)
1640
1641    def _lcm(self, FieldElement other):
1642        """
1643        Return the least common multiple of self and other.
1644        """
1645        if self.is_zero() and other.is_zero():
1646            return self
1647        else:
1648            return self.parent()(1)
1649
1650    def _xgcd(self, FieldElement other):
1651        R = self.parent()
1652        if not self.is_zero():
1653            return R(1), ~self, R(0)
1654        elif not other.is_zero():
1655            return R(1), R(0), ~self
1656        else: # both are 0
1657            return self, self, self
1658
1659
1660    def quo_rem(self, right):
1661        if not isinstance(right, FieldElement) or not (right.parent() is self.parent()):
1662            right = self.parent()(right)
1663        return self/right, 0
1664
1665## def is_FiniteFieldElement(x):
1666##     """
1667##     Return True if x is of type FiniteFieldElement.
1668##     """
1669##     return IS_INSTANCE(x, FiniteFieldElement)
1670
1671cdef class FiniteFieldElement(FieldElement):
1672    pass
1673
1674def is_AlgebraElement(x):
1675    """
1676    Return True if x is of type AlgebraElement.
1677    """
1678    return IS_INSTANCE(x, AlgebraElement)
1679
1680cdef class AlgebraElement(RingElement):
1681    pass
1682
1683def is_CommutativeAlgebraElement(x):
1684    """
1685    Return True if x is of type CommutativeAlgebraElement.
1686    """
1687    return IS_INSTANCE(x, CommutativeAlgebraElement)
1688
1689cdef class CommutativeAlgebraElement(CommutativeRingElement):
1690    pass
1691
1692def is_InfinityElement(x):
1693    """
1694    Return True if x is of type InfinityElement.
1695    """
1696    return IS_INSTANCE(x, InfinityElement)
1697
1698cdef class InfinityElement(RingElement):
1699    pass
1700
1701cdef class PlusInfinityElement(InfinityElement):
1702    pass
1703
1704cdef class MinusInfinityElement(InfinityElement):
1705    pass
1706
1707
1708cdef int have_same_parent(left, right):
1709    """
1710    Return nonzero true value if and only if left and right are
1711    elements and have the same parent.
1712    """
1713    # (We know at least one of the arguments is an Element. So if
1714    # their types are *equal* (fast to check) then they are both
1715    # Elements.  Otherwise use the slower test via PY_TYPE_CHECK.)
1716    if PY_TYPE(left) is PY_TYPE(right):
1717        return (<Element>left)._parent is (<Element>right)._parent
1718
1719    if PY_TYPE_CHECK(right, Element) and PY_TYPE_CHECK(left, Element):
1720        return (<Element>left)._parent is (<Element>right)._parent
1721
1722    return 0
1723
1724
1725
1726
1727
1728
1729#################################################################################
1730#
1731#  Coercion of elements
1732#
1733#################################################################################
1734import __builtin__
1735import operator
1736
1737cimport sage.modules.module
1738import  sage.modules.module
1739
1740#################################################################################
1741# parent
1742#################################################################################
1743cdef parent_c(x):
1744    if PY_TYPE_CHECK(x,Element):
1745        return (<Element>x)._parent
1746    return <object>PY_TYPE(x)
1747
1748def parent(x):
1749    return parent_c(x)
1750
1751#################################################################################
1752# coerce
1753#################################################################################
1754def coerce(Parent p, x):
1755    try:
1756        return p._coerce_c(x)
1757    except AttributeError:
1758        return p(x)
1759
1760
1761#################################################################################
1762# canonical coercion of two ring elements into one of their parents.
1763#################################################################################
1764cdef _verify_canonical_coercion_c(x, y):
1765    if not have_same_parent(x,y):
1766        raise RuntimeError, """There is a bug in the coercion code in SAGE.
1767Both x (=%s) and y (=%s) are supposed to have identical parents but they don't.
1768In fact, x has parent '%s'
1769whereas y has parent '%s'"""%(x,y,parent_c(x),parent_c(y))
1770    return x, y
1771
1772def canonical_coercion(x, y):
1773    """
1774    canonical_coercion(x,y) is what is called before doing an
1775    arithmetic operation between x and y.  It returns a pair (z,w)
1776    such that z is got from x and w from y via canonical coercion and
1777    the parents of z and w are identical.
1778
1779    EXAMPLES:
1780        sage: A = Matrix([[0,1],[1,0]])
1781        sage: canonical_coercion(A,1)
1782        ([0 1]
1783        [1 0], [1 0]
1784        [0 1])
1785    """
1786    return canonical_coercion_c(x,y)
1787
1788cdef canonical_coercion_c(x, y):
1789    cdef int i
1790    xp = parent_c(x)
1791    yp = parent_c(y)
1792    if xp is yp:
1793        return x, y
1794
1795    if PY_IS_NUMERIC(x):
1796        try:
1797            x = yp(x)
1798        except TypeError:
1799            y = x.__class__(y)
1800            return x, y
1801        # Calling this every time incurs overhead -- however, if a mistake
1802        # gets through then one can get infinite loops in C code hence core
1803        # dumps.  And users define _coerce_ and __call__ for rings, which
1804        # can easily have bugs in it, i.e., not really make the element
1805        # have the correct parent.  Thus this check is *crucial*.
1806        return _verify_canonical_coercion_c(x,y)
1807
1808    elif PY_IS_NUMERIC(y):
1809        try:
1810            y = xp(y)
1811        except TypeError:
1812            x = y.__class__(x)
1813            return x, y
1814        return _verify_canonical_coercion_c(x,y)
1815
1816    try:
1817        if xp.has_coerce_map_from(yp):
1818            y = (<Parent>xp)._coerce_c(y)
1819            return _verify_canonical_coercion_c(x,y)
1820    except AttributeError:
1821        pass
1822    try:
1823        if yp.has_coerce_map_from(xp):
1824            x = (<Parent>yp)._coerce_c(x)
1825            return _verify_canonical_coercion_c(x,y)
1826    except AttributeError:
1827        pass
1828    raise TypeError, "no common canonical parent for objects with parents: '%s' and '%s'"%(xp, yp)
1829
1830cdef canonical_base_coercion_c(Element x, Element y):
1831    if not have_same_base(x, y):
1832        if (<Parent> x._parent._base).has_coerce_map_from_c(y._parent._base):
1833            # coerce all elements of y to the base ring of x
1834            y = y.base_extend_c(x._parent._base)
1835        elif (<Parent> y._parent._base).has_coerce_map_from_c(x._parent._base):
1836            # coerce x to have elements in the base ring of y
1837            x = x.base_extend_c(y._parent._base)
1838    return x,y
1839
1840def canonical_base_coercion(x, y):
1841    try:
1842        xb = x.base_ring()
1843    except AttributeError:
1844        #raise TypeError, "unable to find base ring for %s (parent: %s)"%(x,x.parent())
1845        raise TypeError, "unable to find base ring"
1846    try:
1847        yb = y.base_ring()
1848    except AttributeError:
1849        raise TypeError, "unable to find base ring"
1850        #raise TypeError, "unable to find base ring for %s (parent: %s)"%(y,y.parent())
1851    try:
1852        b = canonical_coercion_c(xb(0),yb(0))[0].parent()
1853    except TypeError:
1854        raise TypeError, "unable to find base ring"
1855        #raise TypeError, "unable to find a common base ring for %s (base ring: %s) and %s (base ring %s)"%(x,xb,y,yb)
1856    return x.change_ring(b), y.change_ring(b)
1857
1858
1859D = {'mul':'*', 'add':'+', 'sub':'-', 'div':'/'}
1860cdef arith_error_message(x, y, op):
1861    try:
1862        n = D[op.__name__]
1863    except KeyError:
1864        n = op.__name__
1865    return "unsupported operand parent(s) for '%s': '%s' and '%s'"%(n, parent_c(x), parent_c(y))
1866
1867def bin_op(x, y, op):
1868    return bin_op_c(x,y,op)
1869
1870cdef bin_op_c(x, y, op):
1871    """
1872    Compute x op y, where coercion of x and y works according to
1873    SAGE's coercion rules.
1874    """
1875    # Try canonical element coercion.
1876    try:
1877        x1, y1 = canonical_coercion_c(x, y)
1878        return op(x1,y1)
1879    except TypeError, msg:
1880        # print msg  # this can be useful for debugging.
1881        if not op is operator.mul:
1882            raise TypeError, arith_error_message(x,y,op)
1883
1884    # If the op is multiplication, then some other algebra multiplications
1885    # may be defined
1886
1887    # 2. Try scalar multiplication.
1888    # No way to multiply x and y using the coerce into a canonical
1889    # parent'' rule.
1890    # The next rule to try is scalar multiplication by coercing
1891    # into the base ring.
1892    cdef int x_is_modelt, y_is_modelt
1893
1894    y_is_modelt = PY_TYPE_CHECK(y, ModuleElement)
1895    if y_is_modelt:
1896        # First try to coerce x into the base ring of y if y is an element.
1897        try:
1898            R = (<ModuleElement> y)._parent._base
1899            if R is None:
1900                raise RuntimeError, "base of '%s' must be set to a ring (but it is None)!"%((<ModuleElement> y)._parent)
1901            x = (<Parent>R)._coerce_c(x)
1902            return (<ModuleElement> y)._rmul_c(x)     # the product x * y
1903        except TypeError, msg:
1904            pass
1905
1906    x_is_modelt = PY_TYPE_CHECK(x, ModuleElement)
1907    if x_is_modelt:
1908        # That did not work.  Try to coerce y into the base ring of x.
1909        try:
1910            R = (<ModuleElement> x)._parent._base
1911            if R is None:
1912                raise RuntimeError, "base of '%s' must be set to a ring (but it is None)!"%((<ModuleElement> x)._parent)
1913            y = (<Parent> R)._coerce_c(y)
1914            return (<ModuleElement> x)._lmul_c(y)    # the product x * y
1915        except TypeError:
1916            pass
1917
1918    if y_is_modelt and x_is_modelt:
1919        # 3. Both canonical coercion failed, but both are module elements.
1920        # Try base extending the right object by the parent of the left
1921
1922        ## TODO -- WORRY -- only unambiguous if one succeeds!
1923        if  PY_TYPE_CHECK(x, RingElement):
1924            try:
1925                return x * y.base_extend((<RingElement>x)._parent)
1926            except (TypeError, AttributeError), msg:
1927                pass
1928        # Also try to base extending the left object by the parent of the right
1929        if  PY_TYPE_CHECK(y, RingElement):
1930            try:
1931                return y * x.base_extend((<Element>y)._parent)
1932            except (TypeError, AttributeError), msg:
1933                pass
1934
1935    # 4. Try _l_action or _r_action.
1936    # Test to see if an _r_action or _l_action is
1937    # defined on either side.
1938    try:
1939        return x._l_action(y)
1940    except (AttributeError, TypeError):
1941        pass
1942    try:
1943        return y._r_action(x)
1944    except (AttributeError, TypeError):
1945        pass
1946
1947    raise TypeError, arith_error_message(x,y,op)
1948
1949def coerce_cmp(x,y):
1950    cdef int c
1951    try:
1952        x, y = canonical_coercion_c(x, y)
1953        return cmp(x,y)
1954    except TypeError:
1955        c = cmp(type(x), type(y))
1956        if c == 0: c = -1
1957        return c
1958
1959
1960
1961###############################################################################
1962
1963def lcm(x,y):
1964    from sage.rings.arith import lcm
1965    return lcm(x,y)
1966
1967def gcd(x,y):
1968    from sage.rings.arith import gcd
1969    return gcd(x,y)
1970
1971def xgcd(x,y):
1972    from sage.rings.arith import xgcd
1973    return xgcd(x,y)
1974
1975
1976
1977
1978######################
1979
1980def generic_power(m, dummy):
1981    return generic_power_c(m, dummy)
1982
1983cdef generic_power_c(m, dummy):
1984    cdef int cn
1985
1986    n = int(m)
1987    if n != m:
1988        raise ValueError, "n must be an integer"
1989
1990    if n < 0:
1991        n = -n
1992        a = ~self
1993    else:
1994        a = self
1995
1996    if n < 4:
1997        # These cases will probably be called often
1998        # and don't benifit from the code below
1999        cn = n
2000        if cn == 0:
2001            return (<Element>a)._parent(1)
2002        elif cn == 1:
2003            return a
2004        elif cn == 2:
2005            return a*a
2006        elif cn == 3:
2007            return a*a*a
2008
2009    # One multiplication can be saved by starting with
2010    # the smallest power needed rather than with 1
2011    apow = a
2012    while n&1 == 0:
2013        apow = apow*apow
2014        n = n >> 1
2015    power = apow
2016    n = n >> 1
2017
2018    while n != 0:
2019        apow = apow*apow
2020        if n&1 != 0: power = power*apow
2021        n = n >> 1
2022
2023    return power
2024
Note: See TracBrowser for help on using the repository browser.