Context Navigation

source: sage/rings/integer.pyx @ 5781:0a837401a64f

Visit:

Revision 5781:0a837401a64f, 88.5 KB checked in by William Stein <wstein@…>, 6 years ago (diff)
David Harvey: add new mpz_get_pyintlong() function which returns either python int (fast!) or python long if it doesn't fit; change some Integer methods to use this new function

Line
1	r"""
2	Elements of the ring $\Z$ of integers
3
4	AUTHORS:
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
20	EXAMPLES:
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
49	Multiplication:
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
64	Coercions:
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	#
87	# http://www.gnu.org/licenses/
88	#*****************************************************************************
89
90	doc="""
91	Integers
92	"""
93
94	import operator
95
96	import sys
97
98	include "../ext/gmp.pxi"
99	include "../ext/interrupt.pxi" # ctrl-c interrupt block support
100	include "../ext/stdsage.pxi"
101	include "../ext/python_list.pxi"
102
103	cdef extern from "../ext/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
110	from sage.libs.pari.gen cimport gen as pari_gen
111
112	cdef class Integer(sage.structure.element.EuclideanDomainElement)
113
114	import sage.rings.infinity
115	import sage.libs.pari.all
116
117	cdef mpz_t mpz_tmp
118	mpz_init(mpz_tmp)
119
120	cdef public int set_mpz(Integer self, mpz_t value):
121	mpz_set(self.value, value)
122
123	cdef set_from_Integer(Integer self, Integer other):
124	mpz_set(self.value, other.value)
125
126	cdef set_from_int(Integer self, int other):
127	mpz_set_si(self.value, other)
128
129	cdef public mpz_t* get_value(Integer self):
130	return &self.value
131
132	MAX_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
139	from sage.structure.sage_object cimport SageObject
140	from sage.structure.element cimport EuclideanDomainElement, ModuleElement
141	from sage.structure.element import bin_op
142
143	import integer_ring
144	the_integer_ring = integer_ring.ZZ
145
146	cdef 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
151	set_zero_one_elements()
152
153	def is_Integer(x):
154	return PY_TYPE_CHECK(x, Integer)
155
156	cdef 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, ntl_c_ZZ *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 \code{self % 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	"""
1034	cdef Integer _modulus, _self
1035	_modulus = integer(modulus)
1036	if not _modulus:
1037	raise ZeroDivisionError, "Integer modulo by zero"
1038	_self = integer(self)
1039
1040	cdef Integer x
1041	x = PY_NEW(Integer)
1042
1043	_sig_on
1044	mpz_mod(x.value, _self.value, _modulus.value)
1045	_sig_off
1046
1047	return x
1048
1049
1050	def quo_rem(self, other):
1051	"""
1052	Returns the quotient and the remainder of
1053	self divided by other.
1054
1055	INPUT:
1056	other -- the integer the divisor
1057
1058	OUTPUT:
1059	q -- the quotient of self/other
1060	r -- the remainder of self/other
1061
1062	EXAMPLES:
1063	sage: z = Integer(231)
1064	sage: z.quo_rem(2)
1065	(115, 1)
1066	sage: z.quo_rem(-2)
1067	(-115, 1)
1068	sage: z.quo_rem(0)
1069	Traceback (most recent call last):
1070	...
1071	ZeroDivisionError: other (=0) must be nonzero
1072	"""
1073	cdef Integer _other, _self
1074	_other = integer(other)
1075	if not _other:
1076	raise ZeroDivisionError, "other (=%s) must be nonzero"%other
1077	_self = integer(self)
1078
1079	cdef Integer q, r
1080	q = PY_NEW(Integer)
1081	r = PY_NEW(Integer)
1082
1083	_sig_on
1084	if mpz_sgn(_other.value) == 1:
1085	mpz_fdiv_qr(q.value, r.value, _self.value, _other.value)
1086	else:
1087	mpz_cdiv_qr(q.value, r.value, _self.value, _other.value)
1088	_sig_off
1089
1090	return q, r
1091
1092	def div(self, other):
1093	"""
1094	Returns the quotient of self divided by other.
1095
1096	INPUT:
1097	other -- the integer the divisor
1098
1099	OUTPUT:
1100	q -- the quotient of self/other
1101
1102	EXAMPLES:
1103	sage: z = Integer(231)
1104	sage: z.div(2)
1105	115
1106	sage: z.div(-2)
1107	-115
1108	sage: z.div(0)
1109	Traceback (most recent call last):
1110	...
1111	ZeroDivisionError: other (=0) must be nonzero
1112	"""
1113	cdef Integer _other, _self
1114	_other = integer(other)
1115	if not _other:
1116	raise ZeroDivisionError, "other (=%s) must be nonzero"%other
1117	_self = integer(self)
1118
1119	cdef Integer q, r
1120	q = PY_NEW(Integer)
1121	r = PY_NEW(Integer)
1122
1123	_sig_on
1124	mpz_tdiv_qr(q.value, r.value, _self.value, _other.value)
1125	_sig_off
1126
1127	return q
1128
1129
1130	def powermod(self, exp, mod):
1131	"""
1132	Compute self**exp modulo mod.
1133
1134	EXAMPLES:
1135	sage: z = 2
1136	sage: z.powermod(31,31)
1137	2
1138	sage: z.powermod(0,31)
1139	1
1140	sage: z.powermod(-31,31) == 2^-31 % 31
1141	True
1142
1143	As expected, the following is invalid:
1144	sage: z.powermod(31,0)
1145	Traceback (most recent call last):
1146	...
1147	ZeroDivisionError: cannot raise to a power modulo 0
1148	"""
1149	cdef Integer x, _exp, _mod
1150	_exp = Integer(exp); _mod = Integer(mod)
1151	if mpz_cmp_si(_mod.value,0) == 0:
1152	raise ZeroDivisionError, "cannot raise to a power modulo 0"
1153
1154	x = PY_NEW(Integer)
1155
1156	_sig_on
1157	mpz_powm(x.value, self.value, _exp.value, _mod.value)
1158	_sig_off
1159
1160	return x
1161
1162	def rational_reconstruction(self, Integer m):
1163	import rational
1164	return rational.pyrex_rational_reconstruction(self, m)
1165
1166	def powermodm_ui(self, exp, mod):
1167	r"""
1168	Computes self**exp modulo mod, where exp is an unsigned
1169	long integer.
1170
1171	EXAMPLES:
1172	sage: z = 32
1173	sage: z.powermodm_ui(2, 4)
1174	0
1175	sage: z.powermodm_ui(2, 14)
1176	2
1177	sage: z.powermodm_ui(2^32-2, 14)
1178	2
1179	sage: z.powermodm_ui(2^32-1, 14)
1180	Traceback (most recent call last): # 32-bit
1181	... # 32-bit
1182	OverflowError: exp (=4294967295) must be <= 4294967294 # 32-bit
1183	8 # 64-bit
1184	sage: z.powermodm_ui(2^65, 14)
1185	Traceback (most recent call last):
1186	...
1187	OverflowError: exp (=36893488147419103232) must be <= 4294967294 # 32-bit
1188	OverflowError: exp (=36893488147419103232) must be <= 18446744073709551614 # 64-bit
1189	"""
1190	if exp < 0:
1191	raise ValueError, "exp (=%s) must be nonnegative"%exp
1192	elif exp > MAX_UNSIGNED_LONG:
1193	raise OverflowError, "exp (=%s) must be <= %s"%(exp, MAX_UNSIGNED_LONG)
1194	cdef Integer x, _mod
1195	_mod = Integer(mod)
1196	x = PY_NEW(Integer)
1197
1198	_sig_on
1199	mpz_powm_ui(x.value, self.value, exp, _mod.value)
1200	_sig_off
1201
1202	return x
1203
1204	def __int__(self):
1205	return mpz_get_pyintlong(self.value)
1206
1207	def __long__(self):
1208	return mpz_get_pylong(self.value)
1209
1210	def __float__(self):
1211	return mpz_get_d(self.value)
1212
1213	def __hash__(self):
1214	return mpz_pythonhash(self.value)
1215
1216	def factor(self, algorithm='pari'):
1217	"""
1218	Return the prime factorization of the integer as a list of
1219	pairs $(p,e)$, where $p$ is prime and $e$ is a positive integer.
1220
1221	INPUT:
1222	algorithm -- string
1223	* 'pari' -- (default) use the PARI c library
1224	* 'kash' -- use KASH computer algebra system (requires
1225	the optional kash package be installed)
1226	"""
1227	import sage.rings.integer_ring
1228	return sage.rings.integer_ring.factor(self, algorithm=algorithm)
1229
1230	def coprime_integers(self, m):
1231	"""
1232	Return the positive integers $< m$ that are coprime to self.
1233
1234	EXAMPLES:
1235	sage: n = 8
1236	sage: n.coprime_integers(8)
1237	[1, 3, 5, 7]
1238	sage: n.coprime_integers(11)
1239	[1, 3, 5, 7, 9]
1240	sage: n = 5; n.coprime_integers(10)
1241	[1, 2, 3, 4, 6, 7, 8, 9]
1242	sage: n.coprime_integers(5)
1243	[1, 2, 3, 4]
1244	sage: n = 99; n.coprime_integers(99)
1245	[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]
1246
1247	AUTHORS:
1248	-- Naqi Jaffery (2006-01-24): examples
1249
1250	ALGORITHM: Naive -- compute lots of GCD's. If this isn't good
1251	enough for you, please code something better and submit a
1252	patch.
1253	"""
1254	# TODO -- make VASTLY faster
1255	v = []
1256	for n in range(1,m):
1257	if self.gcd(n) == 1:
1258	v.append(Integer(n))
1259	return v
1260
1261	def divides(self, n):
1262	"""
1263	Return True if self divides n.
1264
1265	EXAMPLES:
1266	sage: Z = IntegerRing()
1267	sage: Z(5).divides(Z(10))
1268	True
1269	sage: Z(0).divides(Z(5))
1270	False
1271	sage: Z(10).divides(Z(5))
1272	False
1273	"""
1274	cdef bint t
1275	cdef Integer _n
1276	_n = Integer(n)
1277	if mpz_sgn(self.value) == 0:
1278	return mpz_sgn(_n.value) == 0
1279	_sig_on
1280	t = mpz_divisible_p(_n.value, self.value)
1281	_sig_off
1282	return t
1283
1284	cdef RingElement _valuation(Integer self, Integer p):
1285	r"""
1286	Return the p-adic valuation of self.
1287
1288	We do not require that p be prime, but it must be at least 2.
1289	For more documentation see \code{valuation}
1290
1291	AUTHOR:
1292	-- David Roe (3/31/07)
1293	"""
1294	if mpz_sgn(self.value) == 0:
1295	return sage.rings.infinity.infinity
1296	if mpz_cmp_ui(p.value, 2) < 0:
1297	raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
1298	cdef Integer v
1299	v = PY_NEW(Integer)
1300	cdef mpz_t u
1301	mpz_init(u)
1302	_sig_on
1303	mpz_set_ui(v.value, mpz_remove(u, self.value, p.value))
1304	_sig_off
1305	mpz_clear(u)
1306	return v
1307
1308	cdef object _val_unit(Integer self, Integer p):
1309	r"""
1310	Returns a pair: the p-adic valuation of self, and the p-adic unit of self.
1311
1312	We do not require the p be prime, but it must be at least 2.
1313	For more documentation see \code{val_unit}
1314
1315	AUTHOR:
1316	-- David Roe (3/31/07)
1317	"""
1318	cdef Integer v, u
1319	if mpz_cmp_ui(p.value, 2) < 0:
1320	raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
1321	if self == 0:
1322	u = ONE
1323	Py_INCREF(ONE)
1324	return (sage.rings.infinity.infinity, u)
1325	v = PY_NEW(Integer)
1326	u = PY_NEW(Integer)
1327	_sig_on
1328	mpz_set_ui(v.value, mpz_remove(u.value, self.value, p.value))
1329	_sig_off
1330	return (v, u)
1331
1332	def valuation(self, p):
1333	"""
1334	Return the p-adic valuation of self.
1335
1336	INPUT:
1337	p -- an integer at least 2.
1338
1339	EXAMPLE:
1340	sage: n = 60
1341	sage: n.valuation(2)
1342	2
1343	sage: n.valuation(3)
1344	1
1345	sage: n.valuation(7)
1346	0
1347	sage: n.valuation(1)
1348	Traceback (most recent call last):
1349	...
1350	ValueError: You can only compute the valuation with respect to a integer larger than 1.
1351
1352	We do not require that p is a prime:
1353	sage: (2^11).valuation(4)
1354	5
1355	"""
1356	return self._valuation(Integer(p))
1357
1358	def val_unit(self, p):
1359	r"""
1360	Returns a pair: the p-adic valuation of self, and the p-adic unit of self.
1361
1362	INPUT:
1363	p -- an integer at least 2.
1364
1365	OUTPUT:
1366	v_p(self) -- the p-adic valuation of self
1367	u_p(self) -- self / p^{v_p(self)}
1368
1369	EXAMPLE:
1370	sage: n = 60
1371	sage: n.val_unit(2)
1372	(2, 15)
1373	sage: n.val_unit(3)
1374	(1, 20)
1375	sage: n.val_unit(7)
1376	(0, 60)
1377	sage: (2^11).val_unit(4)
1378	(5, 2)
1379	sage: 0.val_unit(2)
1380	(+Infinity, 1)
1381	"""
1382	return self._val_unit(Integer(p))
1383
1384	def ord(self, p=None):
1385	"""Synonym for valuation
1386
1387	EXAMPLES:
1388	sage: n=12
1389	sage: n.ord(3)
1390	1
1391	"""
1392	return self.valuation(p)
1393
1394	cdef Integer _divide_knowing_divisible_by(Integer self, Integer right):
1395	r"""
1396	Returns the integer self / right when self is divisible by right.
1397
1398	If self is not divisible by right, the return value is undefined, but seems to be close to self/right.
1399	For more documentation see \code{divide_knowing_divisible_by}
1400
1401	AUTHOR:
1402	-- David Roe (3/31/07)
1403	"""
1404	if mpz_cmp_ui(right.value, 0) == 0:
1405	raise ZeroDivisionError
1406	cdef Integer x
1407	x = PY_NEW(Integer)
1408	if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
1409	# Only use the signal handler (to enable ctrl-c out) when the
1410	# quotient might take a while to compute
1411	_sig_on
1412	mpz_divexact(x.value, self.value, right.value)
1413	_sig_off
1414	else:
1415	mpz_divexact(x.value, self.value, right.value)
1416	return x
1417
1418	def divide_knowing_divisible_by(self, right):
1419	r"""
1420	Returns the integer self / right when self is divisible by right.
1421
1422	If self is not divisible by right, the return value is undefined, but seems to be close to self/right.
1423
1424	EXAMPLES:
1425	sage: a = 8; b = 4
1426	sage: a.divide_knowing_divisible_by(b)
1427	2
1428	"""
1429	return self._divide_knowing_divisible_by(right)
1430
1431	def _lcm(self, Integer n):
1432	"""
1433	Returns the least common multiple of self and $n$.
1434
1435	EXAMPLES:
1436	sage: n = 60
1437	sage: n._lcm(150)
1438	300
1439	"""
1440	cdef mpz_t x
1441
1442	mpz_init(x)
1443
1444	_sig_on
1445	mpz_lcm(x, self.value, n.value)
1446	_sig_off
1447
1448
1449	cdef Integer z
1450	z = PY_NEW(Integer)
1451	mpz_set(z.value,x)
1452	mpz_clear(x)
1453	return z
1454
1455	def denominator(self):
1456	"""
1457	Return the denominator of this integer.
1458
1459	EXAMPLES:
1460	sage: x = 5
1461	sage: x.denominator()
1462	1
1463	sage: x = 0
1464	sage: x.denominator()
1465	1
1466	"""
1467	return ONE
1468
1469	def numerator(self):
1470	"""
1471	Return the numerator of this integer.
1472
1473	EXAMPLE:
1474	sage: x = 5
1475	sage: x.numerator()
1476	5
1477
1478	sage: x = 0
1479	sage: x.numerator()
1480	0
1481	"""
1482	return self
1483
1484	def factorial(self):
1485	"""
1486	Return the factorial $n!=1 \\cdot 2 \\cdot 3 \\cdots n$.
1487	Self must fit in an \\code{unsigned long int}.
1488
1489	EXAMPLES:
1490	sage: for n in srange(7):
1491	... print n, n.factorial()
1492	0 1
1493	1 1
1494	2 2
1495	3 6
1496	4 24
1497	5 120
1498	6 720
1499	"""
1500	if mpz_sgn(self.value) < 0:
1501	raise ValueError, "factorial -- self = (%s) must be nonnegative"%self
1502
1503	if mpz_cmp_ui(self.value,4294967295) > 0:
1504	raise ValueError, "factorial not implemented for n >= 2^32.\nThis is probably OK, since the answer would have billions of digits."
1505
1506	cdef unsigned int n
1507	n = self
1508
1509	cdef mpz_t x
1510	cdef Integer z
1511
1512	mpz_init(x)
1513
1514	_sig_on
1515	mpz_fac_ui(x, n)
1516	_sig_off
1517
1518	z = PY_NEW(Integer)
1519	set_mpz(z, x)
1520	mpz_clear(x)
1521	return z
1522
1523	def floor(self):
1524	"""
1525	Return the floor of self, which is just self since self is an integer.
1526
1527	EXAMPLES:
1528	sage: n = 6
1529	sage: n.floor()
1530	6
1531	"""
1532	return self
1533
1534	def ceil(self):
1535	"""
1536	Return the ceiling of self, which is self since self is an integer.
1537
1538	EXAMPLES:
1539	sage: n = 6
1540	sage: n.ceil()
1541	6
1542	"""
1543	return self
1544
1545	def is_one(self):
1546	r"""
1547	Returns \code{True} if the integers is $1$, otherwise \code{False}.
1548
1549	EXAMPLES:
1550	sage: Integer(1).is_one()
1551	True
1552	sage: Integer(0).is_one()
1553	False
1554	"""
1555	return mpz_cmp_si(self.value, 1) == 0
1556
1557	def __nonzero__(self):
1558	r"""
1559	Returns \code{True} if the integers is not $0$, otherwise \code{False}.
1560
1561	EXAMPLES:
1562	sage: Integer(1).is_zero()
1563	False
1564	sage: Integer(0).is_zero()
1565	True
1566	"""
1567	return mpz_sgn(self.value) != 0
1568
1569	def is_unit(self):
1570	r"""
1571	Returns \code{true} if this integer is a unit, i.e., 1 or $-1$.
1572
1573	EXAMPLES:
1574	sage: for n in srange(-2,3):
1575	... print n, n.is_unit()
1576	-2 False
1577	-1 True
1578	0 False
1579	1 True
1580	2 False
1581	"""
1582	return mpz_cmp_si(self.value, -1) == 0 or mpz_cmp_si(self.value, 1) == 0
1583
1584	def is_square(self):
1585	r"""
1586	Returns \code{True} if self is a perfect square
1587
1588	EXAMPLES:
1589	sage: Integer(4).is_square()
1590	True
1591	sage: Integer(41).is_square()
1592	False
1593	"""
1594	return mpz_perfect_square_p(self.value)
1595
1596	def is_power(self):
1597	r"""
1598	Returns \code{True} if self is a perfect power, ie if there
1599	exist integers a and b, $b > 1$ with $self = a^b$.
1600
1601	EXAMPLES:
1602	sage: Integer(-27).is_power()
1603	True
1604	sage: Integer(12).is_power()
1605	False
1606	"""
1607	return mpz_perfect_power_p(self.value)
1608
1609	cdef bint _is_power_of(Integer self, Integer n):
1610	r"""
1611	Returns a non-zero int if there is an integer b with $self = n^b$.
1612
1613	For more documentation see \code{is_power_of}
1614
1615	AUTHOR:
1616	-- David Roe (3/31/07)
1617	"""
1618	cdef int a
1619	cdef unsigned long b, c
1620	cdef mpz_t u, sabs, nabs
1621	a = mpz_cmp_ui(n.value, 2)
1622	if a <= 0: # n <= 2
1623	if a == 0: # n == 2
1624	if mpz_popcount(self.value) == 1: #number of bits set in self == 1
1625	return 1
1626	else:
1627	return 0
1628	a = mpz_cmp_si(n.value, -2)
1629	if a >= 0: # -2 <= n < 2:
1630	a = mpz_get_si(n.value)
1631	if a == 1: # n == 1
1632	if mpz_cmp_ui(self.value, 1) == 0: # Only 1 is a power of 1
1633	return 1
1634	else:
1635	return 0
1636	elif a == 0: # n == 0
1637	if mpz_cmp_ui(self.value, 0) == 0 or mpz_cmp_ui(self.value, 1) == 0: # 0^0 = 1, 0^x = 0
1638	return 1
1639	else:
1640	return 0
1641	elif a == -1: # n == -1
1642	if mpz_cmp_ui(self.value, 1) == 0 or mpz_cmp_si(self.value, -1) == 0: # 1 and -1 are powers of -1
1643	return 1
1644	else:
1645	return 0
1646	elif a == -2: # n == -2
1647	mpz_init(sabs)
1648	mpz_abs(sabs, self.value)
1649	if mpz_popcount(sabs) == 1: # number of bits set in \|self\| == 1
1650	b = mpz_scan1(sabs, 0) % 2 # b == 1 if \|self\| is an odd power of 2, 0 if \|self\| is an even power
1651	mpz_clear(sabs)
1652	if (b == 1 and mpz_cmp_ui(self.value, 0) < 0) or (b == 0 and mpz_cmp_ui(self.value, 0) > 0):
1653	# An odd power of -2 is negative, an even power must be positive.
1654	return 1
1655	else: # number of bits set in \|self\| is not 1, so self cannot be a power of -2
1656	return 0
1657	else: # \|self\| is not a power of 2, so self cannot be a power of -2
1658	return 0
1659	else: # n < -2
1660	mpz_init(nabs)
1661	mpz_neg(nabs, n.value)
1662	if mpz_popcount(nabs) == 1: # \|n\| = 2^k for k >= 2. We special case this for speed
1663	mpz_init(sabs)
1664	mpz_abs(sabs, self.value)
1665	if mpz_popcount(sabs) == 1: # \|self\| = 2^L for some L >= 0.
1666	b = mpz_scan1(sabs, 0) # the bit that self is set at
1667	c = mpz_scan1(nabs, 0) # the bit that n is set at
1668	# Having obtained b and c, we're done with nabs and sabs (on this branch anyway)
1669	mpz_clear(nabs)
1670	mpz_clear(sabs)
1671	if b % c == 0: # Now we know that \|self\| is a power of \|n\|
1672	b = (b // c) % 2 # Whether b // c is even or odd determines whether (-2^c)^(b // c) is positive or negative
1673	a = mpz_cmp_ui(self.value, 0)
1674	if b == 0 and a > 0 or b == 1 and a < 0:
1675	# These two cases are that b // c is even and self positive, or b // c is odd and self negative
1676	return 1
1677	else: # The sign of self is wrong
1678	return 0
1679	else: # Since \|self\| is not a power of \|n\|, self cannot be a power of n
1680	return 0
1681	else: # self is not a power of 2, and thus cannot be a power of n, which is a power of 2.
1682	mpz_clear(nabs)
1683	mpz_clear(sabs)
1684	return 0
1685	else: # \|n\| is not a power of 2, so we use mpz_remove
1686	mpz_init(u)
1687	_sig_on
1688	b = mpz_remove(u, self.value, nabs)
1689	_sig_off
1690	# Having obtained b and u, we're done with nabs
1691	mpz_clear(nabs)
1692	if mpz_cmp_ui(u, 1) == 0: # self is a power of \|n\|
1693	mpz_clear(u)
1694	if b % 2 == 0: # an even power of \|n\|, and since self > 0, this means that self is a power of n
1695	return 1
1696	else:
1697	return 0
1698	elif mpz_cmp_si(u, -1) == 0: # -self is a power of \|n\|
1699	mpz_clear(u)
1700	if b % 2 == 1: # an odd power of \|n\|, and thus self is a power of n
1701	return 1
1702	else:
1703	return 1
1704	else: # \|self\| is not a power of \|n\|, so self cannot be a power of n
1705	mpz_clear(u)
1706	return 0
1707	elif mpz_popcount(n.value) == 1: # n > 2 and in fact n = 2^k for k >= 2
1708	if mpz_popcount(self.value) == 1: # since n is a power of 2, so must self be.
1709	if mpz_scan1(self.value, 0) % mpz_scan1(n.value, 0) == 0: # log_2(self) is divisible by log_2(n)
1710	return 1
1711	else:
1712	return 0
1713	else: # self is not a power of 2, and thus not a power of n
1714	return 0
1715	else: # n > 2, but not a power of 2, so we use mpz_remove
1716	mpz_init(u)
1717	_sig_on
1718	mpz_remove(u, self.value, n.value)
1719	_sig_off
1720	a = mpz_cmp_ui(u, 1)
1721	mpz_clear(u)
1722	if a == 0:
1723	return 1
1724	else:
1725	return 0
1726
1727	def is_power_of(Integer self, n):
1728	r"""
1729	Returns \code{True} if there is an integer b with $self = n^b$.
1730
1731	EXAMPLES:
1732	sage: Integer(64).is_power_of(4)
1733	True
1734	sage: Integer(64).is_power_of(16)
1735	False
1736
1737	TESTS:
1738	sage: Integer(-64).is_power_of(-4)
1739	True
1740	sage: Integer(-32).is_power_of(-2)
1741	True
1742	sage: Integer(1).is_power_of(1)
1743	True
1744	sage: Integer(-1).is_power_of(-1)
1745	True
1746	sage: Integer(0).is_power_of(1)
1747	False
1748	sage: Integer(0).is_power_of(0)
1749	True
1750	sage: Integer(1).is_power_of(0)
1751	True
1752	sage: Integer(1).is_power_of(8)
1753	True
1754	sage: Integer(-8).is_power_of(2)
1755	False
1756
1757	NOTES:
1758	For large integers self, is_power_of() is faster than is_power(). The following examples gives some indication of how much faster.
1759	sage.: b = lcm(range(1,10000))
1760	sage.: b.exact_log(2)
1761	14446
1762	sage.: time for a in range(2, 1000): k = b.is_power()
1763	CPU times: user 1.68 s, sys: 0.01 s, total: 1.69 s
1764	Wall time: 1.71
1765	sage.: time for a in range(2, 1000): k = b.is_power_of(2)
1766	CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
1767	Wall time: 0.01
1768	sage.: time for a in range(2, 1000): k = b.is_power_of(3)
1769	CPU times: user 0.05 s, sys: 0.00 s, total: 0.05 s
1770	Wall time: 0.05
1771
1772	sage.: b = lcm(range(1, 1000))
1773	sage.: b.exact_log(2)
1774	1437
1775	sage.: time for a in range(2, 10000): k = b.is_power() # note that we change the range from the example above
1776	CPU times: user 0.40 s, sys: 0.00 s, total: 0.40 s
1777	Wall time: 0.40
1778	sage.: for a in range(2, 10000): k = b.is_power_of(TWO)
1779	CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
1780	Wall time: 0.03
1781	sage.: time for a in range(2, 10000): k = b.is_power_of(3)
1782	CPU times: user 0.11 s, sys: 0.01 s, total: 0.12 s
1783	Wall time: 0.13
1784	sage.: time for a in range(2, 10000): k = b.is_power_of(a)
1785	CPU times: user 0.08 s, sys: 0.01 s, total: 0.09 s
1786	Wall time: 0.10
1787	"""
1788	if not PY_TYPE_CHECK(n, Integer):
1789	n = Integer(n)
1790	return self._is_power_of(n)
1791
1792	def is_prime_power(self, flag=0):
1793	r"""
1794	Returns True if $x$ is a prime power, and False otherwise.
1795
1796	INPUT:
1797	flag (for primality testing) -- int
1798	0 (default): use a combination of algorithms.
1799	1: certify primality using the Pocklington-Lehmer Test.
1800	2: certify primality using the APRCL test.
1801
1802	EXAMPLES:
1803	sage: (-10).is_prime_power()
1804	False
1805	sage: (10).is_prime_power()
1806	False
1807	sage: (64).is_prime_power()
1808	True
1809	sage: (3^10000).is_prime_power()
1810	True
1811	sage: (10000).is_prime_power(flag=1)
1812	False
1813	"""
1814	if self.is_zero():
1815	return False
1816	elif self.is_one():
1817	return True
1818	elif mpz_sgn(self.value) < 0:
1819	return False
1820	if self.is_prime():
1821	return True
1822	if not self.is_perfect_power():
1823	return False
1824	k, g = self._pari_().ispower()
1825	if not k:
1826	raise RuntimeError, "inconsistent results between GMP and pari"
1827	return g.isprime(flag=flag)
1828
1829	def is_prime(self):
1830	r"""
1831	Retuns \code{True} if self is prime
1832
1833	EXAMPLES:
1834	sage: z = 2^31 - 1
1835	sage: z.is_prime()
1836	True
1837	sage: z = 2^31
1838	sage: z.is_prime()
1839	False
1840	"""
1841	return bool(self._pari_().isprime())
1842
1843	def is_pseudoprime(self):
1844	r"""
1845	Retuns \code{True} if self is a pseudoprime
1846
1847	EXAMPLES:
1848	sage: z = 2^31 - 1
1849	sage: z.is_pseudoprime()
1850	True
1851	sage: z = 2^31
1852	sage: z.is_pseudoprime()
1853	False
1854	"""
1855	return bool(self._pari_().ispseudoprime())
1856
1857	def is_perfect_power(self):
1858	r"""
1859	Retuns \code{True} if self is a perfect power.
1860
1861	EXAMPLES:
1862	sage: z = 8
1863	sage: z.is_perfect_power()
1864	True
1865	sage: z = 144
1866	sage: z.is_perfect_power()
1867	True
1868	sage: z = 10
1869	sage: z.is_perfect_power()
1870	False
1871	"""
1872	return mpz_perfect_power_p(self.value)
1873
1874	def jacobi(self, b):
1875	r"""
1876	Calculate the Jacobi symbol $\left(\frac{self}{b}\right)$.
1877
1878	EXAMPLES:
1879	sage: z = -1
1880	sage: z.jacobi(17)
1881	1
1882	sage: z.jacobi(19)
1883	-1
1884	sage: z.jacobi(17*19)
1885	-1
1886	sage: (2).jacobi(17)
1887	1
1888	sage: (3).jacobi(19)
1889	-1
1890	sage: (6).jacobi(17*19)
1891	-1
1892	sage: (6).jacobi(33)
1893	0
1894	sage: a = 3; b = 7
1895	sage: a.jacobi(b) == -b.jacobi(a)
1896	True
1897	"""
1898	cdef long tmp
1899	if PY_TYPE_CHECK(b, int):
1900	tmp = b
1901	if (tmp & 1) == 0:
1902	raise ValueError, "Jacobi symbol not defined for even b."
1903	return mpz_kronecker_si(self.value, tmp)
1904	if not PY_TYPE_CHECK(b, Integer):
1905	b = Integer(b)
1906	if mpz_even_p((<Integer>b).value):
1907	raise ValueError, "Jacobi symbol not defined for even b."
1908	return mpz_jacobi(self.value, (<Integer>b).value)
1909
1910	def kronecker(self, b):
1911	r"""
1912	Calculate the Kronecker symbol
1913	$\left(\frac{self}{b}\right)$ with the Kronecker extension
1914	$(self/2)=(2/self)$ when self odd, or $(self/2)=0$ when $self$ even.
1915
1916	EXAMPLES:
1917	EXAMPLES:
1918	sage: z = 5
1919	sage: z.kronecker(41)
1920	1
1921	sage: z.kronecker(43)
1922	-1
1923	sage: z.kronecker(8)
1924	-1
1925	sage: z.kronecker(15)
1926	0
1927	sage: a = 2; b = 5
1928	sage: a.kronecker(b) == b.kronecker(a)
1929	True
1930	"""
1931	if PY_TYPE_CHECK(b, int):
1932	return mpz_kronecker_si(self.value, b)
1933	if not PY_TYPE_CHECK(b, Integer):
1934	b = Integer(b)
1935	return mpz_kronecker(self.value, (<Integer>b).value)
1936
1937	def square_free_part(self):
1938	"""
1939	Return the square free part of $x$, i.e., a divisor z such that $x = z y^2$,
1940	for a perfect square $y^2$.
1941
1942	EXAMPLES:
1943	sage: square_free_part(100)
1944	1
1945	sage: square_free_part(12)
1946	3
1947	sage: square_free_part(173737)
1948	17
1949	sage: square_free_part(-17*32)
1950	-34
1951	sage: square_free_part(1)
1952	1
1953	sage: square_free_part(-1)
1954	-1
1955	sage: square_free_part(-2)
1956	-2
1957	sage: square_free_part(-4)
1958	-1
1959	"""
1960	if self.is_zero():
1961	return self
1962	F = self.factor()
1963	n = Integer(1)
1964	for p, e in F:
1965	if e % 2 != 0:
1966	n = n * p
1967	return n * F.unit()
1968
1969	def next_probable_prime(self):
1970	"""
1971	Returns the next probable prime after self, as determined by PARI.
1972
1973	EXAMPLES:
1974	sage: (-37).next_probable_prime()
1975	2
1976	sage: (100).next_probable_prime()
1977	101
1978	sage: (2^512).next_probable_prime()
1979	13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084171
1980	sage: 0.next_probable_prime()
1981	2
1982	sage: 126.next_probable_prime()
1983	127
1984	sage: 144168.next_probable_prime()
1985	144169
1986	"""
1987	return Integer( (self._pari_()+1).nextprime())
1988
1989	def next_prime(self, proof=True):
1990	r"""
1991	Returns the next prime after self.
1992
1993	INPUT:
1994	proof -- bool (default: True)
1995
1996	EXAMPLES:
1997	sage: Integer(100).next_prime()
1998	101
1999
2000	Use Proof = False, which is way faster:
2001	sage: b = (2^1024).next_prime(proof=False)
2002
2003	sage: Integer(0).next_prime()
2004	2
2005	sage: Integer(1001).next_prime()
2006	1009
2007	"""
2008	if self < 2: # negatives are not prime.
2009	return integer_ring.ZZ(2)
2010	if self == 2:
2011	return integer_ring.ZZ(3)
2012	if not proof:
2013	return Integer( (self._pari_()+1).nextprime())
2014	n = self
2015	if n % 2 == 0:
2016	n += 1
2017	else:
2018	n += 2
2019	while not n.is_prime(): # pari isprime is provably correct
2020	n += 2
2021	return integer_ring.ZZ(n)
2022
2023
2024	def additive_order(self):
2025	"""
2026	Return the additive order of self.
2027
2028	EXAMPLES:
2029	sage: ZZ(0).additive_order()
2030	1
2031	sage: ZZ(1).additive_order()
2032	+Infinity
2033	"""
2034	import sage.rings.infinity
2035	if self.is_zero():
2036	return Integer(1)
2037	else:
2038	return sage.rings.infinity.infinity
2039
2040	def multiplicative_order(self):
2041	r"""
2042	Return the multiplicative order of self.
2043
2044	EXAMPLES:
2045	sage: ZZ(1).multiplicative_order()
2046	1
2047	sage: ZZ(-1).multiplicative_order()
2048	2
2049	sage: ZZ(0).multiplicative_order()
2050	+Infinity
2051	sage: ZZ(2).multiplicative_order()
2052	+Infinity
2053	"""
2054	import sage.rings.infinity
2055	if mpz_cmp_si(self.value, 1) == 0:
2056	return Integer(1)
2057	elif mpz_cmp_si(self.value, -1) == 0:
2058	return Integer(2)
2059	else:
2060	return sage.rings.infinity.infinity
2061
2062	def is_squarefree(self):
2063	"""
2064	Returns True if this integer is not divisible by the square of
2065	any prime and False otherwise.
2066
2067	EXAMPLES:
2068	sage: Integer(100).is_squarefree()
2069	False
2070	sage: Integer(102).is_squarefree()
2071	True
2072	"""
2073	return self._pari_().issquarefree()
2074
2075	def _pari_(self):
2076	"""
2077	Returns the PARI version of this integer.
2078
2079	EXAMPLES:
2080	sage: n = 9390823
2081	sage: m = n._pari_(); m
2082	9390823
2083	sage: type(m)
2084	<type 'sage.libs.pari.gen.gen'>
2085
2086	ALGORITHM: Use base 10 Python string conversion, hence very
2087	very slow for large integers. If you can figure out how to
2088	input a number into PARI in hex, or otherwise optimize this,
2089	please implement it and send me a patch.
2090	"""
2091	#if self._pari is None:
2092	# better to do in hex, but I can't figure out
2093	# how to input/output a number in hex in PARI!!
2094	# TODO: (I could just think carefully about raw bytes and make this all much faster...)
2095	#self._pari = sage.libs.pari.all.pari(str(self))
2096	#return self._pari
2097	return sage.libs.pari.all.pari(str(self))
2098
2099	def _interface_init_(self):
2100	"""
2101	Return canonical string to coerce this integer to any other math
2102	software, i.e., just the string representation of this integer
2103	in base 10.
2104
2105	EXAMPLES:
2106	sage: n = 9390823
2107	sage: n._interface_init_()
2108	'9390823'
2109	"""
2110	return str(self)
2111
2112	def isqrt(self):
2113	r"""
2114	Returns the integer floor of the square root of self, or raises
2115	an \exception{ValueError} if self is negative.
2116
2117	EXAMPLE:
2118	sage: a = Integer(5)
2119	sage: a.isqrt()
2120	2
2121
2122	sage: Integer(-102).isqrt()
2123	Traceback (most recent call last):
2124	...
2125	ValueError: square root of negative integer not defined.
2126	"""
2127	if mpz_sgn(self.value) < 0:
2128	raise ValueError, "square root of negative integer not defined."
2129	cdef Integer x
2130	x = PY_NEW(Integer)
2131
2132	_sig_on
2133	mpz_sqrt(x.value, self.value)
2134	_sig_off
2135
2136	return x
2137
2138	def sqrt_approx(self, prec=None, all=False):
2139	"""
2140	EXAMPLES:
2141	sage: 5.sqrt_approx(prec=200)
2142	2.2360679774997896964091736687312762354406183596115257242709
2143	sage: 5.sqrt_approx()
2144	2.23606797749979
2145	sage: 4.sqrt_approx()
2146	2
2147	"""
2148	try:
2149	return self.sqrt(extend=False,all=all)
2150	except ValueError:
2151	pass
2152	if prec is None:
2153	prec = max(53, 2*(mpz_sizeinbase(self.value, 2)+2))
2154	return self.sqrt(prec=prec, all=all)
2155
2156	def sqrt(self, prec=None, extend=True, all=False):
2157	"""
2158	The square root function.
2159
2160	INPUT:
2161	prec -- integer (default: None): if None, returns an exact
2162	square root; otherwise returns a numerical square
2163	root if necessary, to the given bits of precision.
2164	extend -- bool (default: True); if True, return a square
2165	root in an extension ring, if necessary. Otherwise,
2166	raise a ValueError if the square is not in the base
2167	ring.
2168	all -- bool (default: False); if True, return all square
2169	roots of self, instead of just one.
2170
2171	EXAMPLES:
2172	sage: Integer(144).sqrt()
2173	12
2174	sage: Integer(102).sqrt()
2175	sqrt(102)
2176
2177	sage: n = 2
2178	sage: n.sqrt(all=True)
2179	[sqrt(2), -sqrt(2)]
2180	sage: n.sqrt(prec=10)
2181	1.4
2182	sage: n.sqrt(prec=100)
2183	1.4142135623730950488016887242
2184	sage: n.sqrt(prec=100,all=True)
2185	[1.4142135623730950488016887242, -1.4142135623730950488016887242]
2186	sage: n.sqrt(extend=False)
2187	Traceback (most recent call last):
2188	...
2189	ValueError: square root of 2 not an integer
2190	sage: Integer(144).sqrt(all=True)
2191	[12, -12]
2192	sage: Integer(0).sqrt(all=True)
2193	[0]
2194	"""
2195	if mpz_sgn(self.value) == 0:
2196	return [self] if all else self
2197
2198	if mpz_sgn(self.value) < 0:
2199	if not extend:
2200	raise ValueError, "square root of negative number not an integer"
2201	if prec:
2202	from sage.rings.complex_field import ComplexField
2203	K = ComplexField(prec)
2204	return K(self).sqrt(all=all)
2205	from sage.calculus.calculus import sqrt
2206	return sqrt(self, all=all)
2207
2208
2209	cdef int non_square
2210	cdef Integer z = PY_NEW(Integer)
2211	cdef mpz_t tmp
2212	_sig_on
2213	mpz_init(tmp)
2214	mpz_sqrtrem(z.value, tmp, self.value)
2215	non_square = mpz_sgn(tmp) != 0
2216	mpz_clear(tmp)
2217	_sig_off
2218
2219	if non_square:
2220	if not extend:
2221	raise ValueError, "square root of %s not an integer"%self
2222	if prec:
2223	from sage.rings.real_mpfr import RealField
2224	K = RealField(prec)
2225	return K(self).sqrt(all=all)
2226	from sage.calculus.calculus import sqrt
2227	return sqrt(self, all=all)
2228
2229	if all:
2230	return [z, -z]
2231	return z
2232
2233	def _xgcd(self, Integer n):
2234	r"""
2235	Return a triple $g, s, t \in\Z$ such that
2236	$$
2237	g = s \cdot \mbox{\rm self} + t \cdot n.
2238	$$
2239
2240	EXAMPLES:
2241	sage: n = 6
2242	sage: g, s, t = n._xgcd(8)
2243	sage: s6 + 8t
2244	2
2245	sage: g
2246	2
2247	"""
2248	cdef mpz_t g, s, t
2249	cdef object g0, s0, t0
2250
2251	mpz_init(g)
2252	mpz_init(s)
2253	mpz_init(t)
2254
2255	_sig_on
2256	mpz_gcdext(g, s, t, self.value, n.value)
2257	_sig_off
2258
2259	g0 = PY_NEW(Integer)
2260	s0 = PY_NEW(Integer)
2261	t0 = PY_NEW(Integer)
2262	set_mpz(g0,g)
2263	set_mpz(s0,s)
2264	set_mpz(t0,t)
2265	mpz_clear(g)
2266	mpz_clear(s)
2267	mpz_clear(t)
2268	return g0, s0, t0
2269
2270	cdef _lshift(self, long int n):
2271	"""
2272	Shift self n bits to the left, i.e., quickly multiply by $2^n$.
2273	"""
2274	cdef Integer x
2275	x = <Integer> PY_NEW(Integer)
2276
2277	_sig_on
2278	if n < 0:
2279	mpz_fdiv_q_2exp(x.value, self.value, -n)
2280	else:
2281	mpz_mul_2exp(x.value, self.value, n)
2282	_sig_off
2283	return x
2284
2285	def __lshift__(x,y):
2286	"""
2287	Shift x y bits to the left.
2288
2289	EXAMPLES:
2290	sage: 32 << 2
2291	128
2292	sage: 32 << int(2)
2293	128
2294	sage: int(32) << 2
2295	128
2296	sage: 1 >> 2.5
2297	Traceback (most recent call last):
2298	...
2299	TypeError: unsupported operands for >>
2300	"""
2301	try:
2302	if not PY_TYPE_CHECK(x, Integer):
2303	x = Integer(x)
2304	elif not PY_TYPE_CHECK(y, Integer):
2305	y = Integer(y)
2306	return (<Integer>x)._lshift(long(y))
2307	except TypeError:
2308	raise TypeError, "unsupported operands for <<"
2309	if PY_TYPE_CHECK(x, Integer) and isinstance(y, (Integer, int, long)):
2310	return (<Integer>x)._lshift(long(y))
2311	return bin_op(x, y, operator.lshift)
2312
2313	cdef _rshift(Integer self, long int n):
2314	cdef Integer x
2315	x = <Integer> PY_NEW(Integer)
2316
2317	_sig_on
2318	if n < 0:
2319	mpz_mul_2exp(x.value, self.value, -n)
2320	else:
2321	mpz_fdiv_q_2exp(x.value, self.value, n)
2322	_sig_off
2323	return x
2324
2325	def __rshift__(x, y):
2326	"""
2327	EXAMPLES:
2328	sage: 32 >> 2
2329	8
2330	sage: 32 >> int(2)
2331	8
2332	sage: int(32) >> 2
2333	8
2334	sage: 1<< 2.5
2335	Traceback (most recent call last):
2336	...
2337	TypeError: unsupported operands for <<
2338	"""
2339	try:
2340	if not PY_TYPE_CHECK(x, Integer):
2341	x = Integer(x)
2342	elif not PY_TYPE_CHECK(y, Integer):
2343	y = Integer(y)
2344	return (<Integer>x)._rshift(long(y))
2345	except TypeError:
2346	raise TypeError, "unsupported operands for >>"
2347
2348	#if PY_TYPE_CHECK(x, Integer) and isinstance(y, (Integer, int, long)):
2349	# return (<Integer>x)._rshift(long(y))
2350	#return bin_op(x, y, operator.rshift)
2351
2352	cdef _and(Integer self, Integer other):
2353	cdef Integer x
2354	x = PY_NEW(Integer)
2355	mpz_and(x.value, self.value, other.value)
2356	return x
2357
2358	def __and__(x, y):
2359	if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
2360	return (<Integer>x)._and(y)
2361	return bin_op(x, y, operator.and_)
2362
2363
2364	cdef _or(Integer self, Integer other):
2365	cdef Integer x
2366	x = PY_NEW(Integer)
2367	mpz_ior(x.value, self.value, other.value)
2368	return x
2369
2370	def __or__(x, y):
2371	"""
2372	Return the bitwise or of the integers x and y.
2373
2374	EXAMPLES:
2375	sage: n = 8; m = 4
2376	sage: n.__or__(m)
2377	12
2378	"""
2379	if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
2380	return (<Integer>x)._or(y)
2381	return bin_op(x, y, operator.or_)
2382
2383
2384	def __invert__(self):
2385	"""
2386	Return the multiplicative interse of self, as a rational number.
2387
2388	EXAMPLE:
2389	sage: n = 10
2390	sage: 1/n
2391	1/10
2392	sage: n.__invert__()
2393	1/10
2394	"""
2395	return Integer(1)/self # todo: optimize
2396
2397	def inverse_of_unit(self):
2398	if mpz_cmp_si(self.value, 1) or mpz_cmp_si(self.value, -1):
2399	return self
2400	else:
2401	raise ZeroDivisionError, "Inverse does not exist."
2402
2403	def inverse_mod(self, n):
2404	"""
2405	Returns the inverse of self modulo $n$, if this inverse exists.
2406	Otherwise, raises a \exception{ZeroDivisionError} exception.
2407
2408	INPUT:
2409	self -- Integer
2410	n -- Integer
2411	OUTPUT:
2412	x -- Integer such that x*self = 1 (mod m), or
2413	raises ZeroDivisionError.
2414	IMPLEMENTATION:
2415	Call the mpz_invert GMP library function.
2416
2417	EXAMPLES:
2418	sage: a = Integer(189)
2419	sage: a.inverse_mod(10000)
2420	4709
2421	sage: a.inverse_mod(-10000)
2422	4709
2423	sage: a.inverse_mod(1890)
2424	Traceback (most recent call last):
2425	...
2426	ZeroDivisionError: Inverse does not exist.
2427	sage: a = Integer(19)**100000
2428	sage: b = a*a
2429	sage: c = a.inverse_mod(b)
2430	Traceback (most recent call last):
2431	...
2432	ZeroDivisionError: Inverse does not exist.
2433	"""
2434	cdef mpz_t x
2435	cdef object ans
2436	cdef int r
2437	cdef Integer m
2438	m = Integer(n)
2439
2440	if m == 1:
2441	return Integer(0)
2442
2443	mpz_init(x)
2444
2445	_sig_on
2446	r = mpz_invert(x, self.value, m.value)
2447	_sig_off
2448
2449	if r == 0:
2450	raise ZeroDivisionError, "Inverse does not exist."
2451	ans = PY_NEW(Integer)
2452	set_mpz(ans,x)
2453	mpz_clear(x)
2454	return ans
2455
2456	def gcd(self, n):
2457	"""
2458	Return the greatest common divisor of self and $n$.
2459
2460	EXAMPLE:
2461	sage: gcd(-1,1)
2462	1
2463	sage: gcd(0,1)
2464	1
2465	sage: gcd(0,0)
2466	0
2467	sage: gcd(2,2^6)
2468	2
2469	sage: gcd(21,2^6)
2470	1
2471	"""
2472	cdef mpz_t g
2473	cdef object g0
2474	cdef Integer _n = Integer(n)
2475
2476	mpz_init(g)
2477
2478
2479	_sig_on
2480	mpz_gcd(g, self.value, _n.value)
2481	_sig_off
2482
2483	g0 = PY_NEW(Integer)
2484	set_mpz(g0,g)
2485	mpz_clear(g)
2486	return g0
2487
2488	def crt(self, y, m, n):
2489	"""
2490	Return the unique integer between $0$ and $mn$ that is
2491	congruent to the integer modulo $m$ and to $y$ modulo $n$. We
2492	assume that~$m$ and~$n$ are coprime.
2493	"""
2494	cdef object g, s, t
2495	cdef Integer _y, _m, _n
2496	_y = Integer(y); _m = Integer(m); _n = Integer(n)
2497	g, s, t = _m.xgcd(_n)
2498	if not g.is_one():
2499	raise ArithmeticError, "CRT requires that gcd of moduli is 1."
2500	# Now sm + tn = 1, so the answer is x + (y-x)sm, where x=self.
2501	return (self + (_y-self)s_m) % (_m*_n)
2502
2503	def test_bit(self, index):
2504	r"""
2505	Return the bit at \code{index}.
2506
2507	EXAMPLES:
2508	sage: w = 6
2509	sage: w.str(2)
2510	'110'
2511	sage: w.test_bit(2)
2512	1
2513	sage: w.test_bit(-1)
2514	0
2515	"""
2516	cdef unsigned long int i
2517	i = index
2518	cdef Integer x
2519	x = Integer(self)
2520	return mpz_tstbit(x.value, i)
2521
2522
2523	ONE = Integer(1)
2524
2525
2526	def LCM_list(v):
2527	cdef int i, n
2528	cdef mpz_t z
2529	cdef Integer w
2530
2531	n = len(v)
2532
2533	if n == 0:
2534	return Integer(1)
2535
2536	try:
2537	w = v[0]
2538	mpz_init_set(z, w.value)
2539
2540	_sig_on
2541	for i from 1 <= i < n:
2542	w = v[i]
2543	mpz_lcm(z, z, w.value)
2544	_sig_off
2545	except TypeError:
2546	w = Integer(v[0])
2547	mpz_init_set(z, w.value)
2548
2549	_sig_on
2550	for i from 1 <= i < n:
2551	w = Integer(v[i])
2552	mpz_lcm(z, z, w.value)
2553	_sig_off
2554
2555
2556	w = PY_NEW(Integer)
2557	mpz_set(w.value, z)
2558	mpz_clear(z)
2559	return w
2560
2561
2562
2563	def GCD_list(v):
2564	cdef int i, n
2565	cdef mpz_t z
2566	cdef Integer w
2567
2568	n = len(v)
2569
2570	if n == 0:
2571	return Integer(1)
2572
2573	try:
2574	w = v[0]
2575	mpz_init_set(z, w.value)
2576
2577	_sig_on
2578	for i from 1 <= i < n:
2579	w = v[i]
2580	mpz_gcd(z, z, w.value)
2581	if mpz_cmp_si(z, 1) == 0:
2582	_sig_off
2583	return Integer(1)
2584	_sig_off
2585	except TypeError:
2586	w = Integer(v[0])
2587	mpz_init_set(z, w.value)
2588
2589	_sig_on
2590	for i from 1 <= i < n:
2591	w = Integer(v[i])
2592	mpz_gcd(z, z, w.value)
2593	if mpz_cmp_si(z, 1) == 0:
2594	_sig_off
2595	return Integer(1)
2596	_sig_off
2597
2598
2599	w = PY_NEW(Integer)
2600	mpz_set(w.value, z)
2601	mpz_clear(z)
2602	return w
2603
2604	def make_integer(s):
2605	cdef Integer r
2606	r = PY_NEW(Integer)
2607	r._reduce_set(s)
2608	return r
2609
2610	# TODO: where is this used? Never doctested. Should I call IntegerRing.random_element()?
2611	from random import randint
2612	def random_integer(min=-2, max=2):
2613	cdef Integer x
2614	cdef int _min, _max, r
2615	try:
2616	_min = min
2617	_max = max
2618	x = PY_NEW(Integer)
2619	r = randint() % (_max - _min + 1) + _min
2620	mpz_set_si(x.value, r)
2621	return x
2622	except OverflowError:
2623	return Integer(randint(min,max))
2624
2625
2626	############### INTEGER CREATION CODE #####################
2627	######## There is nothing to see here, move along #######
2628
2629	cdef extern from *:
2630
2631	ctypedef struct RichPyObject "PyObject"
2632
2633	# We need a PyTypeObject with elements so we can
2634	# get and set tp_new, tp_dealloc, tp_flags, and tp_basicsize
2635	ctypedef struct RichPyTypeObject "PyTypeObject":
2636
2637	# We replace this one
2638	PyObject* (* tp_new) ( RichPyTypeObject, PyObject, PyObject*)
2639
2640	# Not used, may be useful to determine correct memory management function
2641	RichPyObject ( tp_alloc) ( RichPyTypeObject*, size_t )
2642
2643	# We replace this one
2644	void (tp_dealloc) ( PyObject)
2645
2646	# Not used, may be useful to determine correct memory management function
2647	void (* tp_free) ( PyObject* )
2648
2649	# sizeof(Object)
2650	size_t tp_basicsize
2651
2652	# We set a flag here to circumvent the memory manager
2653	long tp_flags
2654
2655	cdef long Py_TPFLAGS_HAVE_GC
2656
2657	# We need a PyObject where we can get/set the refcnt directly
2658	# and access the type.
2659	ctypedef struct RichPyObject "PyObject":
2660	int ob_refcnt
2661	RichPyTypeObject* ob_type
2662
2663	# Allocation
2664	RichPyObject* PyObject_MALLOC(int)
2665
2666	# Useful for debugging, see below
2667	void PyObject_INIT(RichPyObject , RichPyTypeObject )
2668
2669	# Free
2670	void PyObject_FREE(PyObject*)
2671
2672
2673	# We need a couple of internal GMP datatypes.
2674
2675	# This may be potentialy very dangerous as it reaches
2676	# deeply into the internal structure of GMP which may not
2677	# be consistant across future versions of GMP.
2678	# See extensive note in the fast_tp_new() function below.
2679
2680	cdef extern from "gmp.h":
2681	ctypedef void* mp_ptr #"mp_ptr"
2682
2683	# We allocate _mp_d directly (mpz_t is typedef of this in GMP)
2684	ctypedef struct __mpz_struct "__mpz_struct":
2685	mp_ptr _mp_d
2686	size_t _mp_alloc
2687	size_t _mp_size
2688
2689	# sets the three free, alloc, and realloc function pointers to the
2690	# memory management functions set in GMP. Accepts NULL pointer.
2691	# Potentially dangerous if changed by calling
2692	# mp_set_memory_functions again after we initialized this module.
2693	void mp_get_memory_functions (void (alloc) (size_t), void (*realloc)(void , size_t, size_t), void (*free) (void , size_t))
2694
2695	# GMP's configuration of how many Bits are stuffed into a limb
2696	cdef int __GMP_BITS_PER_MP_LIMB
2697
2698	# This variable holds the size of any Integer object in bytes.
2699	cdef int sizeof_Integer
2700
2701	# We use a global Integer element to steal all the references
2702	# from. DO NOT INITIALIZE IT AGAIN and DO NOT REFERENCE IT!
2703	cdef Integer global_dummy_Integer
2704	global_dummy_Integer = Integer()
2705
2706	# Accessing the .value attribute of an Integer object causes Pyrex to
2707	# refcount it. This is problematic, because that causes overhead and
2708	# more importantly an infinite loop in the destructor. If you refcount
2709	# in the destructor and the refcount reaches zero (which is true
2710	# everytime) the destructor is called.
2711	#
2712	# To avoid this we calculate the byte offset of the value member and
2713	# remember it in this variable.
2714	#
2715	# Eventually this may be rendered obsolete by a change in SageX allowing
2716	# non-reference counted extension types.
2717	cdef long mpz_t_offset
2718
2719
2720	# stores the GMP alloc function
2721	cdef void * (* mpz_alloc)(size_t)
2722
2723	# stores the GMP free function
2724	cdef void (* mpz_free)(void *, size_t)
2725
2726	# A global pool for performance when integers are rapidly created and destroyed.
2727	# It operates on the following principles:
2728	#
2729	# - The pool starts out empty.
2730	# - When an new integer is needed, one from the pool is returned
2731	# if available, otherwise a new Integer object is created
2732	# - When an integer is collected, it will add it to the pool
2733	# if there is room, otherwise it will be deallocated.
2734
2735	cdef enum:
2736	integer_pool_size = 100 # Pyrex has no way of defining constants
2737
2738	cdef PyObject* integer_pool[integer_pool_size]
2739	cdef int integer_pool_count = 0
2740
2741	# used for profiling the pool
2742	cdef int total_alloc = 0
2743	cdef int use_pool = 0
2744
2745	# The signature of tp_new is
2746	# PyObject* tp_new(RichPyTypeObject t, PyObject a, PyObject *k).
2747	# However we only use t in this implementation.
2748	#
2749	# t in this case is the Integer TypeObject.
2750
2751	cdef PyObject* fast_tp_new(RichPyTypeObject t, PyObject a, PyObject *k):
2752
2753	global integer_pool, integer_pool_count, total_alloc, use_pool
2754
2755	cdef RichPyObject* new
2756
2757	# for profiling pool usage
2758	# total_alloc += 1
2759
2760	# If there is a ready integer in the pool, we will
2761	# decrement the counter and return that.
2762
2763	if integer_pool_count > 0:
2764
2765	# for profiling pool usage
2766	# use_pool += 1
2767
2768	integer_pool_count -= 1
2769	new = <RichPyObject *> integer_pool[integer_pool_count]
2770
2771	# Otherwise, we have to create one.
2772
2773	else:
2774
2775	# allocate enough room for the Integer, sizeof_Integer is
2776	# sizeof(Integer). The use of PyObject_MALLOC directly
2777	# assumes that Integers are not garbage collected, i.e.
2778	# they do not pocess references to other Python
2779	# objects (Aas indicated by the Py_TPFLAGS_HAVE_GC flag).
2780	# See below for a more detailed description.
2781
2782	new = PyObject_MALLOC( sizeof_Integer )
2783
2784	# Now set every member as set in z, the global dummy Integer
2785	# created before this tp_new started to operate.
2786
2787	memcpy(new, (<void*>global_dummy_Integer), sizeof_Integer )
2788
2789	# This line is only needed if Python is compiled in debugging
2790	# mode './configure --with-pydebug'. If that is the case a Python
2791	# object has a bunch of debugging fields which are initialized
2792	# with this macro. For speed reasons, we don't call it if Python
2793	# is not compiled in debug mode. So uncomment the following line
2794	# if you are debugging Python.
2795
2796	#PyObject_INIT(new, (<RichPyObject*>global_dummy_Integer).ob_type)
2797
2798	# We take the address 'new' and move mpz_t_offset bytes (chars)
2799	# to the address of 'value'. We treat that address as a pointer
2800	# to a mpz_t struct and allocate memory for the _mp_d element of
2801	# that struct. We allocate one limb.
2802	#
2803	# What is done here is potentialy very dangerous as it reaches
2804	# deeply into the internal structure of GMP. Consequently things
2805	# may break if a new release of GMP changes some internals. To
2806	# emphazise this, this is what the GMP manual has to say about
2807	# the documentation for the struct we are using:
2808	#
2809	# "This chapter is provided only for informational purposes and the
2810	# various internals described here may change in future GMP releases.
2811	# Applications expecting to be compatible with future releases should use
2812	# only the documented interfaces described in previous chapters."
2813	#
2814	# If this line is used SAGE is not such an application.
2815	#
2816	# The clean version of the following line is:
2817	#
2818	# mpz_init(( <mpz_t>(<char *>new + mpz_t_offset) )
2819	#
2820	# We save time both by avoiding an extra function call and
2821	# because the rest of the mpz struct was already initalized
2822	# fully using the memcpy above.
2823
2824	(<__mpz_struct >( <char >new + mpz_t_offset) )._mp_d = <mp_ptr>mpz_alloc(__GMP_BITS_PER_MP_LIMB >> 3)
2825
2826	# The global_dummy_Integer may have a reference count larger than
2827	# one, but it is expected that newly created objects have a
2828	# reference count of one. This is potentially unneeded if
2829	# everybody plays nice, because the gobal_dummy_Integer has only
2830	# one reference in that case.
2831
2832	# Objects from the pool have reference count zero, so this
2833	# needs to be set in this case.
2834
2835	new.ob_refcnt = 1
2836
2837	return new
2838
2839	cdef void fast_tp_dealloc(PyObject* o):
2840
2841	# If there is room in the pool for a used integer object,
2842	# then put it in rather than deallocating it.
2843
2844	global integer_pool, integer_pool_count
2845
2846	if integer_pool_count < integer_pool_size:
2847
2848	# Here we free any extra memory used by the mpz_t by
2849	# setting it to a single limb.
2850	if (<__mpz_struct >( <char >o + mpz_t_offset))._mp_alloc > 1:
2851	_mpz_realloc(<mpz_t >( <char >o + mpz_t_offset), 1)
2852
2853	# It's cheap to zero out an integer, so do it here.
2854	(<__mpz_struct >( <char >o + mpz_t_offset))._mp_size = 0
2855
2856	# And add it to the pool.
2857	integer_pool[integer_pool_count] = o
2858	integer_pool_count += 1
2859	return
2860
2861	# Again, we move to the mpz_t and clear it. See above, why this is evil.
2862	# The clean version of this line would be:
2863	# mpz_clear(<mpz_t>(<char *>o + mpz_t_offset))
2864
2865	mpz_free((<__mpz_struct >( <char >o + mpz_t_offset) )._mp_d, 0)
2866
2867	# Free the object. This assumes that Py_TPFLAGS_HAVE_GC is not
2868	# set. If it was set another free function would need to be
2869	# called.
2870
2871	PyObject_FREE(o)
2872
2873	hook_fast_tp_functions()
2874
2875	def hook_fast_tp_functions():
2876	"""
2877	"""
2878	global global_dummy_Integer, mpz_t_offset, sizeof_Integer
2879
2880	cdef long flag
2881
2882	cdef RichPyObject* o
2883	o = <RichPyObject*>global_dummy_Integer
2884
2885	# By default every object created in Pyrex is garbage
2886	# collected. This means it may have references to other objects
2887	# the Garbage collector has to look out for. We remove this flag
2888	# as the only reference an Integer has is to the global Integer
2889	# ring. As this object is unique we don't need to garbage collect
2890	# it as we always have a module level reference to it. If another
2891	# attribute is added to the Integer class this flag removal so as
2892	# the alloc and free functions may not be used anymore.
2893	# This object will still be reference counted.
2894	flag = Py_TPFLAGS_HAVE_GC
2895	o.ob_type.tp_flags = <long>(o.ob_type.tp_flags & (~flag))
2896
2897	# calculate the offset of the GMP mpz_t to avoid casting to/from
2898	# an Integer which includes reference counting. Reference counting
2899	# is bad in constructors and destructors as it potentially calls
2900	# the destructor.
2901	# Eventually this may be rendered obsolete by a change in SageX allowing
2902	# non-reference counted extension types.
2903	mpz_t_offset = <char >(&global_dummy_Integer.value) - <char >o
2904
2905	# store how much memory needs to be allocated for an Integer.
2906	sizeof_Integer = o.ob_type.tp_basicsize
2907
2908	# get the functions to do memory management for the GMP elements
2909	# WARNING: if the memory management functions are changed after
2910	# this initialisation, we are/you are doomed.
2911
2912	mp_get_memory_functions(&mpz_alloc, NULL, &mpz_free)
2913
2914	# Finally replace the functions called when an Integer needs
2915	# to be constructed/destructed.
2916	o.ob_type.tp_new = &fast_tp_new
2917	o.ob_type.tp_dealloc = &fast_tp_dealloc
2918
2919	def time_alloc_list(n):
2920	cdef int i
2921	l = []
2922	for i from 0 <= i < n:
2923	l.append(PY_NEW(Integer))
2924
2925	return l
2926
2927	def time_alloc(n):
2928	cdef int i
2929	for i from 0 <= i < n:
2930	z = PY_NEW(Integer)
2931
2932	def pool_stats():
2933	print "Used pool %s / %s times" % (use_pool, total_alloc)
2934	print "Pool contains %s / %s items" % (integer_pool_count, integer_pool_size)
2935
2936	cdef integer(x):
2937	if PY_TYPE_CHECK(x, Integer):
2938	return x
2939	return Integer(x)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: