Context Navigation

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

Visit:

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

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 "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, 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(173737)
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
2026	def additive_order(self):
2027	"""
2028	Return the additive order of self.
2029
2030	EXAMPLES:
2031	sage: ZZ(0).additive_order()
2032	1
2033	sage: ZZ(1).additive_order()
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: s6 + 8t
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 sm + tn = 1, so the answer is x + (y-x)sm, 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
2525	ONE = Integer(1)
2526
2527
2528	def 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
2565	def 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
2606	def 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()?
2613	from random import randint
2614	def 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
2631	cdef 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
2682	cdef 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.
2701	cdef 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!
2705	cdef Integer global_dummy_Integer
2706	global_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.
2719	cdef long mpz_t_offset
2720
2721
2722	# stores the GMP alloc function
2723	cdef void * (* mpz_alloc)(size_t)
2724
2725	# stores the GMP free function
2726	cdef 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
2737	cdef enum:
2738	integer_pool_size = 100 # Pyrex has no way of defining constants
2739
2740	cdef PyObject* integer_pool[integer_pool_size]
2741	cdef int integer_pool_count = 0
2742
2743	# used for profiling the pool
2744	cdef int total_alloc = 0
2745	cdef 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
2753	cdef 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
2841	cdef 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
2875	hook_fast_tp_functions()
2876
2877	def 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
2921	def 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
2929	def time_alloc(n):
2930	cdef int i
2931	for i from 0 <= i < n:
2932	z = PY_NEW(Integer)
2933
2934	def 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
2938	cdef 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.

Download in other formats: