# source:sage/rings/integer.pyx@6031:18e45aa21972

Revision 6031:18e45aa21972, 88.5 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

2.8.3.rc1

Line
1r"""
2Elements of the ring $\Z$ of integers
3
4AUTHORS:
5    -- William Stein (2005): initial version
6    -- Gonzalo Tornaria (2006-03-02): vastly improved python/GMP conversion; hashing
7    -- Didier Deshommes <dfdeshom@gmail.com> (2006-03-06): numerous examples and docstrings
8    -- William Stein (2006-03-31): changes to reflect GMP bug fixes
9    -- William Stein (2006-04-14): added GMP factorial method (since it's
10                                   now very fast).
11    -- David Harvey (2006-09-15): added nth_root, exact_log
12    -- David Harvey (2006-09-16): attempt to optimise Integer constructor
13    -- Rishikesh (2007-02-25): changed quo_rem so that the rem is positive
14    -- David Harvey, Martin Albrecht, Robert Bradshaw (2007-03-01): optimized Integer constructor and pool
15    -- Pablo De Napoli (2007-04-01): multiplicative_order should return +infinity for non zero numbers
16    -- Robert Bradshaw (2007-04-12): is_perfect_power, Jacobi symbol (with Kronecker extension)
17                                     Convert some methods to use GMP directly rather than pari, Integer() -> PY_NEW(Integer)
18    -- David Roe (2007-03-21): sped up valuation and is_square, added val_unit, is_power, is_power_of and divide_knowing_divisible_by
19
20EXAMPLES:
21   Add 2 integers:
22       sage: a = Integer(3) ; b = Integer(4)
23       sage: a + b == 7
24       True
25
26   Add an integer and a real number:
27       sage: a + 4.0
28       7.00000000000000
29
30   Add an integer and a rational number:
31       sage: a + Rational(2)/5
32       17/5
33
34   Add an integer and a complex number:
35       sage: b = ComplexField().0 + 1.5
36       sage: loads((a+b).dumps()) == a+b
37       True
38
39   sage: z = 32
40   sage: -z
41   -32
42   sage: z = 0; -z
43   0
44   sage: z = -0; -z
45   0
46   sage: z = -1; -z
47   1
48
49Multiplication:
50    sage: a = Integer(3) ; b = Integer(4)
51    sage: a * b == 12
52    True
53    sage: loads((a * 4.0).dumps()) == a*b
54    True
55    sage: a * Rational(2)/5
56    6/5
57
58    sage: list([2,3]) * 4
59    [2, 3, 2, 3, 2, 3, 2, 3]
60
61    sage: 'sage'*Integer(3)
62    'sagesagesage'
63
64Coercions:
65    Returns version of this integer in the multi-precision floating
66    real field R.
67
68        sage: n = 9390823
69        sage: RR = RealField(200)
70        sage: RR(n)
71        9390823.0000000000000000000000000000000000000000000000000000
72
73"""
74
75#*****************************************************************************
76#       Copyright (C) 2004 William Stein <wstein@gmail.com>
77#
78#  Distributed under the terms of the GNU General Public License (GPL)
79#
80#    This code is distributed in the hope that it will be useful,
81#    but WITHOUT ANY WARRANTY; without even the implied warranty of
82#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
83#    General Public License for more details.
84#
85#  The full text of the GPL is available at:
86#
88#*****************************************************************************
89
90doc="""
91Integers
92"""
93
94import operator
95
96import sys
97
98include "../ext/gmp.pxi"
99include "../ext/interrupt.pxi"  # ctrl-c interrupt block support
100include "../ext/stdsage.pxi"
101include "../ext/python_list.pxi"
102
103cdef extern from "mpz_pylong.h":
104    cdef mpz_get_pylong(mpz_t src)
105    cdef mpz_get_pyintlong(mpz_t src)
106    cdef int mpz_set_pylong(mpz_t dst, src) except -1
107    cdef long mpz_pythonhash(mpz_t src)
108
109
110from sage.libs.pari.gen cimport gen as pari_gen
111
112cdef class Integer(sage.structure.element.EuclideanDomainElement)
113
114import sage.rings.infinity
115import sage.libs.pari.all
116
117cdef mpz_t mpz_tmp
118mpz_init(mpz_tmp)
119
120cdef public int set_mpz(Integer self, mpz_t value):
121    mpz_set(self.value, value)
122
123cdef set_from_Integer(Integer self, Integer other):
124    mpz_set(self.value, other.value)
125
126cdef set_from_int(Integer self, int other):
127    mpz_set_si(self.value, other)
128
129cdef public mpz_t* get_value(Integer self):
130    return &self.value
131
132MAX_UNSIGNED_LONG = 2 * sys.maxint
133
134# This crashes SAGE:
135#  s = 2003^100300000
136# The problem is related to realloc moving all the memory
137# and returning a pointer to the new block of memory, I think.
138
139from sage.structure.sage_object cimport SageObject
140from sage.structure.element cimport EuclideanDomainElement, ModuleElement
141from sage.structure.element import  bin_op
142
143import integer_ring
144the_integer_ring = integer_ring.ZZ
145
146cdef set_zero_one_elements():
147    global the_integer_ring
148    the_integer_ring._zero_element = Integer(0)
149    the_integer_ring._one_element = Integer(1)
150
151set_zero_one_elements()
152
153def is_Integer(x):
154    return PY_TYPE_CHECK(x, Integer)
155
156cdef class Integer(sage.structure.element.EuclideanDomainElement):
157    r"""
158    The \class{Integer} class represents arbitrary precision
159    integers.  It derives from the \class{Element} class, so
160    integers can be used as ring elements anywhere in SAGE.
161
162    \begin{notice}
163    The class \class{Integer} is implemented in Pyrex, as a wrapper
164    of the GMP \code{mpz_t} integer type.
165    \end{notice}
166    """
167
168    # todo: It would be really nice if we could avoid the __new__ call.
169    # It has python calling conventions, and our timing tests indicate the
170    # overhead can be significant. The difficulty is that then we can't
171    # guarantee that the initialization will be performed exactly once.
172
173    def __new__(self, x=None, unsigned int base=0):
174        mpz_init(self.value)
175        self._parent = <SageObject>the_integer_ring
176
177    def __pyxdoc__init__(self):
178        """
179        You can create an integer from an int, long, string literal, or
180        integer modulo N.
181
182        EXAMPLES:
183            sage: Integer(495)
184            495
185            sage: Integer('495949209809328523')
186            495949209809328523
187            sage: Integer(Mod(3,7))
188            3
189            sage: 2^3
190            8
191        """
192
193    def __init__(self, x=None, unsigned int base=0):
194        """
195        EXAMPLES:
196            sage: a = long(-901824309821093821093812093810928309183091832091)
197            sage: b = ZZ(a); b
198            -901824309821093821093812093810928309183091832091
199            sage: ZZ(b)
200            -901824309821093821093812093810928309183091832091
201            sage: ZZ('-901824309821093821093812093810928309183091832091')
202            -901824309821093821093812093810928309183091832091
203            sage: ZZ(int(-93820984323))
204            -93820984323
205            sage: ZZ(ZZ(-901824309821093821093812093810928309183091832091))
206            -901824309821093821093812093810928309183091832091
207            sage: ZZ(QQ(-901824309821093821093812093810928309183091832091))
208            -901824309821093821093812093810928309183091832091
209            sage: ZZ(pari('Mod(-3,7)'))
210            4
211            sage: ZZ('sage')
212            Traceback (most recent call last):
213            ...
214            TypeError: unable to convert x (=sage) to an integer
215            sage: Integer('zz',36).str(36)
216            'zz'
217            sage: ZZ('0x3b').str(16)
218            '3b'
219            sage: ZZ( ZZ(5).digits(3) , 3)
220            5
221        """
222
223        # TODO: All the code below should somehow be in an external
224        # cdef'd function.  Then e.g., if a matrix or vector or
225        # polynomial is getting filled by mpz_t's, it can use the
226        # rules below to do the fill construction of mpz_t's, but
227        # without the overhead of creating any Python objects at all.
228        # The cdef's function should be of the form
229        #     mpz_init_set_sage(mpz_t y, object x)
230        # Then this function becomes the one liner:
231        #     mpz_init_set_sage(self.value, x)
232
233        cdef Integer tmp
234
235        if x is None:
236            if mpz_sgn(self.value) != 0:
237                mpz_set_si(self.value, 0)
238
239        else:
240            # First do all the type-check versions; these are fast.
241
242            if PY_TYPE_CHECK(x, Integer):
243                set_from_Integer(self, <Integer>x)
244
245            elif PyInt_Check(x):
246                mpz_set_si(self.value, x)
247
248            elif PyLong_Check(x):
249                mpz_set_pylong(self.value, x)
250
251            elif PyString_Check(x):
252                if base < 0 or base > 36:
253                    raise ValueError, "base (=%s) must be between 2 and 36"%base
254                if mpz_set_str(self.value, x, base) != 0:
255                    raise TypeError, "unable to convert x (=%s) to an integer"%x
256
257            # Similarly for "sage.libs.pari.all.pari_gen"
258            elif PY_TYPE_CHECK(x, pari_gen):
259                if x.type() == 't_INTMOD':
260                    x = x.lift()
261                # TODO: figure out how to convert to pari integer in base 16 ?
262
263                # todo: having this "s" variable around here is causing
264                # pyrex to play games with refcount for the None object, which
265                # seems really stupid.
266
267                s = hex(x)
268                if mpz_set_str(self.value, s, 16) != 0:
269                    raise TypeError, "Unable to coerce PARI %s to an Integer."%x
270            elif PyObject_HasAttrString(x, "_integer_"):
271                # todo: Note that PyObject_GetAttrString returns NULL if
272                # the attribute was not found. If we could test for this,
273                # we could skip the double lookup. Unfortunately pyrex doesn't
274                # seem to let us do this; it flags an error if the function
275                # returns NULL, because it can't construct an "object" object
276                # out of the NULL pointer. This really sucks. Perhaps we could
277                # make the function prototype have return type void*, but
278                # then how do we make Pyrex handle the reference counting?
279                set_from_Integer(self, (<object> PyObject_GetAttrString(x, "_integer_"))())
280
281            elif PY_TYPE_CHECK(x, list) and base > 1:
282                b = the_integer_ring(base)
283                tmp = the_integer_ring(0)
284                for i in range(len(x)):
285                    tmp += x[i]*b**i
286                mpz_set(self.value, tmp.value)
287
288            else:
289                raise TypeError, "unable to coerce element to an integer"
290
291
292    def __reduce__(self):
293        # This single line below took me HOURS to figure out.
294        # It is the *trick* needed to pickle pyrex extension types.
295        # The trick is that you must put a pure Python function
296        # as the first argument, and that function must return
297        # the result of unpickling with the argument in the second
298        # tuple as input. All kinds of problems happen
299        # if we don't do this.
300        return sage.rings.integer.make_integer, (self.str(32),)
301
302    def _reduce_set(self, s):
303        mpz_set_str(self.value, s, 32)
304
305    def __index__(self):
306        """
307        Needed so integers can be used as list indices.
308
309        EXAMPLES:
310            sage: v = [1,2,3,4,5]
311            sage: v[Integer(3)]
312            4
313            sage: v[Integer(2):Integer(4)]
314            [3, 4]
315        """
316        return mpz_get_pyintlong(self.value)
317
318    def _im_gens_(self, codomain, im_gens):
319        return codomain._coerce_(self)
320
321    def _xor(Integer self, Integer other):
322        cdef Integer x
323        x = PY_NEW(Integer)
324        mpz_xor(x.value, self.value, other.value)
325        return x
326
327    def __xor__(x, y):
328        """
329        Compute the exclusive or of x and y.
330
331        EXAMPLES:
332            sage: n = ZZ(2); m = ZZ(3)
333            sage: n.__xor__(m)
334            1
335        """
336        if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
337            return x._xor(y)
338        return bin_op(x, y, operator.xor)
339
340
341    def __richcmp__(left, right, int op):
342        return (<sage.structure.element.Element>left)._richcmp(right, op)
343
344    cdef int _cmp_c_impl(left, sage.structure.element.Element right) except -2:
345        cdef int i
346        i = mpz_cmp((<Integer>left).value, (<Integer>right).value)
347        if i < 0: return -1
348        elif i == 0: return 0
349        else: return 1
350
351    def __copy__(self):
352        """
353        Return a copy of the integer.
354
355        EXAMPLES:
356            sage: n = 2
357            sage: copy(n)
358            2
359            sage: copy(n) is n
360            False
361        """
362        cdef Integer z
363        z = PY_NEW(Integer)
364        set_mpz(z,self.value)
365        return z
366
367
368    def list(self):
369        """
370        Return a list with this integer in it, to be
371        compatible with the method for number fields.
372
373        EXAMPLES:
374            sage: m = 5
375            sage: m.list()
376            [5]
377        """
378        return [ self ]
379
380
381    def  __dealloc__(self):
382        mpz_clear(self.value)
383
384    def __repr__(self):
385        return self.str()
386
387    def _latex_(self):
388        return self.str()
389
390    def _mathml_(self):
391        return '<mn>%s</mn>'%self
392
393    def __str_malloc(self, int base=10):
394        r"""
395        Return the string representation of \code{self} in the given
396        base.
397
398        However, self.str() below is nice because we know the size of
399        the string ahead of time, and can work around a bug in GMP
400        nicely.  There seems to be a bug in GMP, where non-2-power
401        base conversion for very large integers > 10 million digits
402        (?) crashes GMP.
403        """
404        _sig_on
405        cdef char *s
406        s = mpz_get_str(NULL, base, self.value)
407        t = str(s)
408        free(s)
409        _sig_off
410        return t
411
412    def str(self, int base=10):
413        r"""
414        Return the string representation of \code{self} in the given
415        base.
416
417        EXAMPLES:
418            sage: Integer(2^10).str(2)
419            '10000000000'
420            sage: Integer(2^10).str(17)
421            '394'
422
423            sage: two=Integer(2)
424            sage: two.str(1)
425            Traceback (most recent call last):
426            ...
427            ValueError: base (=1) must be between 2 and 36
428
429            sage: two.str(37)
430            Traceback (most recent call last):
431            ...
432            ValueError: base (=37) must be between 2 and 36
433
434            sage: big = 10^5000000
435            sage: s = big.str()                 # long time (> 20 seconds)
436            sage: len(s)                        # long time (depends on above defn of s)
437            5000001
438            sage: s[:10]                        # long time (depends on above defn of s)
439            '1000000000'
440        """
441        if base < 2 or base > 36:
442            raise ValueError, "base (=%s) must be between 2 and 36"%base
443        cdef size_t n
444        cdef char *s
445        n = mpz_sizeinbase(self.value, base) + 2
446        s = <char *>PyMem_Malloc(n)
447        if s == NULL:
448            raise MemoryError, "Unable to allocate enough memory for the string representation of an integer."
449        _sig_on
450        mpz_get_str(s, base, self.value)
451        _sig_off
452        k = <object> PyString_FromString(s)
453        PyMem_Free(s)
454        return k
455
456    def __hex__(self):
457        r"""
458        Return the hexadecimal digits of self in lower case.
459
460        \note{'0x' is \emph{not} prepended to the result like is done
461        by the corresponding Python function on int or long.  This is
462        for efficiency sake---adding and stripping the string wastes
463        time; since this function is used for conversions from
464        integers to other C-library structures, it is important that
465        it be fast.}
466
467        EXAMPLES:
468            sage: print hex(Integer(15))
469            f
470            sage: print hex(Integer(16))
471            10
472            sage: print hex(Integer(16938402384092843092843098243))
473            36bb1e3929d1a8fe2802f083
474            sage: print hex(long(16938402384092843092843098243))
475            0x36bb1e3929d1a8fe2802f083L
476        """
477        return self.str(16)
478
479    def binary(self):
480        """
481        Return the binary digits of self as a string.
482
483        EXAMPLES:
484            sage: print Integer(15).binary()
485            1111
486            sage: print Integer(16).binary()
487            10000
488            sage: print Integer(16938402384092843092843098243).binary()
489            1101101011101100011110001110010010100111010001101010001111111000101000000000101111000010000011
490        """
491        return self.str(2)
492
493    def digits(self, int base=2, digits=None):
494        """
495        Return a list of digits for self in the given base in little
496        endian order.
497
498        INPUT:
499            base -- integer (default: 2)
500            digits -- optional indexable object as source for the digits
501
502        EXAMPLE:
503            sage: 5.digits()
504            [1, 0, 1]
505            sage: 5.digits(3)
506            [2, 1]
507            sage: 10.digits(16,'0123456789abcdef')
508            ['a']
509
510        """
511        if base < 2:
512            raise TypeError, "base must be >= 2"
513
514        cdef mpz_t mpz_value
515        cdef mpz_t mpz_base
516        cdef mpz_t mpz_res
517        cdef object l = PyList_New(0)
518        cdef Integer z
519        cdef size_t s
520
521        if base == 2 and digits is None:
522            s = mpz_sizeinbase(self.value, 2)
523            o = the_integer_ring._one_element
524            z = the_integer_ring._zero_element
525            for i from 0<= i < s:
526                if mpz_tstbit(self.value,i):
527                    PyList_Append(l,o)
528                else:
529                    PyList_Append(l,z)
530        else:
531            s = mpz_sizeinbase(self.value, 2)
532            if s > 256:
533                _sig_on
534            mpz_init(mpz_value)
535            mpz_init(mpz_base)
536            mpz_init(mpz_res)
537
538            mpz_set(mpz_value, self.value)
539            mpz_set_ui(mpz_base, base)
540
541            if digits:
542                while mpz_cmp_si(mpz_value,0):
543                    mpz_tdiv_qr(mpz_value, mpz_res, mpz_value, mpz_base)
544                    PyList_Append(l, digits[mpz_get_ui(mpz_res)])
545            else:
546                if base < 10:
547                    for e in reversed(self.str(base)):
548                        PyList_Append(l,the_integer_ring(e))
549                else:
550                    while mpz_cmp_si(mpz_value,0):
551                        # TODO: Use a divide and conquer approach
552                        # here: for base b, compute b^2, b^4, b^8,
553                        # ... (repeated squaring) until you get larger
554                        # than your number; then compute (n // b^256,
555                        # n % b ^ 256) (if b^512 > number) to split
556                        # the number in half and recurse
557                        mpz_tdiv_qr(mpz_value, mpz_res, mpz_value, mpz_base)
558                        z = PY_NEW(Integer)
559                        mpz_set(z.value, mpz_res)
560                        PyList_Append(l, z)
561
562            mpz_clear(mpz_value)
563            mpz_clear(mpz_res)
564            mpz_clear(mpz_base)
565
566            if s > 256:
567                _sig_off
568
569        return l # should we return a tuple?
570
571
572    def set_si(self, signed long int n):
573        """
574        Coerces $n$ to a C signed integer if possible, and sets self
575        equal to $n$.
576
577        EXAMPLES:
578            sage: n= ZZ(54)
579            sage: n.set_si(-43344);n
580            -43344
581            sage: n.set_si(43344);n
582            43344
583
584        Note that an error occurs when we are not dealing with
585        integers anymore
586            sage: n.set_si(2^32);n
587            Traceback (most recent call last):      # 32-bit
588            ...                                     # 32-bit
589            OverflowError: long int too large to convert to int   # 32-bit
590            4294967296       # 64-bit
591            sage: n.set_si(-2^32);n
592            Traceback (most recent call last):      # 32-bit
593            ...                                     # 32-bit
594            OverflowError: long int too large to convert to int     # 32-bit
595            -4294967296      # 64-bit
596        """
597        mpz_set_si(self.value, n)
598
599    def set_str(self, s, base=10):
600        """
601        Set self equal to the number defined by the string $s$ in the
602        given base.
603
604        EXAMPLES:
605            sage: n=100
606            sage: n.set_str('100000',2)
607            sage: n
608            32
609
610        If the number begins with '0X' or '0x', it is converted
611        to an hex number:
612            sage: n.set_str('0x13',0)
613            sage: n
614            19
615            sage: n.set_str('0X13',0)
616            sage: n
617            19
618
619        If the number begins with a '0', it is converted to an octal
620        number:
621            sage: n.set_str('013',0)
622            sage: n
623            11
624
625        '13' is not a valid binary number so the following raises
626        an exception:
627            sage: n.set_str('13',2)
628            Traceback (most recent call last):
629            ...
630            TypeError: unable to convert x (=13) to an integer in base 2
631        """
632        valid = mpz_set_str(self.value, s, base)
633        if valid != 0:
634            raise TypeError, "unable to convert x (=%s) to an integer in base %s"%(s, base)
635
636    cdef void set_from_mpz(Integer self, mpz_t value):
637        mpz_set(self.value, value)
638
639    cdef mpz_t* get_value(Integer self):
640        return &self.value
641
642    cdef void _to_ZZ(self, ZZ_c *z):
643        _sig_on
644        mpz_to_ZZ(z, &self.value)
645        _sig_off
646
647    cdef ModuleElement _add_c_impl(self, ModuleElement right):
648        # self and right are guaranteed to be Integers
649        cdef Integer x
650        x = PY_NEW(Integer)
651        mpz_add(x.value, self.value, (<Integer>right).value)
652        return x
653
654##     def _unsafe_add_in_place(self,  ModuleElement right):
655##         """
656##         Do *not* use this...  unless you really know what you
657##         are doing.
658##         """
659##         if not (right._parent is self._parent):
660##             raise TypeError
661##         mpz_add(self.value, self.value, (<Integer>right).value)
662##     cdef _unsafe_add_in_place_c(self,  ModuleElement right):
663##         """
664##         Do *not* use this...  unless you really know what you
665##         are doing.
666##         """
667##         if not (right._parent is self._parent):
668##             raise TypeError
669##         mpz_add(self.value, self.value, (<Integer>right).value)
670
671    cdef ModuleElement _sub_c_impl(self, ModuleElement right):
672        # self and right are guaranteed to be Integers
673        cdef Integer x
674        x = PY_NEW(Integer)
675        mpz_sub(x.value, self.value, (<Integer>right).value)
676        return x
677
678    cdef ModuleElement _neg_c_impl(self):
679        cdef Integer x
680        x = PY_NEW(Integer)
681        mpz_neg(x.value, self.value)
682        return x
683
684    def _r_action(self, s):
685        """
686        EXAMPLES:
687            sage: 8 * [0]
688            [0, 0, 0, 0, 0, 0, 0, 0]
689            sage: 8 * 'hi'
690            'hihihihihihihihi'
691        """
692        if isinstance(s, (str, list, tuple)):
693            return s*int(self)
694        raise TypeError
695
696    def _l_action(self, s):
697        """
698        EXAMPLES:
699            sage: [0] * 8
700            [0, 0, 0, 0, 0, 0, 0, 0]
701            sage: 'hi' * 8
702            'hihihihihihihihi'
703        """
704        if isinstance(s, (str, list, tuple)):
705            return int(self)*s
706        raise TypeError
707
708    cdef RingElement _mul_c_impl(self, RingElement right):
709        # self and right are guaranteed to be Integers
710        cdef Integer x
711        x = PY_NEW(Integer)
712        if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
713            # We only use the signal handler (to enable ctrl-c out) when the
714            # product might take a while to compute
715            _sig_on
716            mpz_mul(x.value, self.value, (<Integer>right).value)
717            _sig_off
718        else:
719            mpz_mul(x.value, self.value, (<Integer>right).value)
720        return x
721
722    cdef RingElement _div_c_impl(self, RingElement right):
723        r"""
724        Computes \frac{a}{b}
725
726        EXAMPLES:
727            sage: a = Integer(3) ; b = Integer(4)
728            sage: a / b == Rational(3) / 4
729            True
730            sage: Integer(32) / Integer(32)
731            1
732        """
733        # this is vastly faster than doing it here, since here
734        # we can't cimport rationals.
735        return the_integer_ring._div(self, right)
736
737    def __floordiv(Integer self, Integer other):
738        cdef Integer x
739        x = PY_NEW(Integer)
740
741        _sig_on
742        mpz_fdiv_q(x.value, self.value, other.value)
743        _sig_off
744
745        return x
746
747
748    def __floordiv__(x, y):
749        r"""
750        Computes the whole part of $\frac{self}{other}$.
751
752        EXAMPLES:
753            sage: a = Integer(321) ; b = Integer(10)
754            sage: a // b
755            32
756        """
757        if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
758            return x.__floordiv(y)
759        return bin_op(x, y, operator.floordiv)
760
761
762    def __pow__(self, n, dummy):
763        r"""
764        Computes $\text{self}^n$
765
766        EXAMPLES:
767            sage: 2^-6
768            1/64
769            sage: 2^6
770            64
771            sage: 2^0
772            1
773            sage: 2^-0
774            1
775            sage: (-1)^(1/3)
776            -1
777
778        The base need not be an integer (it can be a builtin
779        Python type).
780            sage: int(2)^10
781            1024
782            sage: float(2.5)^10
783            9536.7431640625
784            sage: 'sage'^3
785            'sagesagesage'
786
787
788        The exponent must fit in an unsigned long.
789            sage: x = 2^100000000000000000000000
790            Traceback (most recent call last):
791            ...
792            RuntimeError: exponent must be at most 4294967294  # 32-bit
793            RuntimeError: exponent must be at most 18446744073709551614 # 64-bit
794
795        We raise 2 to various interesting exponents:
796            sage: 2^x                # symbolic x
797            2^x
798            sage: 2^1.5              # real number
799            2.82842712474619
800            sage: 2^I                # complex number
801            2^I
802            sage: f = 2^(sin(x)-cos(x)); f
803            2^(sin(x) - cos(x))
804            sage: f(3)
805            2^(sin(3) - cos(3))
806
807        A symbolic sum:
808            sage: x,y,z = var('x,y,z')
809            sage: 2^(x+y+z)
810            2^(z + y + x)
811            sage: 2^(1/2)
812            sqrt(2)
813            sage: 2^(-1/2)
814            1/sqrt(2)
815        """
816        cdef Integer _n
817        cdef unsigned int _nval
818        if not PY_TYPE_CHECK(self, Integer):
819            if isinstance(self, str):
820                return self * n
821            else:
822                return self.__pow__(int(n))
823
824        try:
825            # todo: should add a fast pathway to deal with n being
826            # an Integer or python int
827            _n = Integer(n)
828        except TypeError:
829            try:
830                s = n.parent()(self)
831                return s**n
832            except AttributeError:
833                raise TypeError, "exponent (=%s) must be an integer.\nCoerce your numbers to real or complex numbers first."%n
834
835        if _n < 0:
836            return Integer(1)/(self**(-_n))
837        if _n > MAX_UNSIGNED_LONG:
838            raise RuntimeError, "exponent must be at most %s"%MAX_UNSIGNED_LONG
839        _self = self
840        cdef Integer x
841        x = PY_NEW(Integer)
842        _nval = _n
843
844        _sig_on
845        mpz_pow_ui(x.value, (<Integer>self).value, _nval)
846        _sig_off
847
848        return x
849
850    def nth_root(self, int n, int report_exact=0):
851        r"""
852        Returns the truncated nth root of self.
853
854        INPUT:
855            n -- integer >= 1 (must fit in C int type)
856            report_exact -- boolean, whether to report if the root extraction
857                          was exact
858
859        OUTPUT:
860           If report_exact is 0 (default), then returns the truncation of the
861           nth root of self (i.e. rounded towards zero).
862
863           If report_exact is 1, then returns the nth root and a boolean
864           indicating whether the root extraction was exact.
865
866        AUTHOR:
867           -- David Harvey (2006-09-15)
868
869        EXAMPLES:
870          sage: Integer(125).nth_root(3)
871          5
872          sage: Integer(124).nth_root(3)
873          4
874          sage: Integer(126).nth_root(3)
875          5
876
877          sage: Integer(-125).nth_root(3)
878          -5
879          sage: Integer(-124).nth_root(3)
880          -4
881          sage: Integer(-126).nth_root(3)
882          -5
883
884          sage: Integer(125).nth_root(2, True)
885          (11, False)
886          sage: Integer(125).nth_root(3, True)
887          (5, True)
888
889          sage: Integer(125).nth_root(-5)
890          Traceback (most recent call last):
891          ...
892          ValueError: n (=-5) must be positive
893
894          sage: Integer(-25).nth_root(2)
895          Traceback (most recent call last):
896          ...
897          ValueError: cannot take even root of negative number
898
899        """
900        if n < 1:
901            raise ValueError, "n (=%s) must be positive" % n
902        if (mpz_sgn(self.value) < 0) and not (n & 1):
903            raise ValueError, "cannot take even root of negative number"
904        cdef Integer x
905        cdef bint is_exact
906        x = PY_NEW(Integer)
907        _sig_on
908        is_exact = mpz_root(x.value, self.value, n)
909        _sig_off
910
911        if report_exact:
912            return x, is_exact
913        else:
914            return x
915
916    def exact_log(self, m):
917        r"""
918        Returns the largest integer $k$ such that $m^k \leq \text{self}$, i.e.,
919        the floor of $\log_m(\text{self})$.
920
921        This is guaranteed to return the correct answer even when the usual
922        log function doesn't have sufficient precision.
923
924        INPUT:
925            m -- integer >= 2
926
927        AUTHOR:
928           -- David Harvey (2006-09-15)
929
930        TODO:
931           -- Currently this is extremely stupid code (although it should
932           always work). Someone needs to think about doing this properly
933           by estimating errors in the log function etc.
934
935        EXAMPLES:
936           sage: Integer(125).exact_log(5)
937           3
938           sage: Integer(124).exact_log(5)
939           2
940           sage: Integer(126).exact_log(5)
941           3
942           sage: Integer(3).exact_log(5)
943           0
944           sage: Integer(1).exact_log(5)
945           0
946
947
948           sage: x = 3^100000
949           sage: RR(log(RR(x), 3))
950           100000.000000000
951           sage: RR(log(RR(x + 100000), 3))
952           100000.000000000
953
954           sage: x.exact_log(3)
955           100000
956           sage: (x+1).exact_log(3)
957           100000
958           sage: (x-1).exact_log(3)
959           99999
960
961           sage: x.exact_log(2.5)
962           Traceback (most recent call last):
963           ...
964           ValueError: base of log must be an integer
965        """
966        _m = int(m)
967        if _m != m:
968            raise ValueError, "base of log must be an integer"
969        m = _m
970        if mpz_sgn(self.value) <= 0:
971            raise ValueError, "self must be positive"
972        if m < 2:
973            raise ValueError, "m must be at least 2"
974        import real_mpfr
975        R = real_mpfr.RealField(53)
976        guess = R(self).log(base = m).floor()
977        power = m ** guess
978
979        while power > self:
980            power = power / m
981            guess = guess - 1
982
983        if power == self:
984            return guess
985
986        while power < self:
987            power = power * m
988            guess = guess + 1
989
990        if power == self:
991            return guess
992        else:
993            return guess - 1
994
995
996    def __pos__(self):
997        """
998        EXAMPLES:
999            sage: z=43434
1000            sage: z.__pos__()
1001            43434
1002        """
1003        return self
1004
1005    def __abs__(self):
1006        """
1007        Computes $|self|$
1008
1009        EXAMPLES:
1010            sage: z = -1
1011            sage: abs(z)
1012            1
1013            sage: abs(z) == abs(1)
1014            True
1015        """
1016        cdef Integer x
1017        x = PY_NEW(Integer)
1018        mpz_abs(x.value, self.value)
1019        return x
1020
1021    def __mod__(self, modulus):
1022        r"""
1023        Returns self modulo the modulus.
1024
1025        EXAMPLES:
1026            sage: z = 43
1027            sage: z % 2
1028            1
1029            sage: z % 0
1030            Traceback (most recent call last):
1031            ...
1032            ZeroDivisionError: Integer modulo by zero
1033            sage: -5 % 7
1034            2
1035         """
1036        cdef Integer _modulus, _self
1037        _modulus = integer(modulus)
1038        if not _modulus:
1039            raise ZeroDivisionError, "Integer modulo by zero"
1040        _self = integer(self)
1041
1042        cdef Integer x
1043        x = PY_NEW(Integer)
1044
1045        _sig_on
1046        mpz_mod(x.value, _self.value, _modulus.value)
1047        _sig_off
1048
1049        return x
1050
1051
1052    def quo_rem(self, other):
1053        """
1054        Returns the quotient and the remainder of
1055        self divided by other.
1056
1057        INPUT:
1058            other -- the integer the divisor
1059
1060        OUTPUT:
1061            q   -- the quotient of self/other
1062            r   -- the remainder of self/other
1063
1064        EXAMPLES:
1065            sage: z = Integer(231)
1066            sage: z.quo_rem(2)
1067            (115, 1)
1068            sage: z.quo_rem(-2)
1069            (-115, 1)
1070            sage: z.quo_rem(0)
1071            Traceback (most recent call last):
1072            ...
1073            ZeroDivisionError: other (=0) must be nonzero
1074        """
1075        cdef Integer _other, _self
1076        _other = integer(other)
1077        if not _other:
1078            raise ZeroDivisionError, "other (=%s) must be nonzero"%other
1079        _self = integer(self)
1080
1081        cdef Integer q, r
1082        q = PY_NEW(Integer)
1083        r = PY_NEW(Integer)
1084
1085        _sig_on
1086        if mpz_sgn(_other.value) == 1:
1087            mpz_fdiv_qr(q.value, r.value, _self.value, _other.value)
1088        else:
1089            mpz_cdiv_qr(q.value, r.value, _self.value, _other.value)
1090        _sig_off
1091
1092        return q, r
1093
1094    def div(self, other):
1095        """
1096        Returns the quotient of self divided by other.
1097
1098        INPUT:
1099            other -- the integer the divisor
1100
1101        OUTPUT:
1102            q   -- the quotient of self/other
1103
1104        EXAMPLES:
1105            sage: z = Integer(231)
1106            sage: z.div(2)
1107            115
1108            sage: z.div(-2)
1109            -115
1110            sage: z.div(0)
1111            Traceback (most recent call last):
1112            ...
1113            ZeroDivisionError: other (=0) must be nonzero
1114        """
1115        cdef Integer _other, _self
1116        _other = integer(other)
1117        if not _other:
1118            raise ZeroDivisionError, "other (=%s) must be nonzero"%other
1119        _self = integer(self)
1120
1121        cdef Integer q, r
1122        q = PY_NEW(Integer)
1123        r = PY_NEW(Integer)
1124
1125        _sig_on
1126        mpz_tdiv_qr(q.value, r.value, _self.value, _other.value)
1127        _sig_off
1128
1129        return q
1130
1131
1132    def powermod(self, exp, mod):
1133        """
1134        Compute self**exp modulo mod.
1135
1136        EXAMPLES:
1137            sage: z = 2
1138            sage: z.powermod(31,31)
1139            2
1140            sage: z.powermod(0,31)
1141            1
1142            sage: z.powermod(-31,31) == 2^-31 % 31
1143            True
1144
1145            As expected, the following is invalid:
1146            sage: z.powermod(31,0)
1147            Traceback (most recent call last):
1148            ...
1149            ZeroDivisionError: cannot raise to a power modulo 0
1150        """
1151        cdef Integer x, _exp, _mod
1152        _exp = Integer(exp); _mod = Integer(mod)
1153        if mpz_cmp_si(_mod.value,0) == 0:
1154            raise ZeroDivisionError, "cannot raise to a power modulo 0"
1155
1156        x = PY_NEW(Integer)
1157
1158        _sig_on
1159        mpz_powm(x.value, self.value, _exp.value, _mod.value)
1160        _sig_off
1161
1162        return x
1163
1164    def rational_reconstruction(self, Integer m):
1165        import rational
1166        return rational.pyrex_rational_reconstruction(self, m)
1167
1168    def powermodm_ui(self, exp, mod):
1169        r"""
1170        Computes self**exp modulo mod, where exp is an unsigned
1171        long integer.
1172
1173        EXAMPLES:
1174            sage: z = 32
1175            sage: z.powermodm_ui(2, 4)
1176            0
1177            sage: z.powermodm_ui(2, 14)
1178            2
1179            sage: z.powermodm_ui(2^32-2, 14)
1180            2
1181            sage: z.powermodm_ui(2^32-1, 14)
1182            Traceback (most recent call last):                              # 32-bit
1183            ...                                                             # 32-bit
1184            OverflowError: exp (=4294967295) must be <= 4294967294          # 32-bit
1185            8              # 64-bit
1186            sage: z.powermodm_ui(2^65, 14)
1187            Traceback (most recent call last):
1188            ...
1189            OverflowError: exp (=36893488147419103232) must be <= 4294967294  # 32-bit
1190            OverflowError: exp (=36893488147419103232) must be <= 18446744073709551614     # 64-bit
1191        """
1192        if exp < 0:
1193            raise ValueError, "exp (=%s) must be nonnegative"%exp
1194        elif exp > MAX_UNSIGNED_LONG:
1195            raise OverflowError, "exp (=%s) must be <= %s"%(exp, MAX_UNSIGNED_LONG)
1196        cdef Integer x, _mod
1197        _mod = Integer(mod)
1198        x = PY_NEW(Integer)
1199
1200        _sig_on
1201        mpz_powm_ui(x.value, self.value, exp, _mod.value)
1202        _sig_off
1203
1204        return x
1205
1206    def __int__(self):
1207        return mpz_get_pyintlong(self.value)
1208
1209    def __long__(self):
1210        return mpz_get_pylong(self.value)
1211
1212    def __float__(self):
1213        return mpz_get_d(self.value)
1214
1215    def __hash__(self):
1216        return mpz_pythonhash(self.value)
1217
1218    def factor(self, algorithm='pari'):
1219        """
1220        Return the prime factorization of the integer as a list of
1221        pairs $(p,e)$, where $p$ is prime and $e$ is a positive integer.
1222
1223        INPUT:
1224            algorithm -- string
1225                 * 'pari' -- (default)  use the PARI c library
1226                 * 'kash' -- use KASH computer algebra system (requires
1227                             the optional kash package be installed)
1228        """
1229        import sage.rings.integer_ring
1230        return sage.rings.integer_ring.factor(self, algorithm=algorithm)
1231
1232    def coprime_integers(self, m):
1233        """
1234        Return the positive integers $< m$ that are coprime to self.
1235
1236        EXAMPLES:
1237            sage: n = 8
1238            sage: n.coprime_integers(8)
1239            [1, 3, 5, 7]
1240            sage: n.coprime_integers(11)
1241            [1, 3, 5, 7, 9]
1242            sage: n = 5; n.coprime_integers(10)
1243            [1, 2, 3, 4, 6, 7, 8, 9]
1244            sage: n.coprime_integers(5)
1245            [1, 2, 3, 4]
1246            sage: n = 99; n.coprime_integers(99)
1247            [1, 2, 4, 5, 7, 8, 10, 13, 14, 16, 17, 19, 20, 23, 25, 26, 28, 29, 31, 32, 34, 35, 37, 38, 40, 41, 43, 46, 47, 49, 50, 52, 53, 56, 58, 59, 61, 62, 64, 65, 67, 68, 70, 71, 73, 74, 76, 79, 80, 82, 83, 85, 86, 89, 91, 92, 94, 95, 97, 98]
1248
1249        AUTHORS:
1250            -- Naqi Jaffery (2006-01-24): examples
1251
1252        ALGORITHM: Naive -- compute lots of GCD's.  If this isn't good
1253        enough for you, please code something better and submit a
1254        patch.
1255        """
1256        # TODO -- make VASTLY faster
1257        v = []
1258        for n in range(1,m):
1259            if self.gcd(n) == 1:
1260                v.append(Integer(n))
1261        return v
1262
1263    def divides(self, n):
1264        """
1265        Return True if self divides n.
1266
1267        EXAMPLES:
1268            sage: Z = IntegerRing()
1269            sage: Z(5).divides(Z(10))
1270            True
1271            sage: Z(0).divides(Z(5))
1272            False
1273            sage: Z(10).divides(Z(5))
1274            False
1275        """
1276        cdef bint t
1277        cdef Integer _n
1278        _n = Integer(n)
1279        if mpz_sgn(self.value) == 0:
1280            return mpz_sgn(_n.value) == 0
1281        _sig_on
1282        t = mpz_divisible_p(_n.value, self.value)
1283        _sig_off
1284        return t
1285
1286    cdef RingElement _valuation(Integer self, Integer p):
1287        r"""
1288        Return the p-adic valuation of self.
1289
1290        We do not require that p be prime, but it must be at least 2.
1291        For more documentation see \code{valuation}
1292
1293        AUTHOR:
1294            -- David Roe (3/31/07)
1295        """
1296        if mpz_sgn(self.value) == 0:
1297            return sage.rings.infinity.infinity
1298        if mpz_cmp_ui(p.value, 2) < 0:
1299            raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
1300        cdef Integer v
1301        v = PY_NEW(Integer)
1302        cdef mpz_t u
1303        mpz_init(u)
1304        _sig_on
1305        mpz_set_ui(v.value, mpz_remove(u, self.value, p.value))
1306        _sig_off
1307        mpz_clear(u)
1308        return v
1309
1310    cdef object _val_unit(Integer self, Integer p):
1311        r"""
1312        Returns a pair: the p-adic valuation of self, and the p-adic unit of self.
1313
1314        We do not require the p be prime, but it must be at least 2.
1315        For more documentation see \code{val_unit}
1316
1317        AUTHOR:
1318            -- David Roe (3/31/07)
1319        """
1320        cdef Integer v, u
1321        if mpz_cmp_ui(p.value, 2) < 0:
1322            raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
1323        if self == 0:
1324            u = ONE
1325            Py_INCREF(ONE)
1326            return (sage.rings.infinity.infinity, u)
1327        v = PY_NEW(Integer)
1328        u = PY_NEW(Integer)
1329        _sig_on
1330        mpz_set_ui(v.value, mpz_remove(u.value, self.value, p.value))
1331        _sig_off
1332        return (v, u)
1333
1334    def valuation(self, p):
1335        """
1336        Return the p-adic valuation of self.
1337
1338        INPUT:
1339            p -- an integer at least 2.
1340
1341        EXAMPLE:
1342        sage: n = 60
1343        sage: n.valuation(2)
1344        2
1345        sage: n.valuation(3)
1346        1
1347        sage: n.valuation(7)
1348        0
1349        sage: n.valuation(1)
1350        Traceback (most recent call last):
1351        ...
1352        ValueError: You can only compute the valuation with respect to a integer larger than 1.
1353
1354        We do not require that p is a prime:
1355        sage: (2^11).valuation(4)
1356        5
1357        """
1358        return self._valuation(Integer(p))
1359
1360    def val_unit(self, p):
1361        r"""
1362        Returns a pair: the p-adic valuation of self, and the p-adic unit of self.
1363
1364        INPUT:
1365            p -- an integer at least 2.
1366
1367        OUTPUT:
1368            v_p(self) -- the p-adic valuation of self
1369            u_p(self) -- self / p^{v_p(self)}
1370
1371        EXAMPLE:
1372        sage: n = 60
1373        sage: n.val_unit(2)
1374        (2, 15)
1375        sage: n.val_unit(3)
1376        (1, 20)
1377        sage: n.val_unit(7)
1378        (0, 60)
1379        sage: (2^11).val_unit(4)
1380        (5, 2)
1381        sage: 0.val_unit(2)
1382        (+Infinity, 1)
1383        """
1384        return self._val_unit(Integer(p))
1385
1386    def ord(self, p=None):
1387        """Synonym for valuation
1388
1389        EXAMPLES:
1390        sage: n=12
1391        sage: n.ord(3)
1392        1
1393        """
1394        return self.valuation(p)
1395
1396    cdef Integer _divide_knowing_divisible_by(Integer self, Integer right):
1397        r"""
1398        Returns the integer self / right when self is divisible by right.
1399
1400        If self is not divisible by right, the return value is undefined, but seems to be close to self/right.
1401        For more documentation see \code{divide_knowing_divisible_by}
1402
1403        AUTHOR:
1404            -- David Roe (3/31/07)
1405        """
1406        if mpz_cmp_ui(right.value, 0) == 0:
1407            raise ZeroDivisionError
1408        cdef Integer x
1409        x = PY_NEW(Integer)
1410        if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
1411            # Only use the signal handler (to enable ctrl-c out) when the
1412            # quotient might take a while to compute
1413            _sig_on
1414            mpz_divexact(x.value, self.value, right.value)
1415            _sig_off
1416        else:
1417            mpz_divexact(x.value, self.value, right.value)
1418        return x
1419
1420    def divide_knowing_divisible_by(self, right):
1421        r"""
1422        Returns the integer self / right when self is divisible by right.
1423
1424        If self is not divisible by right, the return value is undefined, but seems to be close to self/right.
1425
1426        EXAMPLES:
1427        sage: a = 8; b = 4
1428        sage: a.divide_knowing_divisible_by(b)
1429        2
1430        """
1431        return self._divide_knowing_divisible_by(right)
1432
1433    def _lcm(self, Integer n):
1434        """
1435        Returns the least common multiple of self and $n$.
1436
1437        EXAMPLES:
1438            sage: n = 60
1439            sage: n._lcm(150)
1440            300
1441        """
1442        cdef mpz_t x
1443
1444        mpz_init(x)
1445
1446        _sig_on
1447        mpz_lcm(x, self.value, n.value)
1448        _sig_off
1449
1450
1451        cdef Integer z
1452        z = PY_NEW(Integer)
1453        mpz_set(z.value,x)
1454        mpz_clear(x)
1455        return z
1456
1457    def denominator(self):
1458        """
1459        Return the denominator of this integer.
1460
1461        EXAMPLES:
1462            sage: x = 5
1463            sage: x.denominator()
1464            1
1465            sage: x = 0
1466            sage: x.denominator()
1467            1
1468        """
1469        return ONE
1470
1471    def numerator(self):
1472        """
1473        Return the numerator of this integer.
1474
1475        EXAMPLE:
1476            sage: x = 5
1477            sage: x.numerator()
1478            5
1479
1480            sage: x = 0
1481            sage: x.numerator()
1482            0
1483        """
1484        return self
1485
1486    def factorial(self):
1487        """
1488        Return the factorial $n!=1 \\cdot 2 \\cdot 3 \\cdots n$.
1489        Self must fit in an \\code{unsigned long int}.
1490
1491        EXAMPLES:
1492            sage: for n in srange(7):
1493            ...    print n, n.factorial()
1494            0 1
1495            1 1
1496            2 2
1497            3 6
1498            4 24
1499            5 120
1500            6 720
1501        """
1502        if mpz_sgn(self.value) < 0:
1503            raise ValueError, "factorial -- self = (%s) must be nonnegative"%self
1504
1505        if mpz_cmp_ui(self.value,4294967295) > 0:
1506            raise ValueError, "factorial not implemented for n >= 2^32.\nThis is probably OK, since the answer would have billions of digits."
1507
1508        cdef unsigned int n
1509        n = self
1510
1511        cdef mpz_t x
1512        cdef Integer z
1513
1514        mpz_init(x)
1515
1516        _sig_on
1517        mpz_fac_ui(x, n)
1518        _sig_off
1519
1520        z = PY_NEW(Integer)
1521        set_mpz(z, x)
1522        mpz_clear(x)
1523        return z
1524
1525    def floor(self):
1526        """
1527        Return the floor of self, which is just self since self is an integer.
1528
1529        EXAMPLES:
1530            sage: n = 6
1531            sage: n.floor()
1532            6
1533        """
1534        return self
1535
1536    def ceil(self):
1537        """
1538        Return the ceiling of self, which is self since self is an integer.
1539
1540        EXAMPLES:
1541            sage: n = 6
1542            sage: n.ceil()
1543            6
1544        """
1545        return self
1546
1547    def is_one(self):
1548        r"""
1549        Returns \code{True} if the integers is $1$, otherwise \code{False}.
1550
1551        EXAMPLES:
1552            sage: Integer(1).is_one()
1553            True
1554            sage: Integer(0).is_one()
1555            False
1556        """
1557        return mpz_cmp_si(self.value, 1) == 0
1558
1559    def __nonzero__(self):
1560        r"""
1561        Returns \code{True} if the integers is not $0$, otherwise \code{False}.
1562
1563        EXAMPLES:
1564            sage: Integer(1).is_zero()
1565            False
1566            sage: Integer(0).is_zero()
1567            True
1568        """
1569        return mpz_sgn(self.value) != 0
1570
1571    def is_unit(self):
1572        r"""
1573        Returns \code{true} if this integer is a unit, i.e., 1 or $-1$.
1574
1575        EXAMPLES:
1576        sage: for n in srange(-2,3):
1577        ...    print n, n.is_unit()
1578        -2 False
1579        -1 True
1580        0 False
1581        1 True
1582        2 False
1583        """
1584        return mpz_cmp_si(self.value, -1) == 0 or mpz_cmp_si(self.value, 1) == 0
1585
1586    def is_square(self):
1587        r"""
1588        Returns \code{True} if self is a perfect square
1589
1590        EXAMPLES:
1591        sage: Integer(4).is_square()
1592        True
1593        sage: Integer(41).is_square()
1594        False
1595        """
1596        return mpz_perfect_square_p(self.value)
1597
1598    def is_power(self):
1599        r"""
1600        Returns \code{True} if self is a perfect power, ie if there
1601        exist integers a and b, $b > 1$ with $self = a^b$.
1602
1603        EXAMPLES:
1604        sage: Integer(-27).is_power()
1605        True
1606        sage: Integer(12).is_power()
1607        False
1608        """
1609        return mpz_perfect_power_p(self.value)
1610
1611    cdef bint _is_power_of(Integer self, Integer n):
1612        r"""
1613        Returns a non-zero int if there is an integer b with $self = n^b$.
1614
1615        For more documentation see \code{is_power_of}
1616
1617        AUTHOR:
1618            -- David Roe (3/31/07)
1619        """
1620        cdef int a
1621        cdef unsigned long b, c
1622        cdef mpz_t u, sabs, nabs
1623        a = mpz_cmp_ui(n.value, 2)
1624        if a <= 0: # n <= 2
1625            if a == 0: # n == 2
1626                if mpz_popcount(self.value) == 1: #number of bits set in self == 1
1627                    return 1
1628                else:
1629                    return 0
1630            a = mpz_cmp_si(n.value, -2)
1631            if a >= 0: # -2 <= n < 2:
1632                a = mpz_get_si(n.value)
1633                if a == 1: # n == 1
1634                    if mpz_cmp_ui(self.value, 1) == 0: # Only 1 is a power of 1
1635                        return 1
1636                    else:
1637                        return 0
1638                elif a == 0: # n == 0
1639                    if mpz_cmp_ui(self.value, 0) == 0 or mpz_cmp_ui(self.value, 1) == 0: # 0^0 = 1, 0^x = 0
1640                        return 1
1641                    else:
1642                        return 0
1643                elif a == -1: # n == -1
1644                    if mpz_cmp_ui(self.value, 1) == 0 or mpz_cmp_si(self.value, -1) == 0: # 1 and -1 are powers of -1
1645                        return 1
1646                    else:
1647                        return 0
1648                elif a == -2: # n == -2
1649                    mpz_init(sabs)
1650                    mpz_abs(sabs, self.value)
1651                    if mpz_popcount(sabs) == 1: # number of bits set in |self| == 1
1652                        b = mpz_scan1(sabs, 0) % 2 # b == 1 if |self| is an odd power of 2, 0 if |self| is an even power
1653                        mpz_clear(sabs)
1654                        if (b == 1 and mpz_cmp_ui(self.value, 0) < 0) or (b == 0 and mpz_cmp_ui(self.value, 0) > 0):
1655                            # An odd power of -2 is negative, an even power must be positive.
1656                            return 1
1657                        else: # number of bits set in |self| is not 1, so self cannot be a power of -2
1658                            return 0
1659                    else: # |self| is not a power of 2, so self cannot be a power of -2
1660                        return 0
1661            else: # n < -2
1662                mpz_init(nabs)
1663                mpz_neg(nabs, n.value)
1664                if mpz_popcount(nabs) == 1: # |n| = 2^k for k >= 2.  We special case this for speed
1665                    mpz_init(sabs)
1666                    mpz_abs(sabs, self.value)
1667                    if mpz_popcount(sabs) == 1: # |self| = 2^L for some L >= 0.
1668                        b = mpz_scan1(sabs, 0) # the bit that self is set at
1669                        c = mpz_scan1(nabs, 0) # the bit that n is set at
1670                        # Having obtained b and c, we're done with nabs and sabs (on this branch anyway)
1671                        mpz_clear(nabs)
1672                        mpz_clear(sabs)
1673                        if b % c == 0: # Now we know that |self| is a power of |n|
1674                            b = (b // c) % 2 # Whether b // c is even or odd determines whether (-2^c)^(b // c) is positive or negative
1675                            a = mpz_cmp_ui(self.value, 0)
1676                            if b == 0 and a > 0 or b == 1 and a < 0:
1677                                # These two cases are that b // c is even and self positive, or b // c is odd and self negative
1678                                return 1
1679                            else: # The sign of self is wrong
1680                                return 0
1681                        else: # Since |self| is not a power of |n|, self cannot be a power of n
1682                            return 0
1683                    else: # self is not a power of 2, and thus cannot be a power of n, which is a power of 2.
1684                        mpz_clear(nabs)
1685                        mpz_clear(sabs)
1686                        return 0
1687                else: # |n| is not a power of 2, so we use mpz_remove
1688                    mpz_init(u)
1689                    _sig_on
1690                    b = mpz_remove(u, self.value, nabs)
1691                    _sig_off
1692                    # Having obtained b and u, we're done with nabs
1693                    mpz_clear(nabs)
1694                    if mpz_cmp_ui(u, 1) == 0: # self is a power of |n|
1695                        mpz_clear(u)
1696                        if b % 2 == 0: # an even power of |n|, and since self > 0, this means that self is a power of n
1697                            return 1
1698                        else:
1699                            return 0
1700                    elif mpz_cmp_si(u, -1) == 0: # -self is a power of |n|
1701                        mpz_clear(u)
1702                        if b % 2 == 1: # an odd power of |n|, and thus self is a power of n
1703                            return 1
1704                        else:
1705                            return 1
1706                    else: # |self| is not a power of |n|, so self cannot be a power of n
1707                        mpz_clear(u)
1708                        return 0
1709        elif mpz_popcount(n.value) == 1: # n > 2 and in fact n = 2^k for k >= 2
1710            if mpz_popcount(self.value) == 1: # since n is a power of 2, so must self be.
1711                if mpz_scan1(self.value, 0) % mpz_scan1(n.value, 0) == 0: # log_2(self) is divisible by log_2(n)
1712                    return 1
1713                else:
1714                    return 0
1715            else: # self is not a power of 2, and thus not a power of n
1716                return 0
1717        else: # n > 2, but not a power of 2, so we use mpz_remove
1718            mpz_init(u)
1719            _sig_on
1720            mpz_remove(u, self.value, n.value)
1721            _sig_off
1722            a = mpz_cmp_ui(u, 1)
1723            mpz_clear(u)
1724            if a == 0:
1725                return 1
1726            else:
1727                return 0
1728
1729    def is_power_of(Integer self, n):
1730        r"""
1731        Returns \code{True} if there is an integer b with $self = n^b$.
1732
1733        EXAMPLES:
1734        sage: Integer(64).is_power_of(4)
1735        True
1736        sage: Integer(64).is_power_of(16)
1737        False
1738
1739        TESTS:
1740        sage: Integer(-64).is_power_of(-4)
1741        True
1742        sage: Integer(-32).is_power_of(-2)
1743        True
1744        sage: Integer(1).is_power_of(1)
1745        True
1746        sage: Integer(-1).is_power_of(-1)
1747        True
1748        sage: Integer(0).is_power_of(1)
1749        False
1750        sage: Integer(0).is_power_of(0)
1751        True
1752        sage: Integer(1).is_power_of(0)
1753        True
1754        sage: Integer(1).is_power_of(8)
1755        True
1756        sage: Integer(-8).is_power_of(2)
1757        False
1758
1759        NOTES:
1760        For large integers self, is_power_of() is faster than is_power().  The following examples gives some indication of how much faster.
1761        sage.: b = lcm(range(1,10000))
1762        sage.: b.exact_log(2)
1763        14446
1764        sage.: time for a in range(2, 1000): k = b.is_power()
1765        CPU times: user 1.68 s, sys: 0.01 s, total: 1.69 s
1766        Wall time: 1.71
1767        sage.: time for a in range(2, 1000): k = b.is_power_of(2)
1768        CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
1769        Wall time: 0.01
1770        sage.: time for a in range(2, 1000): k = b.is_power_of(3)
1771        CPU times: user 0.05 s, sys: 0.00 s, total: 0.05 s
1772        Wall time: 0.05
1773
1774        sage.: b = lcm(range(1, 1000))
1775        sage.: b.exact_log(2)
1776        1437
1777        sage.: time for a in range(2, 10000): k = b.is_power() # note that we change the range from the example above
1778        CPU times: user 0.40 s, sys: 0.00 s, total: 0.40 s
1779        Wall time: 0.40
1780        sage.: for a in range(2, 10000): k = b.is_power_of(TWO)
1781        CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
1782        Wall time: 0.03
1783        sage.: time for a in range(2, 10000): k = b.is_power_of(3)
1784        CPU times: user 0.11 s, sys: 0.01 s, total: 0.12 s
1785        Wall time: 0.13
1786        sage.: time for a in range(2, 10000): k = b.is_power_of(a)
1787        CPU times: user 0.08 s, sys: 0.01 s, total: 0.09 s
1788        Wall time: 0.10
1789        """
1790        if not PY_TYPE_CHECK(n, Integer):
1791            n = Integer(n)
1792        return self._is_power_of(n)
1793
1794    def is_prime_power(self, flag=0):
1795        r"""
1796        Returns True if $x$ is a prime power, and False otherwise.
1797
1798        INPUT:
1799            flag (for primality testing) -- int
1800                    0 (default): use a combination of algorithms.
1801                    1: certify primality using the Pocklington-Lehmer Test.
1802                    2: certify primality using the APRCL test.
1803
1804        EXAMPLES:
1805            sage: (-10).is_prime_power()
1806            False
1807            sage: (10).is_prime_power()
1808            False
1809            sage: (64).is_prime_power()
1810            True
1811            sage: (3^10000).is_prime_power()
1812            True
1813            sage: (10000).is_prime_power(flag=1)
1814            False
1815        """
1816        if self.is_zero():
1817            return False
1818        elif self.is_one():
1819            return True
1820        elif mpz_sgn(self.value) < 0:
1821            return False
1822        if self.is_prime():
1823            return True
1824        if not self.is_perfect_power():
1825            return False
1826        k, g = self._pari_().ispower()
1827        if not k:
1828            raise RuntimeError, "inconsistent results between GMP and pari"
1829        return g.isprime(flag=flag)
1830
1831    def is_prime(self):
1832        r"""
1833        Retuns \code{True} if self is prime
1834
1835        EXAMPLES:
1836            sage: z = 2^31 - 1
1837            sage: z.is_prime()
1838            True
1839            sage: z = 2^31
1840            sage: z.is_prime()
1841            False
1842        """
1843        return bool(self._pari_().isprime())
1844
1845    def is_pseudoprime(self):
1846        r"""
1847        Retuns \code{True} if self is a pseudoprime
1848
1849        EXAMPLES:
1850            sage: z = 2^31 - 1
1851            sage: z.is_pseudoprime()
1852            True
1853            sage: z = 2^31
1854            sage: z.is_pseudoprime()
1855            False
1856        """
1857        return bool(self._pari_().ispseudoprime())
1858
1859    def is_perfect_power(self):
1860        r"""
1861        Retuns \code{True} if self is a perfect power.
1862
1863        EXAMPLES:
1864            sage: z = 8
1865            sage: z.is_perfect_power()
1866            True
1867            sage: z = 144
1868            sage: z.is_perfect_power()
1869            True
1870            sage: z = 10
1871            sage: z.is_perfect_power()
1872            False
1873        """
1874        return mpz_perfect_power_p(self.value)
1875
1876    def jacobi(self, b):
1877        r"""
1878        Calculate the Jacobi symbol $\left(\frac{self}{b}\right)$.
1879
1880        EXAMPLES:
1881            sage: z = -1
1882            sage: z.jacobi(17)
1883            1
1884            sage: z.jacobi(19)
1885            -1
1886            sage: z.jacobi(17*19)
1887            -1
1888            sage: (2).jacobi(17)
1889            1
1890            sage: (3).jacobi(19)
1891            -1
1892            sage: (6).jacobi(17*19)
1893            -1
1894            sage: (6).jacobi(33)
1895            0
1896            sage: a = 3; b = 7
1897            sage: a.jacobi(b) == -b.jacobi(a)
1898            True
1899        """
1900        cdef long tmp
1901        if PY_TYPE_CHECK(b, int):
1902            tmp = b
1903            if (tmp & 1) == 0:
1904                raise ValueError, "Jacobi symbol not defined for even b."
1905            return mpz_kronecker_si(self.value, tmp)
1906        if not PY_TYPE_CHECK(b, Integer):
1907            b = Integer(b)
1908        if mpz_even_p((<Integer>b).value):
1909            raise ValueError, "Jacobi symbol not defined for even b."
1910        return mpz_jacobi(self.value, (<Integer>b).value)
1911
1912    def kronecker(self, b):
1913        r"""
1914        Calculate the Kronecker symbol
1915        $\left(\frac{self}{b}\right)$ with the Kronecker extension
1916        $(self/2)=(2/self)$ when self odd, or $(self/2)=0$ when $self$ even.
1917
1918        EXAMPLES:
1919        EXAMPLES:
1920            sage: z = 5
1921            sage: z.kronecker(41)
1922            1
1923            sage: z.kronecker(43)
1924            -1
1925            sage: z.kronecker(8)
1926            -1
1927            sage: z.kronecker(15)
1928            0
1929            sage: a = 2; b = 5
1930            sage: a.kronecker(b) == b.kronecker(a)
1931            True
1932        """
1933        if PY_TYPE_CHECK(b, int):
1934            return mpz_kronecker_si(self.value, b)
1935        if not PY_TYPE_CHECK(b, Integer):
1936            b = Integer(b)
1937        return mpz_kronecker(self.value, (<Integer>b).value)
1938
1939    def square_free_part(self):
1940        """
1941        Return the square free part of $x$, i.e., a divisor z such that $x = z y^2$,
1942        for a perfect square $y^2$.
1943
1944        EXAMPLES:
1945            sage: square_free_part(100)
1946            1
1947            sage: square_free_part(12)
1948            3
1949            sage: square_free_part(17*37*37)
1950            17
1951            sage: square_free_part(-17*32)
1952            -34
1953            sage: square_free_part(1)
1954            1
1955            sage: square_free_part(-1)
1956            -1
1957            sage: square_free_part(-2)
1958            -2
1959            sage: square_free_part(-4)
1960            -1
1961        """
1962        if self.is_zero():
1963            return self
1964        F = self.factor()
1965        n = Integer(1)
1966        for p, e in F:
1967            if e % 2 != 0:
1968                n = n * p
1969        return n * F.unit()
1970
1971    def next_probable_prime(self):
1972        """
1973        Returns the next probable prime after self, as determined by PARI.
1974
1975        EXAMPLES:
1976            sage: (-37).next_probable_prime()
1977            2
1978            sage: (100).next_probable_prime()
1979            101
1980            sage: (2^512).next_probable_prime()
1981            13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084171
1982            sage: 0.next_probable_prime()
1983            2
1984            sage: 126.next_probable_prime()
1985            127
1986            sage: 144168.next_probable_prime()
1987            144169
1988        """
1989        return Integer( (self._pari_()+1).nextprime())
1990
1991    def next_prime(self, proof=True):
1992        r"""
1993        Returns the next prime after self.
1994
1995        INPUT:
1996            proof -- bool (default: True)
1997
1998        EXAMPLES:
1999            sage: Integer(100).next_prime()
2000            101
2001
2002        Use Proof = False, which is way faster:
2003            sage: b = (2^1024).next_prime(proof=False)
2004
2005            sage: Integer(0).next_prime()
2006            2
2007            sage: Integer(1001).next_prime()
2008            1009
2009        """
2010        if self < 2:   # negatives are not prime.
2011            return integer_ring.ZZ(2)
2012        if self == 2:
2013            return integer_ring.ZZ(3)
2014        if not proof:
2015            return Integer( (self._pari_()+1).nextprime())
2016        n = self
2017        if n % 2 == 0:
2018            n += 1
2019        else:
2020            n += 2
2021        while not n.is_prime():  # pari isprime is provably correct
2022            n += 2
2023        return integer_ring.ZZ(n)
2024
2025
2027        """
2028        Return the additive order of self.
2029
2030        EXAMPLES:
2032            1
2034            +Infinity
2035        """
2036        import sage.rings.infinity
2037        if self.is_zero():
2038            return Integer(1)
2039        else:
2040            return sage.rings.infinity.infinity
2041
2042    def multiplicative_order(self):
2043        r"""
2044        Return the multiplicative order of self.
2045
2046        EXAMPLES:
2047            sage: ZZ(1).multiplicative_order()
2048            1
2049            sage: ZZ(-1).multiplicative_order()
2050            2
2051            sage: ZZ(0).multiplicative_order()
2052            +Infinity
2053            sage: ZZ(2).multiplicative_order()
2054            +Infinity
2055        """
2056        import sage.rings.infinity
2057        if  mpz_cmp_si(self.value, 1) == 0:
2058                return Integer(1)
2059        elif mpz_cmp_si(self.value, -1) == 0:
2060                return Integer(2)
2061        else:
2062                return sage.rings.infinity.infinity
2063
2064    def is_squarefree(self):
2065        """
2066        Returns True if this integer is not divisible by the square of
2067        any prime and False otherwise.
2068
2069        EXAMPLES:
2070            sage: Integer(100).is_squarefree()
2071            False
2072            sage: Integer(102).is_squarefree()
2073            True
2074        """
2075        return self._pari_().issquarefree()
2076
2077    def _pari_(self):
2078        """
2079        Returns the PARI version of this integer.
2080
2081        EXAMPLES:
2082            sage: n = 9390823
2083            sage: m = n._pari_(); m
2084            9390823
2085            sage: type(m)
2086            <type 'sage.libs.pari.gen.gen'>
2087
2088        ALGORITHM: Use base 10 Python string conversion, hence very
2089        very slow for large integers. If you can figure out how to
2090        input a number into PARI in hex, or otherwise optimize this,
2091        please implement it and send me a patch.
2092        """
2093        #if self._pari is None:
2094            # better to do in hex, but I can't figure out
2095            # how to input/output a number in hex in PARI!!
2096            # TODO: (I could just think carefully about raw bytes and make this all much faster...)
2097            #self._pari = sage.libs.pari.all.pari(str(self))
2098        #return self._pari
2099        return sage.libs.pari.all.pari(str(self))
2100
2101    def _interface_init_(self):
2102        """
2103        Return canonical string to coerce this integer to any other math
2104        software, i.e., just the string representation of this integer
2105        in base 10.
2106
2107        EXAMPLES:
2108            sage: n = 9390823
2109            sage: n._interface_init_()
2110            '9390823'
2111        """
2112        return str(self)
2113
2114    def isqrt(self):
2115        r"""
2116        Returns the integer floor of the square root of self, or raises
2117        an \exception{ValueError} if self is negative.
2118
2119        EXAMPLE:
2120            sage: a = Integer(5)
2121            sage: a.isqrt()
2122            2
2123
2124            sage: Integer(-102).isqrt()
2125            Traceback (most recent call last):
2126            ...
2127            ValueError: square root of negative integer not defined.
2128        """
2129        if mpz_sgn(self.value) < 0:
2130            raise ValueError, "square root of negative integer not defined."
2131        cdef Integer x
2132        x = PY_NEW(Integer)
2133
2134        _sig_on
2135        mpz_sqrt(x.value, self.value)
2136        _sig_off
2137
2138        return x
2139
2140    def sqrt_approx(self, prec=None, all=False):
2141        """
2142        EXAMPLES:
2143            sage: 5.sqrt_approx(prec=200)
2144            2.2360679774997896964091736687312762354406183596115257242709
2145            sage: 5.sqrt_approx()
2146            2.23606797749979
2147            sage: 4.sqrt_approx()
2148            2
2149        """
2150        try:
2151            return self.sqrt(extend=False,all=all)
2152        except ValueError:
2153            pass
2154        if prec is None:
2155            prec = max(53, 2*(mpz_sizeinbase(self.value, 2)+2))
2156        return self.sqrt(prec=prec, all=all)
2157
2158    def sqrt(self, prec=None, extend=True, all=False):
2159        """
2160        The square root function.
2161
2162        INPUT:
2163            prec -- integer (default: None): if None, returns an exact
2164                 square root; otherwise returns a numerical square
2165                 root if necessary, to the given bits of precision.
2166            extend -- bool (default: True); if True, return a square
2167                 root in an extension ring, if necessary. Otherwise,
2168                 raise a ValueError if the square is not in the base
2169                 ring.
2170            all -- bool (default: False); if True, return all square
2171                 roots of self, instead of just one.
2172
2173        EXAMPLES:
2174            sage: Integer(144).sqrt()
2175            12
2176            sage: Integer(102).sqrt()
2177            sqrt(102)
2178
2179            sage: n = 2
2180            sage: n.sqrt(all=True)
2181            [sqrt(2), -sqrt(2)]
2182            sage: n.sqrt(prec=10)
2183            1.4
2184            sage: n.sqrt(prec=100)
2185            1.4142135623730950488016887242
2186            sage: n.sqrt(prec=100,all=True)
2187            [1.4142135623730950488016887242, -1.4142135623730950488016887242]
2188            sage: n.sqrt(extend=False)
2189            Traceback (most recent call last):
2190            ...
2191            ValueError: square root of 2 not an integer
2192            sage: Integer(144).sqrt(all=True)
2193            [12, -12]
2194            sage: Integer(0).sqrt(all=True)
2195            [0]
2196        """
2197        if mpz_sgn(self.value) == 0:
2198            return [self] if all else self
2199
2200        if mpz_sgn(self.value) < 0:
2201            if not extend:
2202                raise ValueError, "square root of negative number not an integer"
2203            if prec:
2204                from sage.rings.complex_field import ComplexField
2205                K = ComplexField(prec)
2206                return K(self).sqrt(all=all)
2207            from sage.calculus.calculus import sqrt
2208            return sqrt(self, all=all)
2209
2210
2211        cdef int non_square
2212        cdef Integer z = PY_NEW(Integer)
2213        cdef mpz_t tmp
2214        _sig_on
2215        mpz_init(tmp)
2216        mpz_sqrtrem(z.value, tmp, self.value)
2217        non_square = mpz_sgn(tmp) != 0
2218        mpz_clear(tmp)
2219        _sig_off
2220
2221        if non_square:
2222            if not extend:
2223                raise ValueError, "square root of %s not an integer"%self
2224            if prec:
2225                from sage.rings.real_mpfr import RealField
2226                K = RealField(prec)
2227                return K(self).sqrt(all=all)
2228            from sage.calculus.calculus import sqrt
2229            return sqrt(self, all=all)
2230
2231        if all:
2232           return [z, -z]
2233        return z
2234
2235    def _xgcd(self, Integer n):
2236        r"""
2237        Return a triple $g, s, t \in\Z$ such that
2238        $$2239 g = s \cdot \mbox{\rm self} + t \cdot n. 2240$$
2241
2242        EXAMPLES:
2243            sage: n = 6
2244            sage: g, s, t = n._xgcd(8)
2245            sage: s*6 + 8*t
2246            2
2247            sage: g
2248            2
2249        """
2250        cdef mpz_t g, s, t
2251        cdef object g0, s0, t0
2252
2253        mpz_init(g)
2254        mpz_init(s)
2255        mpz_init(t)
2256
2257        _sig_on
2258        mpz_gcdext(g, s, t, self.value, n.value)
2259        _sig_off
2260
2261        g0 = PY_NEW(Integer)
2262        s0 = PY_NEW(Integer)
2263        t0 = PY_NEW(Integer)
2264        set_mpz(g0,g)
2265        set_mpz(s0,s)
2266        set_mpz(t0,t)
2267        mpz_clear(g)
2268        mpz_clear(s)
2269        mpz_clear(t)
2270        return g0, s0, t0
2271
2272    cdef _lshift(self, long int n):
2273        """
2274        Shift self n bits to the left, i.e., quickly multiply by $2^n$.
2275        """
2276        cdef Integer x
2277        x = <Integer> PY_NEW(Integer)
2278
2279        _sig_on
2280        if n < 0:
2281            mpz_fdiv_q_2exp(x.value, self.value, -n)
2282        else:
2283            mpz_mul_2exp(x.value, self.value, n)
2284        _sig_off
2285        return x
2286
2287    def __lshift__(x,y):
2288        """
2289        Shift x y bits to the left.
2290
2291        EXAMPLES:
2292            sage: 32 << 2
2293            128
2294            sage: 32 << int(2)
2295            128
2296            sage: int(32) << 2
2297            128
2298            sage: 1 >> 2.5
2299            Traceback (most recent call last):
2300            ...
2301            TypeError: unsupported operands for >>
2302        """
2303        try:
2304            if not PY_TYPE_CHECK(x, Integer):
2305                x = Integer(x)
2306            elif not PY_TYPE_CHECK(y, Integer):
2307                y = Integer(y)
2308            return (<Integer>x)._lshift(long(y))
2309        except TypeError:
2310            raise TypeError, "unsupported operands for <<"
2311        if PY_TYPE_CHECK(x, Integer) and isinstance(y, (Integer, int, long)):
2312            return (<Integer>x)._lshift(long(y))
2313        return bin_op(x, y, operator.lshift)
2314
2315    cdef _rshift(Integer self, long int n):
2316        cdef Integer x
2317        x = <Integer> PY_NEW(Integer)
2318
2319        _sig_on
2320        if n < 0:
2321            mpz_mul_2exp(x.value, self.value, -n)
2322        else:
2323            mpz_fdiv_q_2exp(x.value, self.value, n)
2324        _sig_off
2325        return x
2326
2327    def __rshift__(x, y):
2328        """
2329        EXAMPLES:
2330            sage: 32 >> 2
2331            8
2332            sage: 32 >> int(2)
2333            8
2334            sage: int(32) >> 2
2335            8
2336            sage: 1<< 2.5
2337            Traceback (most recent call last):
2338            ...
2339            TypeError: unsupported operands for <<
2340        """
2341        try:
2342            if not PY_TYPE_CHECK(x, Integer):
2343                x = Integer(x)
2344            elif not PY_TYPE_CHECK(y, Integer):
2345                y = Integer(y)
2346            return (<Integer>x)._rshift(long(y))
2347        except TypeError:
2348            raise TypeError, "unsupported operands for >>"
2349
2350        #if PY_TYPE_CHECK(x, Integer) and isinstance(y, (Integer, int, long)):
2351        #    return (<Integer>x)._rshift(long(y))
2352        #return bin_op(x, y, operator.rshift)
2353
2354    cdef _and(Integer self, Integer other):
2355        cdef Integer x
2356        x = PY_NEW(Integer)
2357        mpz_and(x.value, self.value, other.value)
2358        return x
2359
2360    def __and__(x, y):
2361        if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
2362            return (<Integer>x)._and(y)
2363        return bin_op(x, y, operator.and_)
2364
2365
2366    cdef _or(Integer self, Integer other):
2367        cdef Integer x
2368        x = PY_NEW(Integer)
2369        mpz_ior(x.value, self.value, other.value)
2370        return x
2371
2372    def __or__(x, y):
2373        """
2374        Return the bitwise or of the integers x and y.
2375
2376        EXAMPLES:
2377            sage: n = 8; m = 4
2378            sage: n.__or__(m)
2379            12
2380        """
2381        if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
2382            return (<Integer>x)._or(y)
2383        return bin_op(x, y, operator.or_)
2384
2385
2386    def __invert__(self):
2387        """
2388        Return the multiplicative interse of self, as a rational number.
2389
2390        EXAMPLE:
2391            sage: n = 10
2392            sage: 1/n
2393            1/10
2394            sage: n.__invert__()
2395            1/10
2396        """
2397        return Integer(1)/self    # todo: optimize
2398
2399    def inverse_of_unit(self):
2400        if mpz_cmp_si(self.value, 1) or mpz_cmp_si(self.value, -1):
2401            return self
2402        else:
2403            raise ZeroDivisionError, "Inverse does not exist."
2404
2405    def inverse_mod(self, n):
2406        """
2407        Returns the inverse of self modulo $n$, if this inverse exists.
2408        Otherwise, raises a \exception{ZeroDivisionError} exception.
2409
2410        INPUT:
2411           self -- Integer
2412           n -- Integer
2413        OUTPUT:
2414           x -- Integer such that x*self = 1 (mod m), or
2415                raises ZeroDivisionError.
2416        IMPLEMENTATION:
2417           Call the mpz_invert GMP library function.
2418
2419        EXAMPLES:
2420            sage: a = Integer(189)
2421            sage: a.inverse_mod(10000)
2422            4709
2423            sage: a.inverse_mod(-10000)
2424            4709
2425            sage: a.inverse_mod(1890)
2426            Traceback (most recent call last):
2427            ...
2428            ZeroDivisionError: Inverse does not exist.
2429            sage: a = Integer(19)**100000
2430            sage: b = a*a
2431            sage: c = a.inverse_mod(b)
2432            Traceback (most recent call last):
2433            ...
2434            ZeroDivisionError: Inverse does not exist.
2435        """
2436        cdef mpz_t x
2437        cdef object ans
2438        cdef int r
2439        cdef Integer m
2440        m = Integer(n)
2441
2442        if m == 1:
2443            return Integer(0)
2444
2445        mpz_init(x)
2446
2447        _sig_on
2448        r = mpz_invert(x, self.value, m.value)
2449        _sig_off
2450
2451        if r == 0:
2452            raise ZeroDivisionError, "Inverse does not exist."
2453        ans = PY_NEW(Integer)
2454        set_mpz(ans,x)
2455        mpz_clear(x)
2456        return ans
2457
2458    def gcd(self, n):
2459        """
2460        Return the greatest common divisor of self and $n$.
2461
2462        EXAMPLE:
2463            sage: gcd(-1,1)
2464            1
2465            sage: gcd(0,1)
2466            1
2467            sage: gcd(0,0)
2468            0
2469            sage: gcd(2,2^6)
2470            2
2471            sage: gcd(21,2^6)
2472            1
2473        """
2474        cdef mpz_t g
2475        cdef object g0
2476        cdef Integer _n = Integer(n)
2477
2478        mpz_init(g)
2479
2480
2481        _sig_on
2482        mpz_gcd(g, self.value, _n.value)
2483        _sig_off
2484
2485        g0 = PY_NEW(Integer)
2486        set_mpz(g0,g)
2487        mpz_clear(g)
2488        return g0
2489
2490    def crt(self, y, m, n):
2491        """
2492        Return the unique integer between $0$ and $mn$ that is
2493        congruent to the integer modulo $m$ and to $y$ modulo $n$.  We
2494        assume that~$m$ and~$n$ are coprime.
2495        """
2496        cdef object g, s, t
2497        cdef Integer _y, _m, _n
2498        _y = Integer(y); _m = Integer(m); _n = Integer(n)
2499        g, s, t = _m.xgcd(_n)
2500        if not g.is_one():
2501            raise ArithmeticError, "CRT requires that gcd of moduli is 1."
2502        # Now s*m + t*n = 1, so the answer is x + (y-x)*s*m, where x=self.
2503        return (self + (_y-self)*s*_m) % (_m*_n)
2504
2505    def test_bit(self, index):
2506        r"""
2507        Return the bit at \code{index}.
2508
2509        EXAMPLES:
2510            sage: w = 6
2511            sage: w.str(2)
2512            '110'
2513            sage: w.test_bit(2)
2514            1
2515            sage: w.test_bit(-1)
2516            0
2517        """
2518        cdef unsigned long int i
2519        i = index
2520        cdef Integer x
2521        x = Integer(self)
2522        return mpz_tstbit(x.value, i)
2523
2524
2525ONE = Integer(1)
2526
2527
2528def LCM_list(v):
2529    cdef int i, n
2530    cdef mpz_t z
2531    cdef Integer w
2532
2533    n = len(v)
2534
2535    if n == 0:
2536        return Integer(1)
2537
2538    try:
2539        w = v[0]
2540        mpz_init_set(z, w.value)
2541
2542        _sig_on
2543        for i from 1 <= i < n:
2544            w = v[i]
2545            mpz_lcm(z, z, w.value)
2546        _sig_off
2547    except TypeError:
2548        w = Integer(v[0])
2549        mpz_init_set(z, w.value)
2550
2551        _sig_on
2552        for i from 1 <= i < n:
2553            w = Integer(v[i])
2554            mpz_lcm(z, z, w.value)
2555        _sig_off
2556
2557
2558    w = PY_NEW(Integer)
2559    mpz_set(w.value, z)
2560    mpz_clear(z)
2561    return w
2562
2563
2564
2565def GCD_list(v):
2566    cdef int i, n
2567    cdef mpz_t z
2568    cdef Integer w
2569
2570    n = len(v)
2571
2572    if n == 0:
2573        return Integer(1)
2574
2575    try:
2576        w = v[0]
2577        mpz_init_set(z, w.value)
2578
2579        _sig_on
2580        for i from 1 <= i < n:
2581            w = v[i]
2582            mpz_gcd(z, z, w.value)
2583            if mpz_cmp_si(z, 1) == 0:
2584                _sig_off
2585                return Integer(1)
2586        _sig_off
2587    except TypeError:
2588        w = Integer(v[0])
2589        mpz_init_set(z, w.value)
2590
2591        _sig_on
2592        for i from 1 <= i < n:
2593            w = Integer(v[i])
2594            mpz_gcd(z, z, w.value)
2595            if mpz_cmp_si(z, 1) == 0:
2596                _sig_off
2597                return Integer(1)
2598        _sig_off
2599
2600
2601    w = PY_NEW(Integer)
2602    mpz_set(w.value, z)
2603    mpz_clear(z)
2604    return w
2605
2606def make_integer(s):
2607    cdef Integer r
2608    r = PY_NEW(Integer)
2609    r._reduce_set(s)
2610    return r
2611
2612# TODO: where is this used? Never doctested. Should I call IntegerRing.random_element()?
2613from random import randint
2614def random_integer(min=-2, max=2):
2615    cdef Integer x
2616    cdef int _min, _max, r
2617    try:
2618        _min = min
2619        _max = max
2620        x = PY_NEW(Integer)
2621        r = randint() % (_max - _min + 1) + _min
2622        mpz_set_si(x.value, r)
2623        return x
2624    except OverflowError:
2625        return Integer(randint(min,max))
2626
2627
2628############### INTEGER CREATION CODE #####################
2629######## There is nothing to see here, move along   #######
2630
2631cdef extern from *:
2632
2633    ctypedef struct RichPyObject "PyObject"
2634
2635    # We need a PyTypeObject with elements so we can
2636    # get and set tp_new, tp_dealloc, tp_flags, and tp_basicsize
2637    ctypedef struct RichPyTypeObject "PyTypeObject":
2638
2639        # We replace this one
2640        PyObject*      (*    tp_new) ( RichPyTypeObject*, PyObject*, PyObject*)
2641
2642        # Not used, may be useful to determine correct memory management function
2643        RichPyObject *(*   tp_alloc) ( RichPyTypeObject*, size_t )
2644
2645        # We replace this one
2646        void           (*tp_dealloc) ( PyObject*)
2647
2648        # Not used, may be useful to determine correct memory management function
2649        void          (*    tp_free) ( PyObject* )
2650
2651        # sizeof(Object)
2652        size_t tp_basicsize
2653
2654        # We set a flag here to circumvent the memory manager
2655        long tp_flags
2656
2657    cdef long Py_TPFLAGS_HAVE_GC
2658
2659    # We need a PyObject where we can get/set the refcnt directly
2660    # and access the type.
2661    ctypedef struct RichPyObject "PyObject":
2662        int ob_refcnt
2663        RichPyTypeObject* ob_type
2664
2665    # Allocation
2666    RichPyObject* PyObject_MALLOC(int)
2667
2668    # Useful for debugging, see below
2669    void PyObject_INIT(RichPyObject *, RichPyTypeObject *)
2670
2671    # Free
2672    void PyObject_FREE(PyObject*)
2673
2674
2675# We need a couple of internal GMP datatypes.
2676
2677# This may be potentialy very dangerous as it reaches
2678# deeply into the internal structure of GMP which may not
2679# be consistant across future versions of GMP.
2680# See extensive note in the fast_tp_new() function below.
2681
2682cdef extern from "gmp.h":
2683    ctypedef void* mp_ptr #"mp_ptr"
2684
2685    # We allocate _mp_d directly (mpz_t is typedef of this in GMP)
2686    ctypedef struct __mpz_struct "__mpz_struct":
2687        mp_ptr _mp_d
2688        size_t _mp_alloc
2689        size_t _mp_size
2690
2691    # sets the three free, alloc, and realloc function pointers to the
2692    # memory management functions set in GMP. Accepts NULL pointer.
2693    # Potentially dangerous if changed by calling
2694    # mp_set_memory_functions again after we initialized this module.
2695    void mp_get_memory_functions (void *(**alloc) (size_t), void *(**realloc)(void *, size_t, size_t), void (**free) (void *, size_t))
2696
2697    # GMP's configuration of how many Bits are stuffed into a limb
2698    cdef int __GMP_BITS_PER_MP_LIMB
2699
2700# This variable holds the size of any Integer object in bytes.
2701cdef int sizeof_Integer
2702
2703# We use a global Integer element to steal all the references
2704# from.  DO NOT INITIALIZE IT AGAIN and DO NOT REFERENCE IT!
2705cdef Integer global_dummy_Integer
2706global_dummy_Integer = Integer()
2707
2708# Accessing the .value attribute of an Integer object causes Pyrex to
2709# refcount it. This is problematic, because that causes overhead and
2710# more importantly an infinite loop in the destructor. If you refcount
2711# in the destructor and the refcount reaches zero (which is true
2712# everytime) the destructor is called.
2713#
2714# To avoid this we calculate the byte offset of the value member and
2715# remember it in this variable.
2716#
2717# Eventually this may be rendered obsolete by a change in SageX allowing
2718# non-reference counted extension types.
2719cdef long mpz_t_offset
2720
2721
2722# stores the GMP alloc function
2723cdef void * (* mpz_alloc)(size_t)
2724
2725# stores the GMP free function
2726cdef void (* mpz_free)(void *, size_t)
2727
2728# A global  pool for performance when integers are rapidly created and destroyed.
2729# It operates on the following principles:
2730#
2731# - The pool starts out empty.
2732# - When an new integer is needed, one from the pool is returned
2733#   if available, otherwise a new Integer object is created
2734# - When an integer is collected, it will add it to the pool
2735#   if there is room, otherwise it will be deallocated.
2736
2737cdef enum:
2738    integer_pool_size = 100 # Pyrex has no way of defining constants
2739
2740cdef PyObject* integer_pool[integer_pool_size]
2741cdef int integer_pool_count = 0
2742
2743# used for profiling the pool
2744cdef int total_alloc = 0
2745cdef int use_pool = 0
2746
2747# The signature of tp_new is
2748# PyObject* tp_new(RichPyTypeObject *t, PyObject *a, PyObject *k).
2749# However we only use t in this implementation.
2750#
2751# t in this case is the Integer TypeObject.
2752
2753cdef PyObject* fast_tp_new(RichPyTypeObject *t, PyObject *a, PyObject *k):
2754
2755    global integer_pool, integer_pool_count, total_alloc, use_pool
2756
2757    cdef RichPyObject* new
2758
2759    # for profiling pool usage
2760    # total_alloc += 1
2761
2762    # If there is a ready integer in the pool, we will
2763    # decrement the counter and return that.
2764
2765    if integer_pool_count > 0:
2766
2767        # for profiling pool usage
2768        # use_pool += 1
2769
2770        integer_pool_count -= 1
2771        new = <RichPyObject *> integer_pool[integer_pool_count]
2772
2773    # Otherwise, we have to create one.
2774
2775    else:
2776
2777        # allocate enough room for the Integer, sizeof_Integer is
2778        # sizeof(Integer). The use of PyObject_MALLOC directly
2779        # assumes that Integers are not garbage collected, i.e.
2780        # they do not pocess references to other Python
2781        # objects (Aas indicated by the Py_TPFLAGS_HAVE_GC flag).
2782        # See below for a more detailed description.
2783
2784        new = PyObject_MALLOC( sizeof_Integer )
2785
2786        # Now set every member as set in z, the global dummy Integer
2787        # created before this tp_new started to operate.
2788
2789        memcpy(new, (<void*>global_dummy_Integer), sizeof_Integer )
2790
2791        # This line is only needed if Python is compiled in debugging
2792        # mode './configure --with-pydebug'. If that is the case a Python
2793        # object has a bunch of debugging fields which are initialized
2794        # with this macro. For speed reasons, we don't call it if Python
2795        # is not compiled in debug mode. So uncomment the following line
2796        # if you are debugging Python.
2797
2798        #PyObject_INIT(new, (<RichPyObject*>global_dummy_Integer).ob_type)
2799
2800        # We take the address 'new' and move mpz_t_offset bytes (chars)
2801        # to the address of 'value'. We treat that address as a pointer
2802        # to a mpz_t struct and allocate memory for the _mp_d element of
2803        # that struct. We allocate one limb.
2804        #
2805        # What is done here is potentialy very dangerous as it reaches
2806        # deeply into the internal structure of GMP. Consequently things
2807        # may break if a new release of GMP changes some internals. To
2808        # emphazise this, this is what the GMP manual has to say about
2809        # the documentation for the struct we are using:
2810        #
2811        #  "This chapter is provided only for informational purposes and the
2812        #  various internals described here may change in future GMP releases.
2813        #  Applications expecting to be compatible with future releases should use
2814        #  only the documented interfaces described in previous chapters."
2815        #
2816        # If this line is used SAGE is not such an application.
2817        #
2818        # The clean version of the following line is:
2819        #
2820        #  mpz_init(( <mpz_t>(<char *>new + mpz_t_offset) )
2821        #
2822        # We save time both by avoiding an extra function call and
2823        # because the rest of the mpz struct was already initalized
2824        # fully using the memcpy above.
2825
2826        (<__mpz_struct *>( <char *>new + mpz_t_offset) )._mp_d = <mp_ptr>mpz_alloc(__GMP_BITS_PER_MP_LIMB >> 3)
2827
2828    # The global_dummy_Integer may have a reference count larger than
2829    # one, but it is expected that newly created objects have a
2830    # reference count of one. This is potentially unneeded if
2831    # everybody plays nice, because the gobal_dummy_Integer has only
2832    # one reference in that case.
2833
2834    # Objects from the pool have reference count zero, so this
2835    # needs to be set in this case.
2836
2837    new.ob_refcnt = 1
2838
2839    return new
2840
2841cdef void fast_tp_dealloc(PyObject* o):
2842
2843    # If there is room in the pool for a used integer object,
2844    # then put it in rather than deallocating it.
2845
2846    global integer_pool, integer_pool_count
2847
2848    if integer_pool_count < integer_pool_size:
2849
2850        # Here we free any extra memory used by the mpz_t by
2851        # setting it to a single limb.
2852        if (<__mpz_struct *>( <char *>o + mpz_t_offset))._mp_alloc > 1:
2853            _mpz_realloc(<mpz_t *>( <char *>o + mpz_t_offset), 1)
2854
2855        # It's cheap to zero out an integer, so do it here.
2856        (<__mpz_struct *>( <char *>o + mpz_t_offset))._mp_size = 0
2857
2858        # And add it to the pool.
2859        integer_pool[integer_pool_count] = o
2860        integer_pool_count += 1
2861        return
2862
2863    # Again, we move to the mpz_t and clear it. See above, why this is evil.
2864    # The clean version of this line would be:
2865    #   mpz_clear(<mpz_t>(<char *>o + mpz_t_offset))
2866
2867    mpz_free((<__mpz_struct *>( <char *>o + mpz_t_offset) )._mp_d, 0)
2868
2869    # Free the object. This assumes that Py_TPFLAGS_HAVE_GC is not
2870    # set. If it was set another free function would need to be
2871    # called.
2872
2873    PyObject_FREE(o)
2874
2875hook_fast_tp_functions()
2876
2877def hook_fast_tp_functions():
2878    """
2879    """
2880    global global_dummy_Integer, mpz_t_offset, sizeof_Integer
2881
2882    cdef long flag
2883
2884    cdef RichPyObject* o
2885    o = <RichPyObject*>global_dummy_Integer
2886
2887    # By default every object created in Pyrex is garbage
2888    # collected. This means it may have references to other objects
2889    # the Garbage collector has to look out for. We remove this flag
2890    # as the only reference an Integer has is to the global Integer
2891    # ring. As this object is unique we don't need to garbage collect
2892    # it as we always have a module level reference to it. If another
2893    # attribute is added to the Integer class this flag removal so as
2894    # the alloc and free functions may not be used anymore.
2895    # This object will still be reference counted.
2896    flag = Py_TPFLAGS_HAVE_GC
2897    o.ob_type.tp_flags = <long>(o.ob_type.tp_flags & (~flag))
2898
2899    # calculate the offset of the GMP mpz_t to avoid casting to/from
2900    # an Integer which includes reference counting. Reference counting
2901    # is bad in constructors and destructors as it potentially calls
2902    # the destructor.
2903    # Eventually this may be rendered obsolete by a change in SageX allowing
2904    # non-reference counted extension types.
2905    mpz_t_offset = <char *>(&global_dummy_Integer.value) - <char *>o
2906
2907    # store how much memory needs to be allocated for an Integer.
2908    sizeof_Integer = o.ob_type.tp_basicsize
2909
2910    # get the functions to do memory management for the GMP elements
2911    # WARNING: if the memory management functions are changed after
2912    # this initialisation, we are/you are doomed.
2913
2914    mp_get_memory_functions(&mpz_alloc, NULL, &mpz_free)
2915
2916    # Finally replace the functions called when an Integer needs
2917    # to be constructed/destructed.
2918    o.ob_type.tp_new = &fast_tp_new
2919    o.ob_type.tp_dealloc = &fast_tp_dealloc
2920
2921def time_alloc_list(n):
2922    cdef int i
2923    l = []
2924    for i from 0 <= i < n:
2925        l.append(PY_NEW(Integer))
2926
2927    return l
2928
2929def time_alloc(n):
2930    cdef int i
2931    for i from 0 <= i < n:
2932        z = PY_NEW(Integer)
2933
2934def pool_stats():
2935    print "Used pool %s / %s times" % (use_pool, total_alloc)
2936    print "Pool contains %s / %s items" % (integer_pool_count, integer_pool_size)
2937
2938cdef integer(x):
2939    if PY_TYPE_CHECK(x, Integer):
2940        return x
2941    return Integer(x)
Note: See TracBrowser for help on using the repository browser.