Context Navigation

source: sage/ext/integer.pyx @ 5:33ee5f0ec62e

Visit:

Revision 5:33ee5f0ec62e, 28.9 KB checked in by Gonzalo Tornaria <tornaria@…>, 7 years ago (diff)
[project @ patch to sage-1.0.2]

Line
1	r"""
2	Elements of the ring $\Z$ of integers
3	"""
4
5	#*****************************************************************************
6	# Copyright (C) 2004 William Stein <wstein@ucsd.edu>
7	#
8	# Distributed under the terms of the GNU General Public License (GPL)
9	#
10	# This code is distributed in the hope that it will be useful,
11	# but WITHOUT ANY WARRANTY; without even the implied warranty of
12	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13	# General Public License for more details.
14	#
15	# The full text of the GPL is available at:
16	#
17	# http://www.gnu.org/licenses/
18	#*****************************************************************************
19
20	doc="""
21	Integers
22	"""
23
24	import operator
25
26	include "gmp.pxi"
27	include "interrupt.pxi" # ctrl-c interrupt block support
28
29	cdef class Integer(element.EuclideanDomainElement)
30
31	import sage.rings.integer_ring
32	import sage.rings.coerce
33	import sage.rings.infinity
34	import sage.rings.complex_field
35	import rational as rational
36	import sage.libs.pari.all
37	import mpfr
38
39	cdef public int set_mpz(Integer self, mpz_t value):
40	mpz_set(self.value, value)
41
42	cdef set_from_Integer(Integer self, Integer other):
43	mpz_set(self.value, other.value)
44
45	cdef set_from_int(Integer self, int other):
46	mpz_set_si(self.value, other)
47
48	cdef public mpz_t* get_value(Integer self):
49	return &self.value
50
51	# This crashes SAGE:
52	# s = 2003^100300000
53	# The problem is related to realloc moving all the memory
54	# and returning a pointer to the new block of memory, I think.
55
56	cdef extern from "stdlib.h":
57	void abort()
58
59	cdef void* pymem_realloc(void *ptr, size_t old_size, size_t new_size):
60	return PyMem_Realloc(ptr, new_size)
61	#print "hello-realloc: ", <long> ptr, old_size, new_size
62	#cdef void* p
63	#p = PyMem_Realloc(ptr, new_size)
64	#print "p = ", <long> p
65	#if <void> p != <void> ptr:
66	# abort()
67	#return p
68
69
70	cdef void pymem_free(void *ptr, size_t size):
71	PyMem_Free(ptr)
72
73	cdef void* pymem_malloc(size_t size):
74	return PyMem_Malloc(size)
75
76	cdef extern from "gmp.h":
77	void mp_set_memory_functions (void (alloc_func_ptr) (size_t), void (realloc_func_ptr) (void , size_t, size_t), void (free_func_ptr) (void *, size_t))
78
79
80	def pmem_malloc():
81	"""
82	Use the Python heap for GMP memory management.
83	"""
84	mp_set_memory_functions(PyMem_Malloc, pymem_realloc, pymem_free)
85	#mp_set_memory_functions(pymem_malloc, pymem_realloc, pymem_free)
86
87	pmem_malloc()
88
89	cdef class Integer(element.EuclideanDomainElement):
90	"""
91	The \\class{Integer} class represents arbitrary precision
92	integers. It derives from the \\class{Element} class, so
93	integers can be used as ring elements anywhere in SAGE.
94
95	\\begin{notice}
96	The class \\class{Integer} is implemented in Pyrex, as a wrapper
97	of the GMP \\code{mpz_t} integer type.
98	\\end{notice}
99	"""
100	def __new__(self, x=None):
101	mpz_init(self.value)
102
103	def __pyxdoc__init__(self):
104	"""
105	You can create an integer from an int, long, string literal, or
106	integer modulo N.
107
108	EXAMPLES:
109	sage: Integer(495)
110	495
111	sage: Integer('495949209809328523')
112	495949209809328523
113	sage: Integer(Mod(3,7))
114	3
115
116	Integers also support the standard arithmetic operations, such
117	as +,-,*,/,^, \code{abs}, \code{mod}, \code{float}:
118	sage: 2^3
119	8
120	"""
121	def __init__(self, x=None):
122	if not (x is None):
123	if isinstance(x, Integer):
124	set_from_Integer(self, x)
125
126	elif hasattr(x, "_integer_"):
127	set_from_Integer(self, x._integer_())
128
129	elif isinstance(x, int):
130	mpz_set_si(self.value, x)
131
132	elif isinstance(x, long):
133	s = "%x"%x
134	mpz_set_str(self.value, s, 16)
135
136	elif isinstance(x, str):
137	if mpz_set_str(self.value, x, 0) != 0:
138	raise TypeError, "unable to convert x (=%s) to an integer"%x
139
140	elif isinstance(x, rational.Rational):
141	if x.denominator() != 1:
142	raise TypeError, "Unable to coerce rational (=%s) to an Integer."%x
143	set_from_Integer(self, x.numer())
144
145
146	elif isinstance(x, sage.libs.pari.all.gen):
147	if x.type() == 't_INTMOD':
148	x = x.lift()
149	# TODO: figure out how to convert to pari integer in base 16 ?
150	s = str(x)
151	if mpz_set_str(self.value, s, 0) != 0:
152	raise TypeError, "Unable to coerce PARI %s to an Integer."%x
153
154	else:
155	raise TypeError, "Unable to coerce %s (of type %s) to an Integer."%(x,type(x))
156
157
158	def __reduce__(self):
159	# This single line below took me HOURS to figure out.
160	# It is the trick needed to pickle pyrex extension types.
161	# The trick is that you must put a pure Python function
162	# as the first argument, and that function must return
163	# the result of unpickling with the argument in the second
164	# tuple as input. All kinds of problems happen
165	# if we don't do this.
166	return sage.rings.integer.make_integer, (self.str(32),)
167
168	def _reduce_set(self, s):
169	mpz_set_str(self.value, s, 32)
170
171	def _im_gens_(self, codomain, im_gens):
172	return codomain._coerce_(self)
173
174	#def __reduce__(self):
175	# return sage.rings.integer_ring.Z, (long(self),)
176
177	cdef int cmp(self, Integer x):
178	cdef int i
179	i = mpz_cmp(self.value, x.value)
180	if i < 0:
181	return -1
182	elif i == 0:
183	return 0
184	else:
185	return 1
186
187	def __richcmp__(Integer self, right, int op):
188	cdef int n
189	if not isinstance(right, Integer):
190	try:
191	n = sage.rings.coerce.cmp(self, right)
192	except TypeError:
193	n = -1
194	else:
195	n = self.cmp(right)
196	return self._rich_to_bool(op, n)
197
198	def __cmp__(self, x):
199	return self.cmp(x)
200
201
202	def copy(self):
203	"""
204	Return a copy of the integer.
205	"""
206	cdef Integer z
207	z = Integer()
208	set_mpz(z,self.value)
209	return z
210
211	def __dealloc__(self):
212	mpz_clear(self.value)
213
214	def __repr__(self):
215	return self.str()
216
217	def _latex_(self):
218	return self.str()
219
220	def __str_malloc(self, int base=10):
221	"""
222	Return the string representation of \\code{self} in the given
223	base. (Uses malloc then PyMem. This is actually slightly
224	faster than self.str() below, but it is unpythonic to use
225	malloc.) However, self.str() below is nice because we know
226	the size of the string ahead of time, and can work around a
227	bug in GMP nicely. There seems to be a bug in GMP, where
228	non-2-power base conversion for very large integers > 10
229	million digits (?) crashes GMP.
230	"""
231	_sig_on
232	cdef char *s
233	s = mpz_get_str(NULL, base, self.value)
234	t = str(s)
235	free(s)
236	_sig_off
237	return t
238
239	def str(self, int base=10):
240	r"""
241	Return the string representation of \code{self} in the given
242	base.
243
244	\begin{notice}
245	String representation of integers with more than about {\em four
246	million digits} is not available in bases other than a power of
247	2. This is because of a bug in GMP. To obtain the base-$b$
248	expansion of an integer $n$ use \code{n.str(b)}.
249	\end{notice}
250	"""
251	if base < 2 or base > 36:
252	raise ValueError, "base (=%s) must be between 2 and 36"%base
253	cdef size_t n
254	cdef char *s
255	n = mpz_sizeinbase(self.value, base) + 2
256	if n > 4200000 and base != 2 and base != 4 and base != 8 and base != 16 and base != 32:
257	raise RuntimeError, "String representation of integers with more than 4200000 digits is not available in bases other than a power of 2. This is because of a bug in GMP. To obtain the base-b expansion of an integer n use n.str(b)."
258	s = <char *>PyMem_Malloc(n)
259	if s == NULL:
260	raise MemoryError, "Unable to allocate enough memory for the string representation of an integer."
261	_sig_on
262	mpz_get_str(s, base, self.value)
263	_sig_off
264	k = PyString_FromString(s)
265	PyMem_Free(s)
266	return k
267
268	def __hex__(self):
269	r"""
270	Return the hexadecimal digits of self in lower case.
271
272	\note{'0x' is \emph{not} prepended to the result like is done
273	by the corresponding Python function on int or long. This is
274	for efficiency sake---adding and stripping the string wastes
275	time; since this function is used for conversions from
276	integers to other C-library structures, it is important that
277	it be fast.}
278
279	EXAMPLES:
280	sage: print hex(Integer(15))
281	f
282	sage: print hex(Integer(16))
283	10
284	sage: print hex(Integer(16938402384092843092843098243))
285	36bb1e3929d1a8fe2802f083
286	sage: print hex(long(16938402384092843092843098243))
287	0x36BB1E3929D1A8FE2802F083L
288	"""
289	return self.str(16)
290
291	def binary(self):
292	"""
293	Return the binary digits of self as a string.
294
295	EXAMPLES:
296	sage: print Integer(15).binary()
297	1111
298	sage: print Integer(16).binary()
299	10000
300	sage: print Integer(16938402384092843092843098243).binary()
301	1101101011101100011110001110010010100111010001101010001111111000101000000000101111000010000011
302	"""
303	return self.str(2)
304
305
306
307	def set_si(self, signed long int n):
308	"""
309	Coerces $n$ to a C signed integer if possible, and sets self
310	equal to $n$.
311
312	"""
313	mpz_set_si(self.value, n)
314
315	def set_str(self, s, base=10):
316	"""
317	Set self equal to the number defined by the string s in the
318	given base.
319	"""
320	valid = mpz_set_str(self.value, s, base)
321	if valid != 0:
322	raise TypeError, "unable to convert x (=%s) to an integer"%s
323
324	cdef void set_from_mpz(Integer self, mpz_t value):
325	mpz_set(self.value, value)
326
327	cdef mpz_t* get_value(Integer self):
328	return &self.value
329
330	def __add_(Integer self, Integer other):
331	cdef Integer x
332	x = Integer()
333	mpz_add(x.value, self.value, other.value)
334	return x
335
336	def __add__(x, y):
337	if isinstance(x, Integer) and isinstance(y, Integer):
338	return x.__add_(y)
339	return sage.rings.coerce.bin_op(x, y, operator.add)
340
341	def __sub_(Integer self, Integer other):
342	cdef Integer x
343	x = Integer()
344	mpz_sub(x.value, self.value, other.value)
345	return x
346
347	def __sub__(x, y):
348	"""
349	EXAMPLES:
350	sage: a = Integer(5); b = Integer(3)
351	sage: b - a
352	-2
353	"""
354	if isinstance(x, Integer) and isinstance(y, Integer):
355	return x.__sub_(y)
356	return sage.rings.coerce.bin_op(x, y, operator.sub)
357
358	def __mul_(Integer self, Integer other):
359	cdef Integer x
360	x = Integer()
361
362	_sig_on
363	mpz_mul(x.value, self.value, other.value)
364	_sig_off
365
366	return x
367
368	def __mul__(x, y):
369	if isinstance(x, Integer) and isinstance(y, Integer):
370	return x.__mul_(y)
371	if isinstance(x, list):
372	return x * int(y)
373	return sage.rings.coerce.bin_op(x, y, operator.mul)
374
375	def __div_(Integer self, Integer other):
376	cdef Integer x
377	#if mpz_divisible_p(self.value, other.value):
378	# x = Integer()
379	# mpz_divexact(x.value, self.value, other.value)
380	# return x
381	#else:
382	return rational.Rational(self)/rational.Rational(other)
383	#raise ArithmeticError, "Exact division impossible."
384
385	def __div__(x, y):
386	if isinstance(x, Integer) and isinstance(y, Integer):
387	return x.__div_(y)
388	return sage.rings.coerce.bin_op(x, y, operator.div)
389
390	def __floordiv(Integer self, Integer other):
391	cdef Integer x
392	x = Integer()
393
394
395	_sig_on
396	mpz_fdiv_q(x.value, self.value, other.value)
397	_sig_off
398
399
400	return x
401
402
403	def __floordiv__(x, y):
404	if isinstance(x, Integer) and isinstance(y, Integer):
405	return x.__floordiv(y)
406	return sage.rings.coerce.bin_op(x, y, operator.floordiv)
407
408
409	def __pow__(self, n, dummy):
410	cdef Integer _self, _n
411	cdef unsigned int _nval
412	if not isinstance(self, Integer):
413	return self.__pow__(int(n))
414	_n = Integer(n)
415	if _n < 0:
416	return Integer(1)/(self**(-_n))
417	_self = integer(self)
418	cdef Integer x
419	x = Integer()
420	_nval = _n
421
422	_sig_on
423	mpz_pow_ui(x.value, _self.value, _nval)
424	_sig_off
425
426	return x
427
428	def __pos__(self):
429	return self
430
431	def __neg__(self):
432	cdef Integer x
433	x = Integer()
434	mpz_neg(x.value, self.value)
435	return x
436
437	def __abs__(self):
438	"""
439	"""
440	cdef Integer x
441	x = Integer()
442	mpz_abs(x.value, self.value)
443	return x
444
445	def __mod__(self, modulus):
446	cdef Integer _modulus, _self
447	_modulus = integer(modulus)
448	if not _modulus:
449	raise ZeroDivisionError, "Integer modulo by zero"
450	_self = integer(self)
451
452	cdef Integer x
453	x = Integer()
454
455	_sig_on
456	mpz_mod(x.value, _self.value, _modulus.value)
457	_sig_off
458
459	return x
460
461
462	def quo_rem(self, other):
463	cdef Integer _other, _self
464	_other = integer(other)
465	if not _other:
466	raise ZeroDivisionError, "other (=%s) must be nonzero"%other
467	_self = integer(self)
468
469	cdef Integer q, r
470	q = Integer()
471	r = Integer()
472
473	_sig_on
474	mpz_tdiv_qr(q.value, r.value, _self.value, _other.value)
475	_sig_off
476
477	return q, r
478
479
480
481	def powermod(self, exp, mod):
482	"""
483	powermod(self, Integer exp, Integer mod):
484	Compute self**exp modulo mod.
485	"""
486	cdef Integer x, _exp, _mod
487	_exp = Integer(exp); _mod = Integer(mod)
488	x = Integer()
489
490	_sig_on
491	mpz_powm(x.value, self.value, _exp.value, _mod.value)
492	_sig_off
493
494	return x
495
496	def powermodm_ui(self, unsigned long int exp, mod):
497	cdef Integer x, _mod
498	_mod = Integer(mod)
499	x = Integer()
500
501	_sig_on
502	mpz_powm_ui(x.value, self.value, exp, _mod.value)
503	_sig_off
504
505	return x
506
507	def __int__(self):
508	cdef char *s
509	s = mpz_get_str(NULL, 16, self.value)
510	n = int(s,16)
511	free(s)
512	return n
513
514	def __nonzero__(self):
515	return not self.is_zero()
516
517	def __long__(self):
518	cdef char *s
519	s = mpz_get_str(NULL, 16, self.value)
520	n = long(s,16)
521	free(s)
522	return n
523
524	def __float__(self):
525	return mpz_get_d(self.value)
526
527	def __hash__(self):
528	return hash(self.str(16))
529
530	def factor(self, algorithm='pari'):
531	"""
532	Return the prime factorization of the integer as a list of
533	pairs $(p,e)$, where $p$ is prime and $e$ is a positive integer.
534
535	INPUT:
536	algorithm -- string
537	* 'pari' -- (default) use the PARI c library
538	* 'kash' -- use KASH computer algebra system (requires
539	the optional kash package be installed)
540	"""
541	return sage.rings.integer_ring.factor(self, algorithm=algorithm)
542
543	def coprime_integers(self, m):
544	"""
545	Return the positive integers $< m$ that are coprime to self.
546
547	EXAMPLES:
548	sage: 8.coprime_integers(8)
549	[1, 3, 5, 7]
550	sage: 8.coprime_integers(11)
551	[1, 3, 5, 7, 9]
552	sage: 5.coprime_integers(10)
553	[1, 2, 3, 4, 6, 7, 8, 9]
554	sage: 5.coprime_integers(5)
555	[1, 2, 3, 4]
556	sage: 99.coprime_integers(99)
557	[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]
558
559	AUTHORS:
560	-- Naqi Jaffery (2006-01-24): examples
561	"""
562	# TODO -- make VASTLY faster
563	v = []
564	for n in range(1,m):
565	if self.gcd(n) == 1:
566	v.append(Integer(n))
567	return v
568
569	def divides(self, n):
570	"""
571	Return True if self divides n.
572
573	EXAMPLES:
574	sage: Z = IntegerRing()
575	sage: Z(5).divides(Z(10))
576	True
577	sage: Z(0).divides(Z(5))
578	False
579	sage: Z(10).divides(Z(5))
580	False
581	"""
582	cdef int t
583	cdef Integer _n
584	_n = Integer(n)
585	if mpz_cmp_si(self.value, 0) == 0:
586	return bool(mpz_cmp_si(_n.value, 0) == 0)
587	_sig_on
588	t = mpz_divisible_p(_n.value, self.value)
589	_sig_off
590	return bool(t)
591
592
593	def valuation(self, p):
594	if self == 0:
595	return sage.rings.infinity.infinity
596	cdef int k
597	k = 0
598	while self % p == 0:
599	k = k + 1
600	self = self.__floordiv__(p)
601	return Integer(k)
602
603	def _lcm(self, Integer n):
604	"""
605	Returns the least common multiple of self and $n$.
606	"""
607	cdef mpz_t x
608
609	mpz_init(x)
610
611	_sig_on
612	mpz_lcm(x, self.value, n.value)
613	_sig_off
614
615
616	cdef Integer z
617	z = Integer()
618	mpz_set(z.value,x)
619	mpz_clear(x)
620	return z
621
622	def denominator(self):
623	"""
624	Return the denominator of this integer.
625
626	EXAMPLES:
627	sage: x = 5
628	sage: x.denominator()
629	1
630	sage: x = 0
631	sage: x.denominator()
632	1
633	"""
634	return ONE
635
636	def numerator(self):
637	"""
638	Return the numerator of this integer.
639
640	EXAMPLE:
641	sage: x = 5
642	sage: x.numerator()
643	5
644
645	sage: x = 0
646	sage: x.numerator()
647	0
648	"""
649	return self
650
651	def is_one(self):
652	"""
653	Returns \\code{True} if the integers is $1$, otherwise \\code{False}.
654	"""
655	return mpz_cmp_si(self.value, 1) == 0
656
657	def is_zero(self):
658	"""
659	Returns \\code{True} if the integers is $0$, otherwise \\code{False}.
660	"""
661	return mpz_cmp_si(self.value, 0) == 0
662
663	def is_unit(self):
664	"""
665	Returns \\code{true} if this integer is a unit, i.e., 1 or $-1$.
666	"""
667	return bool(mpz_cmp_si(self.value, -1) == 0 or mpz_cmp_si(self.value, 1) == 0)
668
669	def is_square(self):
670	return bool(self._pari_().issquare())
671
672	def is_prime(self):
673	return bool(self._pari_().isprime())
674
675	def square_free_part(self):
676	"""
677	Return the square free part of $x$, i.e., a divisor z such that $x = z y^2$,
678	for a perfect square $y^2$.
679
680	EXAMPLES:
681	sage: square_free_part(100)
682	1
683	sage: square_free_part(12)
684	3
685	sage: square_free_part(173737)
686	17
687	sage: square_free_part(-17*32)
688	-34
689	sage: square_free_part(1)
690	1
691	sage: square_free_part(-1)
692	-1
693	sage: square_free_part(-2)
694	-2
695	sage: square_free_part(-4)
696	-1
697	"""
698	if self.is_zero():
699	return self
700	F = self.factor()
701	n = Integer(1)
702	for p, e in F:
703	if e % 2 != 0:
704	n = n * p
705	return n * F.unit()
706
707	def next_prime(self):
708	return Integer( (self._pari_()+1).nextprime())
709
710	def additive_order(self):
711	"""
712	Return the additive order of self.
713
714	EXAMPLES:
715	sage: ZZ(0).additive_order()
716	1
717	sage: ZZ(1).additive_order()
718	Infinity
719	"""
720	import sage.rings.infinity
721	if self.is_zero():
722	return Integer(1)
723	else:
724	return sage.rings.infinity.infinity
725
726	def multiplicative_order(self):
727	r"""
728	Return the multiplicative order of self, if self is a unit, or raise
729	\code{ArithmeticError} otherwise.
730
731	EXAMPLES:
732	sage: ZZ(1).multiplicative_order()
733	1
734	sage: ZZ(-1).multiplicative_order()
735	2
736	sage: ZZ(0).multiplicative_order()
737	Traceback (most recent call last):
738	...
739	ArithmeticError: no power of 0 is a unit
740	sage: ZZ(2).multiplicative_order()
741	Traceback (most recent call last):
742	...
743	ArithmeticError: no power of 2 is a unit
744	"""
745	if mpz_cmp_si(self.value, 1) == 0:
746	return Integer(1)
747	elif mpz_cmp_si(self.value, -1) == 0:
748	return Integer(2)
749	else:
750	raise ArithmeticError, "no power of %s is a unit"%self
751
752	def is_square_free(self):
753	"""
754	Returns True if this integer is not divisible by the square of
755	any prime and False otherwise.
756	"""
757	return self._pari_().issquarefree()
758
759	def _pari_(self):
760	if self._pari is None:
761	# better to do in hex, but I can't figure out
762	# how to input/output a number in hex in PARI!!
763	# TODO: (I could just think carefully about raw bytes and make this all much faster...)
764	self._pari = sage.libs.pari.all.pari(str(self))
765	return self._pari
766
767	def parent(self):
768	"""
769	Return the ring $\\Z$ of integers.
770	"""
771	return sage.rings.integer_ring.Z
772
773	def isqrt(self):
774	"""
775	Returns the integer floor of the square root of self, or raises
776	an \\exception{ValueError} if self is negative.
777
778	EXAMPLE:
779	sage: a = Integer(5)
780	sage: a.isqrt()
781	2
782	"""
783	if self < 0:
784	raise ValueError, "square root of negative number not defined."
785	cdef Integer x
786	x = Integer()
787
788	_sig_on
789	mpz_sqrt(x.value, self.value)
790	_sig_off
791
792	return x
793
794	def _mpfr_(self, R):
795	return R(str(self)) # TODO: use base 16 or something (!)
796
797
798	def sqrt(self, bits=None):
799	r"""
800	Returns the positive square root of self, possibly as a
801	\emph{a real number} if self is not a perfect integer
802	square.
803
804	INPUT:
805	bits -- number of bits of precision.
806	If bits is not specified, the number of
807	bits of precision is at least twice the
808	number of bits of self (the precision
809	is always at least 53 bits if not specified).
810	OUTPUT:
811	integer, real number, or complex number.
812
813	For the guaranteed integer square root of a perfect square
814	(with error checking), use \code{self.square_root()}.
815
816	EXAMPLE:
817	sage: Z = IntegerRing()
818	sage: Z(2).sqrt(53)
819	1.4142135623730951
820	sage: Z(2).sqrt(100)
821	1.4142135623730950488016887242092
822	sage: 39188072418583779289.square_root()
823	6260037733
824	sage: (100^100).sqrt()
825	10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
826	sage: (-1).sqrt()
827	1.0000000000000000*I
828	sage: sqrt(-2)
829	1.4142135623730949*I
830	sage: sqrt(97)
831	9.8488578017961039
832	sage: 97.sqrt(200)
833	9.8488578017961047217462114149176244816961362874427641717231516
834	"""
835	if bits is None:
836	bits = max(53, 2*(mpz_sizeinbase(self.value, 2)+2))
837	if self < 0:
838	x = sage.rings.complex_field.ComplexField(bits)(self)
839	return x.sqrt()
840	else:
841	try:
842	return self.square_root()
843	except ValueError:
844	pass
845	R = mpfr.RealField(bits)
846	return self._mpfr_(R).sqrt()
847
848	def square_root(self):
849	"""
850	Return the positive integer square root of self, or raises a ValueError
851	if self is not a perfect square.
852	"""
853	n = self.isqrt()
854	if n * n == self:
855	return n
856	raise ValueError, "self (=%s) is not a perfect square"%self
857
858
859	def _xgcd(self, Integer n):
860	"""
861	Return a triple $g, s, t \\in\\Z$ such that
862	$$
863	g = s \\cdot \\mbox{\\rm self} + t \\cdot n.
864	$$
865	"""
866	cdef mpz_t g, s, t
867	cdef object g0, s0, t0
868
869	mpz_init(g)
870	mpz_init(s)
871	mpz_init(t)
872
873	_sig_on
874	mpz_gcdext(g, s, t, self.value, n.value)
875	_sig_off
876
877	g0 = Integer()
878	s0 = Integer()
879	t0 = Integer()
880	set_mpz(g0,g)
881	set_mpz(s0,s)
882	set_mpz(t0,t)
883	mpz_clear(g)
884	mpz_clear(s)
885	mpz_clear(t)
886	return g0, s0, t0
887
888	def __lshift__(Integer self, n):
889	"""
890	"""
891	n = int(n)
892	return self(Integer(2)*n) # todo: optimize
893
894	def __rshift__(Integer self, n):
895	"""
896	"""
897	n = int(n)
898	return self.__floordiv(Integer(2)**n) # todo: optimize
899
900	def _and(Integer self, Integer other):
901	cdef Integer x
902	x = Integer()
903	mpz_and(x.value, self.value, other.value)
904	return x
905
906	def __and__(x, y):
907	if isinstance(x, Integer) and isinstance(y, Integer):
908	return x._and(y)
909	return sage.rings.coerce.bin_op(x, y, operator.and_)
910
911
912	def _or(Integer self, Integer other):
913	cdef Integer x
914	x = Integer()
915	mpz_ior(x.value, self.value, other.value)
916	return x
917
918	def __or__(x, y):
919	if isinstance(x, Integer) and isinstance(y, Integer):
920	return x._or(y)
921	return sage.rings.coerce.bin_op(x, y, operator.or_)
922
923
924	def __invert__(self):
925	"""
926	"""
927	return Integer(1)/self # todo: optimize
928
929
930	def inverse_mod(self, n):
931	"""
932	Returns the inverse of self modulo $n$, if this inverse exists.
933	Otherwise, raises a \exception{ZeroDivisionError} exception.
934
935	INPUT:
936	self -- Integer
937	n -- Integer
938	OUTPUT:
939	x -- Integer such that x*self = 1 (mod m), or
940	raises ZeroDivisionError.
941	IMPLEMENTATION:
942	Call the mpz_invert GMP library function.
943	EXAMPLES:
944	sage: a = Integer(189)
945	sage: a.inverse_mod(10000)
946	4709
947	sage: a.inverse_mod(-10000)
948	4709
949	sage: a.inverse_mod(1890)
950	Traceback (most recent call last):
951	...
952	ZeroDivisionError: Inverse does not exist.
953	sage: a = Integer(19)**100000
954	sage: b = a*a
955	sage: c=a.inverse_mod(b)
956	Traceback (most recent call last):
957	...
958	ZeroDivisionError: Inverse does not exist.
959	"""
960	cdef mpz_t x
961	cdef object ans
962	cdef int r
963	cdef Integer m
964	m = Integer(n)
965
966	if m == 1:
967	return Integer(0)
968
969	mpz_init(x)
970
971	_sig_on
972	r = mpz_invert(x, self.value, m.value)
973	_sig_off
974
975	if r == 0:
976	raise ZeroDivisionError, "Inverse does not exist."
977	ans = Integer()
978	set_mpz(ans,x)
979	mpz_clear(x)
980	return ans
981
982	def _gcd(self, Integer n):
983	"""
984	Return the greatest common divisor of self and $n$.
985	"""
986	cdef mpz_t g
987	cdef object g0
988
989	mpz_init(g)
990
991
992	_sig_on
993	mpz_gcd(g, self.value, n.value)
994	_sig_off
995
996	g0 = Integer()
997	set_mpz(g0,g)
998	mpz_clear(g)
999	return g0
1000
1001	def crt(self, y, m, n):
1002	"""
1003	Return the unique integer between $0$ and $mn$ that is
1004	congruent to the integer modulo $m$ and to $y$ modulo $n$. We
1005	assume that~$m$ and~$n$ are coprime.
1006	"""
1007	cdef object g, s, t
1008	cdef Integer _y, _m, _n
1009	_y = Integer(y); _m = Integer(m); _n = Integer(n)
1010	g, s, t = _m.xgcd(_n)
1011	if not g.is_one():
1012	raise ArithmeticError, "CRT requires that gcd of moduli is 1."
1013	# Now sm + tn = 1, so the answer is x + (y-x)sm, where x=self.
1014	return (self + (_y-self)s_m) % (_m*_n)
1015
1016
1017
1018	ONE = Integer(1)
1019
1020	def integer(x):
1021	if isinstance(x, Integer):
1022	return x
1023	return Integer(x)
1024
1025	def factorial(unsigned long int n):
1026	"""
1027	Return the factorial $n!=1 \\cdot 2 \\cdot 3 \\cdots n$.
1028	The input integer $n$ must fit in an \\code{unsigned long int}.
1029	"""
1030	cdef mpz_t x
1031	cdef Integer z
1032
1033	mpz_init(x)
1034
1035	_sig_on
1036	mpz_fac_ui(x, n)
1037	_sig_off
1038
1039	z = Integer()
1040	set_mpz(z, x)
1041	mpz_clear(x)
1042	return z
1043

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

Download in other formats: