Context Navigation

source: sage/rings/integer_mod.pyx @ 5462:d156f768015e

Visit:

Revision 5462:d156f768015e, 77.3 KB checked in by Robert Bradshaw <robertwb@…>, 6 years ago (diff)
fast integer_mod and rational homomorphisms

Line
1	r"""
2	Elements of $\Z/n\Z$
3
4	An element of the integers modulo $n$.
5
6	There are three types of integer_mod classes, depending on the
7	size of the modulus. The range is capped such that a single
8	arithmetic operation (e.g. multiplication) will not overflow.
9
10	- IntegerMod_int stores its value in a int_fast32_t (typically an int)
11	- IntegerMod_int64 stores its value in a int_fast64_t (typically a long long)
12	- IntegerMod_gmp stores its value in a mpz_t
13
14	All extend IntegerMod_abstract.
15
16	For efficency reasons, it stores the modulus (in all three forms, if possible) in a common (cdef) class NativeIntStruct rather than in the parent.
17
18
19	AUTHORS:
20	-- Robert Bradshaw (most of the work)
21	-- Didier Deshommes (bit shifting)
22	-- William Stein (editing and polishing; new arith architecture)
23	-- Robert Bradshaw (implement native is_square and square_root)
24	-- William Stein (sqrt)
25
26	TESTS:
27	sage: R = Integers(101^3)
28	sage: a = R(824362); b = R(205942)
29	sage: a * b
30	851127
31	"""
32
33	#################################################################################
34	# Copyright (C) 2006 Robert Bradshaw <robertwb@math.washington.edu>
35	# 2006 William Stein <wstein@gmail.com>
36	#
37	# Distributed under the terms of the GNU General Public License (GPL)
38	#
39	# http://www.gnu.org/licenses/
40	#*****************************************************************************
41
42
43	include "../ext/interrupt.pxi" # ctrl-c interrupt block support
44	include "../ext/stdsage.pxi"
45
46	cdef extern from "math.h":
47	double log(double)
48	int ceil(double)
49
50	import operator
51
52
53	import integer_mod_ring
54	import arith
55	import rational
56	from sage.libs.all import pari, PariError
57	import integer_ring
58
59	import commutative_ring_element
60	import sage.interfaces.all
61
62	import sage.rings.integer
63	cimport sage.rings.integer
64	from sage.rings.integer cimport Integer
65
66	import sage.structure.element
67	cimport sage.structure.element
68	from sage.structure.element cimport RingElement, ModuleElement, Element
69	from sage.categories.morphism cimport Morphism
70
71	from sage.structure.parent cimport Parent
72
73
74	def Mod(n, m, parent=None):
75	"""
76	Return the equivalence class of n modulo m as an element of
77	$\Z/m\Z$.
78
79	EXAMPLES:
80	sage: x = Mod(12345678, 32098203845329048)
81	sage: x
82	12345678
83	sage: x^100
84	1017322209155072
85
86	You can also use the lowercase version:
87	sage: mod(12,5)
88	2
89	"""
90	cdef IntegerMod_abstract x
91	x = IntegerMod(integer_mod_ring.IntegerModRing(m), n)
92	if parent is None:
93	return x
94	x._parent = parent
95	return x
96
97
98	mod = Mod
99
100
101	def IntegerMod(parent, value):
102	"""
103	Create an integer modulo $n$ with the given parent.
104
105	This is mainly for internal use.
106	"""
107	cdef NativeIntStruct modulus
108	cdef Py_ssize_t res
109	modulus = parent._pyx_order
110	if modulus.table is not None:
111	if PY_TYPE_CHECK(value, sage.rings.integer.Integer) or PY_TYPE_CHECK(value, int) or PY_TYPE_CHECK(value, long):
112	res = value % modulus.int64
113	if res < 0:
114	res = res + modulus.int64
115	a = modulus.lookup(res)
116	if (<Element>a)._parent is not parent:
117	(<Element>a)._parent = parent
118	# print (<Element>a)._parent, " is not ", parent
119	return a
120	if modulus.int32 != -1:
121	return IntegerMod_int(parent, value)
122	elif modulus.int64 != -1:
123	return IntegerMod_int64(parent, value)
124	else:
125	return IntegerMod_gmp(parent, value)
126
127	def is_IntegerMod(x):
128	"""
129	Return try if and only if x is an integer modulo n.
130
131	EXAMPLES:
132	sage: is_IntegerMod(5)
133	False
134	sage: is_IntegerMod(Mod(5,10))
135	True
136	"""
137	return isinstance(x, IntegerMod_abstract)
138
139	def makeNativeIntStruct(sage.rings.integer.Integer z):
140	return NativeIntStruct(z)
141
142	cdef class NativeIntStruct:
143	"""
144	We store the various forms of the modulus here rather than in the
145	parent for efficiency reasons.
146
147	We may also store a cached table of all elements of a given
148	ring in this class.
149	"""
150	def __init__(NativeIntStruct self, sage.rings.integer.Integer z):
151	self.int64 = -1
152	self.int32 = -1
153	self.table = None # NULL
154	self.sageInteger = z
155	if mpz_cmp_si(z.value, INTEGER_MOD_INT64_LIMIT) < 0:
156	self.int64 = mpz_get_si(z.value)
157	if self.int64 < INTEGER_MOD_INT32_LIMIT:
158	self.int32 = self.int64
159
160	def __reduce__(NativeIntStruct self):
161	return sage.rings.integer_mod.makeNativeIntStruct, (self.sageInteger, )
162
163	def precompute_table(NativeIntStruct self, parent, inverses=True):
164	self.table = PyList_New(self.int64)
165	cdef Py_ssize_t i
166	if self.int32 != -1:
167	for i from 0 <= i < self.int32:
168	z = IntegerMod_int(parent, i)
169	Py_INCREF(z); PyList_SET_ITEM(self.table, i, z)
170	else:
171	for i from 0 <= i < self.int64:
172	z = IntegerMod_int64(parent, i)
173	Py_INCREF(z); PyList_SET_ITEM(self.table, i, z)
174
175	if inverses:
176	tmp = [None] * self.int64
177	for i from 1 <= i < self.int64:
178	try:
179	tmp[i] = ~self.table[i]
180	except ZeroDivisionError:
181	pass
182	self.inverses = tmp
183
184	cdef lookup(NativeIntStruct self, Py_ssize_t value):
185	return <object>PyList_GET_ITEM(self.table, value)
186
187
188	cdef class IntegerMod_abstract(sage.structure.element.CommutativeRingElement):
189
190	def __init__(self, parent):
191	"""
192	EXAMPLES:
193	sage: a = Mod(10,30^10); a
194	10
195	sage: loads(a.dumps()) == a
196	True
197	"""
198	self._parent = parent
199	self.__modulus = parent._pyx_order
200
201
202	cdef _new_c_from_long(self, long value):
203	cdef IntegerMod_abstract x
204	x = <IntegerMod_abstract>PY_NEW(<object>PY_TYPE(self))
205	if PY_TYPE_CHECK(x, IntegerMod_gmp):
206	mpz_init((<IntegerMod_gmp>x).value) # should be done by the new method
207	x._parent = self._parent
208	x.__modulus = self.__modulus
209	x.set_from_long(value)
210	return x
211
212	cdef void set_from_mpz(self, mpz_t value):
213	raise NotImplementedError, "Must be defined in child class."
214
215	cdef void set_from_long(self, long value):
216	raise NotImplementedError, "Must be defined in child class."
217
218	def __abs__(self):
219	"""
220	Raise an error message, since abs(x) makes no sense when x is an
221	integer modulo n.
222
223	EXAMPLES:
224	sage: abs(Mod(2,3))
225	Traceback (most recent call last):
226	...
227	ArithmeticError: absolute valued not defined on integers modulo n.
228	"""
229	raise ArithmeticError, "absolute valued not defined on integers modulo n."
230
231	def __reduce__(IntegerMod_abstract self):
232	"""
233	EXAMPLES:
234	sage: a = Mod(4,5); a
235	4
236	sage: loads(a.dumps()) == a
237	True
238	sage: a = Mod(-1,5^30)^25;
239	sage: loads(a.dumps()) == a
240	True
241	"""
242	return sage.rings.integer_mod.mod, (self.lift(), self.modulus(), self.parent())
243
244	def is_nilpotent(self):
245	r"""
246	Return True if self is nilpotent, i.e., some power of self is zero.
247
248	EXAMPLES:
249	sage: a = Integers(90384098234^3)
250	sage: factor(a.order())
251	2^3 * 191^3 * 236607587^3
252	sage: b = a(2*191)
253	sage: b.is_nilpotent()
254	False
255	sage: b = a(2191236607587)
256	sage: b.is_nilpotent()
257	True
258
259	ALGORITHM: Let $m \geq \log_2(n)$, where $n$ is the modulus.
260	Then $x \in \ZZ/n\ZZ$ is nilpotent if and only if $x^m = 0$.
261
262	PROOF: This is clear if you reduce to the prime power case,
263	which you can do via the Chinese Remainder Theorem.
264
265	We could alternatively factor n and check to see if the prime
266	divisors of n all divide x. This is asymptotically slower :-).
267	"""
268	if self.is_zero():
269	return True
270	m = self.__modulus.sageInteger.exact_log(2) + 1
271	return (self**m).is_zero()
272
273	#################################################################
274	# Interfaces
275	#################################################################
276	def _pari_init_(self):
277	return 'Mod(%s,%s)'%(str(self), self.__modulus.sageInteger)
278
279	def pari(self):
280	return pari(self._pari_init_()) # TODO: is this called implicitly anywhere?
281
282	def _gap_init_(self):
283	r"""
284	Return string representation of corresponding GAP object.
285
286	This can be slow since non-prime GAP finite field elements are
287	represented as powers of a generator for the multiplicative
288	group, so the discrete log problem must be solved.
289
290	\note{This function will create a meaningless GAP object if the
291	modulus is not a power of a prime. Also, the modulus must
292	be $\leq 65536$.}
293
294	EXAMPLES:
295	sage: a = Mod(2,19)
296	sage: gap(a)
297	Z(19)
298	sage: a._gap_(gap)
299	Z(19)
300	sage: gap(a).Int()
301	2
302	sage: b = Mod(0,25)
303	sage: gap(b)
304	0*Z(5)
305	"""
306	R = self.parent()
307	m = self.__modulus.sageInteger
308
309	if m > 65536:
310	raise ValueError, "order must be at most 65536."
311
312	if self == 0:
313	return '0*Z(%s)'%m
314
315	# I couldn't find a guarentee in the GAP docs that the
316	# root of unity they use must be the smallest. This
317	# was not the case in MAGMA once, so who knows, especially
318	# when the order of the ring is not prime. So we make
319	# no such dangerous assumptions (for now).
320
321	# Find the root of unity used by Gap.
322	from sage.interfaces.all import gap # here to reduce dependencies
323	g = int(gap.eval('Int(Z(%s))'%m))
324	n = self.log(R(g))
325	return 'Z(%s)^%s'%(m, n)
326
327	def _magma_init_(self):
328	"""
329	Coercion to Magma.
330
331	EXAMPLES:
332	sage: a = Integers(15)(4)
333	sage: b = magma(a) # optional
334	sage: b.Type() # optional
335	RngIntResElt
336	sage: b^2 # optional
337	1
338	"""
339	return '%s!%s'%(self.parent()._magma_init_(), self)
340
341	def log(self, b=None):
342	r"""
343	Return integer $x$ such that $b^x = a$, where $a$ is self.
344
345	INPUT:
346	self -- unit modulo N
347	b -- a generator of the multiplicative group modulo N.
348	If b is not given, R.multiplicative_generator() is used,
349	where R is the parent of self.
350
351	OUTPUT:
352	Integer $x$ such that $b^x = a$.
353
354	NOTE: The base must not be too big or the current
355	implementation, which is in PARI, will fail.
356
357	EXAMPLES:
358	sage: r = Integers(125)
359	sage: b = r.multiplicative_generator()^3
360	sage: a = b^17
361	sage: a.log(b)
362	17
363	sage: a.log()
364	63
365
366	A bigger example.
367	sage: FF = FiniteField(2^32+61)
368	sage: c = FF(4294967356)
369	sage: x = FF(2)
370	sage: a = c.log(x)
371	sage: a
372	2147483678
373	sage: x^a
374	4294967356
375
376	Things that can go wrong. E.g., if the base is not a
377	generator for the multiplicative group, or not even a unit.
378	You can sometimes use the function \code{discrete_log_generic}
379	in general, but don't expect it to be very fast.
380
381	sage: a = Mod(9, 100); b = Mod(3,100)
382	sage: a.log(b)
383	Traceback (most recent call last):
384	...
385	ValueError: base (=3) for discrete log must generate multiplicative group
386	sage: discrete_log_generic(b^2,b)
387	2
388	sage: a = Mod(16, 100); b = Mod(4,100)
389	sage: a.log(b)
390	Traceback (most recent call last):
391	...
392	ValueError: (8)
393	PARI failed to compute discrete log (perhaps base is not a generator or is too large)
394	sage: discrete_log_generic(a,b)
395	Traceback (most recent call last):
396	...
397	ArithmeticError: multiplicative order of 4 not defined since it is not a unit modulo 100
398
399	AUTHOR:
400	-- David Joyner and William Stein (2005-11)
401	-- William Stein (2007-01-27): update to use PARI as requested by
402	David Kohel.
403	"""
404	if b is None:
405	b = self._parent.multiplicative_generator()
406	else:
407	b = self._parent(b)
408	cmd = 'b=Mod(%s,%s); if(znorder(b)!=eulerphi(%s),-1,znlog(%s,b))'%(b, self.__modulus.sageInteger,
409	self.__modulus.sageInteger, self)
410	try:
411	n = Integer(pari(cmd))
412	if n == -1:
413	raise ValueError, "base (=%s) for discrete log must generate multiplicative group"%b
414	return n
415	except PariError, msg:
416	raise ValueError, "%s\nPARI failed to compute discrete log (perhaps base is not a generator or is too large)"%msg
417
418
419	def modulus(IntegerMod_abstract self):
420	"""
421	EXAMPLES:
422	sage: Mod(3,17).modulus()
423	17
424	"""
425	return self.__modulus.sageInteger
426
427	def charpoly(self, var):
428	"""
429	Returns the characteristic polynomial of this element.
430
431	EXAMPLES:
432	sage: k = GF(3)
433	sage: a = k.gen()
434	sage: a.charpoly('x')
435	x + 2
436	sage: a + 2
437	0
438
439	AUTHOR:
440	-- Craig Citro
441	"""
442	from polynomial.polynomial_ring import PolynomialRing
443	R = PolynomialRing(self._parent, var)
444	return R([-self,1])
445
446	def norm(self):
447	"""
448	Returns the norm of this element, which is itself. (This
449	is here for compatibility with higher order finite fields.)
450
451	EXAMPLES:
452	sage: k = GF(691)
453	sage: a = k(389)
454	sage: a.norm()
455	389
456
457	AUTHOR:
458	-- Craig Citro
459	"""
460	return self
461
462	def trace(self):
463	"""
464	Returns the trace of this element, which is itself. (This
465	is here for compatibility with higher order finite fields.)
466
467	EXAMPLES:
468	sage: k = GF(691)
469	sage: a = k(389)
470	sage: a.trace()
471	389
472
473	AUTHOR:
474	-- Craig Citro
475	"""
476	return self
477
478
479	def is_square(self):
480	r"""
481	EXAMPLES:
482	sage: Mod(3,17).is_square()
483	False
484	sage: Mod(9,17).is_square()
485	True
486	sage: Mod(9,17*19^2).is_square()
487	True
488	sage: Mod(-1,17^30).is_square()
489	True
490	sage: Mod(1/9, next_prime(2^40)).is_square()
491	True
492	sage: Mod(1/25, next_prime(2^90)).is_square()
493	True
494
495	TESTS:
496	sage: Mod(1/25, 2^8).is_square()
497	True
498	sage: Mod(1/25, 2^40).is_square()
499	True
500
501	ALGORITHM:
502	Calculate the Jacobi symbol $(self/p)$ at each prime $p$ dividing $n$
503	It must be 1 or 0 for each prime, and if it is 0 mod $p$,
504	where $p^k \|\| n$, then $ord_p(self)$ must be even or greater than $k$.
505
506	$p = 2$ handled seperatly.
507
508	AUTHOR:
509	-- Robert Bradshaw
510	"""
511	return self.is_square_c()
512
513	cdef bint is_square_c(self) except -2:
514	if self.is_zero() or self.is_one():
515	return 1
516	moduli = self.parent().factored_order()
517	cdef int val, e
518	lift = self.lift()
519	if len(moduli) == 1:
520	p, e = moduli[0]
521	if e == 1:
522	return lift.jacobi(p) != -1
523	elif p == 2:
524	return self.pari().issquare() # TODO: implement directly
525	elif self % p == 0:
526	val = lift.valuation(p)
527	return val >= e or (val % 2 == 0 and (lift // p**val).jacobi(p) != -1)
528	else:
529	return lift.jacobi(p) != -1
530	else:
531	for p, e in moduli:
532	if p == 2:
533	if e > 1 and not self.pari().issquare(): # TODO: implement directly
534	return 0
535	elif e > 1 and lift % p == 0:
536	val = lift.valuation(p)
537	if val < e and (val % 2 == 1 or (lift // p**val).jacobi(p) == -1):
538	return 0
539	elif lift.jacobi(p) == -1:
540	return 0
541	return 1
542
543	def sqrt(self, extend=True, all=False):
544	r"""
545	Returns square root or square roots of self modulo n.
546
547	INPUT:
548	extend -- bool (default: True); if True, return a square
549	root in an extension ring, if necessary. Otherwise,
550	raise a ValueError if the square is not in the base ring.
551	all -- bool (default: False); if True, return all square
552	roots of self, instead of just one.
553
554	ALGORITHM: Calculates the square roots mod $p$ for each of the
555	primes $p$ dividing the order of the ring, then lifts them
556	p-adically and uses the CRT to find a square root mod $n$.
557
558	See also \code{square_root_mod_prime_power} and
559	\code{square_root_mod_prime} (in this module) for more
560	algorithmic details.
561
562	EXAMPLES:
563	sage: mod(-1, 17).sqrt()
564	4
565	sage: mod(5, 389).sqrt()
566	86
567	sage: mod(7, 18).sqrt()
568	5
569	sage: a = mod(14, 5^60).sqrt()
570	sage: a*a
571	14
572	sage: mod(15, 389).sqrt(extend=False)
573	Traceback (most recent call last):
574	...
575	ValueError: self must be a square
576	sage: Mod(1/9, next_prime(2^40)).sqrt()^(-2)
577	9
578	sage: Mod(1/25, next_prime(2^90)).sqrt()^(-2)
579	25
580
581	sage: a = Mod(3,5); a
582	3
583	sage: x = Mod(-1, 360)
584	sage: x.sqrt(extend=False)
585	Traceback (most recent call last):
586	...
587	ValueError: self must be a square
588	sage: y = x.sqrt(); y
589	sqrt359
590	sage: y.parent()
591	Univariate Quotient Polynomial Ring in sqrt359 over Ring of integers modulo 360 with modulus x^2 + 1
592	sage: y^2
593	359
594
595	We compute all square roots in several cases:
596	sage: R = Integers(52^33^2); R
597	Ring of integers modulo 360
598	sage: R(40).sqrt(all=True)
599	[20, 160, 200, 340]
600	sage: [x for x in R if x^2 == 40] # Brute force verification
601	[20, 160, 200, 340]
602	sage: R(1).sqrt(all=True)
603	[1, 19, 71, 89, 91, 109, 161, 179, 181, 199, 251, 269, 271, 289, 341, 359]
604	sage: R(0).sqrt(all=True)
605	[0, 60, 120, 180, 240, 300]
606
607	sage: R = Integers(513^337); R
608	Ring of integers modulo 406445
609	sage: v = R(-1).sqrt(all=True); v
610	[78853, 111808, 160142, 193097, 213348, 246303, 294637, 327592]
611	sage: [x^2 for x in v]
612	[406444, 406444, 406444, 406444, 406444, 406444, 406444, 406444]
613	sage: v = R(169).sqrt(all=True); min(v), -max(v), len(v)
614	(13, 13, 104)
615	sage: all([x^2==169 for x in v])
616	True
617
618	Modulo a power of 2:
619	sage: R = Integers(2^7); R
620	Ring of integers modulo 128
621	sage: a = R(17)
622	sage: a.sqrt()
623	23
624	sage: a.sqrt(all=True)
625	[23, 41, 87, 105]
626	sage: [x for x in R if x^2==17]
627	[23, 41, 87, 105]
628
629	"""
630	if self.is_one():
631	if all:
632	return list(self.parent().square_roots_of_one())
633	else:
634	return self
635
636	if not self.is_square_c():
637	if extend:
638	y = 'sqrt%s'%self
639	R = self.parent()['x']
640	modulus = R.gen()**2 - R(self)
641	if self._parent.is_field():
642	import finite_field
643	Q = finite_field.FiniteField(self.__modulus.sageInteger**2, y, modulus)
644	else:
645	R = self.parent()['x']
646	Q = R.quotient(modulus, names=(y,))
647	z = Q.gen()
648	if all:
649	# TODO
650	raise NotImplementedError
651	return z
652	raise ValueError, "self must be a square"
653
654	F = self._parent.factored_order()
655	cdef long e, exp, val
656	if len(F) == 1:
657	p, e = F[0]
658
659	if all and e > 1 and not self.is_unit():
660	if self.is_zero():
661	# All multiples of p^ciel(e/2) vanish
662	return [self._parent(x) for x in xrange(0, self.__modulus.sageInteger, p**((e+1)/2))]
663	else:
664	z = self.lift()
665	val = z.valuation(p)/2 # square => valuation is even
666	from integer_mod_ring import IntegerModRing
667	# Find the unit part (mod the ring with appropriate precision)
668	u = IntegerModRing(p(e-val))(z // p(2*val))
669	# will add multiples of p^exp
670	exp = e - val
671	if p == 2:
672	exp -= 1 # note the factor of 2 below
673	if 2*exp < e:
674	exp = (e+1)/2
675	# For all a^2 = u and all integers b
676	# (ap^val + bp^exp) ^ 2
677	# = up^(2val) + 2abp^(val+exp) + b^2p^(2*exp)
678	# = up^(2val) mod p^e
679	# whenever min(val+exp, 2*exp) > e
680	p_val = p**val
681	p_exp = p**exp
682	w = [self._parent(a.lift() * p_val + b)
683	for a in u.sqrt(all=True)
684	for b in xrange(0, self.__modulus.sageInteger, p_exp)]
685	if p == 2:
686	w = list(set(w))
687	w.sort()
688	return w
689
690	if e > 1:
691	x = square_root_mod_prime_power(mod(self, p**e), p, e)
692	else:
693	x = square_root_mod_prime(self, p)
694	x = x._balanced_abs()
695
696	if not all:
697	return x
698
699	v = list(set([x*a for a in self._parent.square_roots_of_one()]))
700	v.sort()
701	return v
702
703	else:
704	if not all:
705	# Use CRT to combine together a square root modulo each prime power
706	sqrts = [square_root_mod_prime(mod(self, p), p) for p, e in F if e == 1] + \
707	[square_root_mod_prime_power(mod(self, p**e), p, e) for p, e in F if e != 1]
708
709	x = sqrts.pop()
710	for y in sqrts:
711	x = x.crt(y)
712	return x._balanced_abs()
713	else:
714	# Use CRT to combine together all square roots modulo each prime power
715	vmod = []
716	moduli = []
717	P = self.parent()
718	from integer_mod_ring import IntegerModRing
719	for p, e in F:
720	k = p**e
721	R = IntegerModRing(p**e)
722	w = [P(x) for x in R(self).sqrt(all=True)]
723	vmod.append(w)
724	moduli.append(k)
725	# Now combine in all possible ways using the CRT
726	from arith import CRT_basis
727	basis = CRT_basis(moduli)
728	from sage.misc.mrange import cartesian_product_iterator
729	v = []
730	for x in cartesian_product_iterator(vmod):
731	# x is a specific choice of roots modulo each prime power divisor
732	a = sum([basis[i]*x[i] for i in range(len(x))])
733	v.append(a)
734	v.sort()
735	return v
736
737	def square_root(self):
738	return self.sqrt()
739
740	def _balanced_abs(self):
741	"""
742	This function returns x or -x, whichever has a positive
743	representative in -p/2 < x < p/2.
744
745	This is used so that the same square root is always returned,
746	despite the possibly probabalistic nature of the underlying
747	algorithm.
748	"""
749	if self.lift() > self.__modulus.sageInteger >> 1:
750	return -self
751	else:
752	return self
753
754
755	def rational_reconstruction(self):
756	"""
757	EXAMPLES:
758	sage: R = IntegerModRing(97)
759	sage: a = R(2) / R(3)
760	sage: a
761	33
762	sage: a.rational_reconstruction()
763	2/3
764	"""
765	return self.lift().rational_reconstruction(self.modulus())
766
767	def crt(IntegerMod_abstract self, IntegerMod_abstract other):
768	"""
769	Use the Chinese Remainder Theorem to find an element of the
770	integers modulo the product of the moduli that reduces to self
771	and to other. The modulus of other must be coprime to the
772	modulus of self.
773	EXAMPLES:
774	sage: a = mod(3,5)
775	sage: b = mod(2,7)
776	sage: a.crt(b)
777	23
778
779	sage: a = mod(37,10^8)
780	sage: b = mod(9,3^8)
781	sage: a.crt(b)
782	125900000037
783
784	AUTHOR:
785	-- Robert Bradshaw
786	"""
787	cdef int_fast64_t new_modulus
788	if not isinstance(self, IntegerMod_gmp) and not isinstance(other, IntegerMod_gmp):
789
790	if other.__modulus.int64 == 1: return self
791	new_modulus = self.__modulus.int64 * other.__modulus.int64
792	if new_modulus < INTEGER_MOD_INT32_LIMIT:
793	return self.__crt(other)
794
795	elif new_modulus < INTEGER_MOD_INT64_LIMIT:
796	if not isinstance(self, IntegerMod_int64):
797	self = IntegerMod_int64(self._parent, self.lift())
798	if not isinstance(other, IntegerMod_int64):
799	other = IntegerMod_int64(other._parent, other.lift())
800	return self.__crt(other)
801
802	if not isinstance(self, IntegerMod_gmp):
803	self = IntegerMod_gmp(self._parent, self.lift())
804
805	if not isinstance(other, IntegerMod_gmp):
806	other = IntegerMod_gmp(other._parent, other.lift())
807
808	return self.__crt(other)
809
810
811	def additive_order(self):
812	r"""
813	Returns the additive order of self.
814
815	This is the same as \code{self.order()}.
816
817	EXAMPLES:
818	sage: Integers(20)(2).additive_order()
819	10
820	sage: Integers(20)(7).additive_order()
821	20
822	sage: Integers(90308402384902)(2).additive_order()
823	45154201192451
824	"""
825	n = self.__modulus.sageInteger
826	return sage.rings.integer.Integer(n.__floordiv__(self.lift().gcd(n)))
827
828	def multiplicative_order(self):
829	"""
830	Returns the multiplicative order of self.
831
832	EXAMPLES:
833	sage: Mod(-1,5).multiplicative_order()
834	2
835	sage: Mod(1,5).multiplicative_order()
836	1
837	sage: Mod(0,5).multiplicative_order()
838	Traceback (most recent call last):
839	...
840	ArithmeticError: multiplicative order of 0 not defined since it is not a unit modulo 5
841	"""
842	try:
843	return sage.rings.integer.Integer(self.pari().order()) # pari's "order" is by default multiplicative
844	except PariError:
845	raise ArithmeticError, "multiplicative order of %s not defined since it is not a unit modulo %s"%(
846	self, self.__modulus.sageInteger)
847
848	def _repr_(self):
849	return str(self.lift())
850
851	def _latex_(self):
852	return str(self)
853
854	def _integer_(self):
855	return self.lift()
856
857	def _rational_(self):
858	return rational.Rational(self.lift())
859
860
861
862
863	######################################################################
864	# class IntegerMod_gmp
865	######################################################################
866
867
868	cdef class IntegerMod_gmp(IntegerMod_abstract):
869	"""
870	Elements of $\Z/n\Z$ for n not small enough to be operated on in word size
871	AUTHORS:
872	-- Robert Bradshaw (2006-08-24)
873	"""
874
875	def __init__(IntegerMod_gmp self, parent, value, empty=False):
876	"""
877	EXAMPLES:
878	sage: a = mod(5,14^20)
879	sage: type(a)
880	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
881	sage: loads(dumps(a)) == a
882	True
883	"""
884	mpz_init(self.value)
885	IntegerMod_abstract.__init__(self, parent)
886	if empty:
887	return
888	cdef sage.rings.integer.Integer z
889	if PY_TYPE_CHECK(value, sage.rings.integer.Integer):
890	z = value
891	elif PY_TYPE_CHECK(value, rational.Rational):
892	z = value % self.__modulus.sageInteger
893	elif PY_TYPE_CHECK(value, int):
894	self.set_from_long(value)
895	return
896	else:
897	z = sage.rings.integer_ring.Z(value)
898	self.set_from_mpz(z.value)
899
900	cdef IntegerMod_gmp _new_c(self):
901	cdef IntegerMod_gmp x
902	x = PY_NEW(IntegerMod_gmp)
903	mpz_init(x.value)
904	x.__modulus = self.__modulus
905	x._parent = self._parent
906	return x
907
908	def __dealloc__(self):
909	mpz_clear(self.value)
910
911	cdef void set_from_mpz(self, mpz_t value):
912	cdef sage.rings.integer.Integer modulus
913	modulus = self.__modulus.sageInteger
914	if mpz_sgn(value) == -1 or mpz_cmp(value, modulus.value) >= 0:
915	mpz_mod(self.value, value, modulus.value)
916	else:
917	mpz_set(self.value, value)
918
919	cdef void set_from_long(self, long value):
920	cdef sage.rings.integer.Integer modulus
921	mpz_set_si(self.value, value)
922	if value < 0 or mpz_cmp_si(self.__modulus.sageInteger.value, value) >= 0:
923	mpz_mod(self.value, self.value, self.__modulus.sageInteger.value)
924
925	cdef mpz_t* get_value(IntegerMod_gmp self):
926	return &self.value
927
928	def __lshift__(IntegerMod_gmp self, int right):
929	r"""
930	Multiply self by $2^\text{right}$ very quickly via bit shifting.
931
932	EXAMPLES:
933	sage: e = Mod(19, 10^10)
934	sage: e << 102
935	9443608576
936	"""
937	cdef IntegerMod_gmp x
938	x = self._new_c()
939	mpz_mul_2exp(x.value, self.value, right)
940	mpz_fdiv_r(x.value, x.value, self.__modulus.sageInteger.value)
941	return x
942
943	cdef int _cmp_c_impl(left, Element right) except -2:
944	"""
945	EXAMPLES:
946	sage: mod(5,13^20) == mod(5,13^20)
947	True
948	sage: mod(5,13^20) == mod(-5,13^20)
949	False
950	sage: mod(5,13^20) == mod(-5,13)
951	False
952	"""
953	cdef int i
954	i = mpz_cmp((<IntegerMod_gmp>left).value, (<IntegerMod_gmp>right).value)
955	if i < 0:
956	return -1
957	elif i == 0:
958	return 0
959	else:
960	return 1
961
962	def __richcmp__(left, right, int op):
963	return (<Element>left)._richcmp(right, op)
964
965
966	def is_one(IntegerMod_gmp self):
967	"""
968	Returns \code{True} if this is $1$, otherwise \code{False}.
969
970	EXAMPLES:
971	sage: mod(1,5^23).is_one()
972	True
973	sage: mod(0,5^23).is_one()
974	False
975	"""
976	return mpz_cmp_si(self.value, 1) == 0
977
978	def __nonzero__(IntegerMod_gmp self):
979	"""
980	Returns \code{True} if this is not $0$, otherwise \code{False}.
981
982	EXAMPLES:
983	sage: mod(13,5^23).is_zero()
984	False
985	sage: (mod(25,5^23)^23).is_zero()
986	True
987	"""
988	return mpz_cmp_si(self.value, 0) != 0
989
990	def is_unit(self):
991	return self.lift().gcd(self.modulus()) == 1
992
993	def __crt(IntegerMod_gmp self, IntegerMod_gmp other):
994	cdef IntegerMod_gmp lift, x
995	cdef sage.rings.integer.Integer modulus, other_modulus
996
997	modulus = self.__modulus.sageInteger
998	other_modulus = other.__modulus.sageInteger
999	lift = IntegerMod_gmp(integer_mod_ring.IntegerModRing(modulus*other_modulus), None, empty=True)
1000	try:
1001	if mpz_cmp(self.value, other.value) > 0:
1002	x = (other - IntegerMod_gmp(other._parent, self.lift())) / IntegerMod_gmp(other._parent, modulus)
1003	mpz_mul(lift.value, x.value, modulus.value)
1004	mpz_add(lift.value, lift.value, self.value)
1005	else:
1006	x = (self - IntegerMod_gmp(self._parent, other.lift())) / IntegerMod_gmp(self._parent, other_modulus)
1007	mpz_mul(lift.value, x.value, other_modulus.value)
1008	mpz_add(lift.value, lift.value, other.value)
1009	return lift
1010	except ZeroDivisionError:
1011	raise ZeroDivisionError, "moduli must be coprime"
1012
1013
1014	def __copy__(IntegerMod_gmp self):
1015	cdef IntegerMod_gmp x
1016	x = self._new_c()
1017	mpz_set(x.value, self.value)
1018
1019	cdef ModuleElement _add_c_impl(self, ModuleElement right):
1020	"""
1021	EXAMPLES:
1022	sage: R = Integers(10^10)
1023	sage: R(7) + R(8)
1024	15
1025	"""
1026	cdef IntegerMod_gmp x
1027	x = self._new_c()
1028	mpz_add(x.value, self.value, (<IntegerMod_gmp>right).value)
1029	if mpz_cmp(x.value, self.__modulus.sageInteger.value) >= 0:
1030	mpz_sub(x.value, x.value, self.__modulus.sageInteger.value)
1031	return x;
1032
1033	cdef ModuleElement _sub_c_impl(self, ModuleElement right):
1034	"""
1035	EXAMPLES:
1036	sage: R = Integers(10^10)
1037	sage: R(7) - R(8)
1038	9999999999
1039	"""
1040	cdef IntegerMod_gmp x
1041	x = self._new_c()
1042	mpz_sub(x.value, self.value, (<IntegerMod_gmp>right).value)
1043	if mpz_sgn(x.value) == -1:
1044	mpz_add(x.value, x.value, self.__modulus.sageInteger.value)
1045	return x;
1046
1047	cdef ModuleElement _neg_c_impl(self):
1048	"""
1049	EXAMPLES:
1050	sage: -mod(5,10^10)
1051	9999999995
1052	sage: -mod(0,10^10)
1053	0
1054	"""
1055	if mpz_cmp_si(self.value, 0) == 0:
1056	return self
1057	cdef IntegerMod_gmp x
1058	x = self._new_c()
1059	mpz_sub(x.value, self.__modulus.sageInteger.value, self.value)
1060	return x
1061
1062	cdef RingElement _mul_c_impl(self, RingElement right):
1063	"""
1064	EXAMPLES:
1065	sage: R = Integers(10^11)
1066	sage: R(700000) * R(800000)
1067	60000000000
1068	"""
1069	cdef IntegerMod_gmp x
1070	x = self._new_c()
1071	mpz_mul(x.value, self.value, (<IntegerMod_gmp>right).value)
1072	mpz_fdiv_r(x.value, x.value, self.__modulus.sageInteger.value)
1073	return x
1074
1075	cdef RingElement _div_c_impl(self, RingElement right):
1076	"""
1077	EXAMPLES:
1078	sage: R = Integers(10^11)
1079	sage: R(3) / R(7)
1080	71428571429
1081	"""
1082	return self._mul_c(~right)
1083
1084	def __int__(self):
1085	return int(self.lift())
1086
1087	def __index__(self):
1088	"""
1089	Needed so integers modulo n can be used as list indices.
1090
1091	EXAMPLES:
1092	sage: v = [1,2,3,4,5]
1093	sage: v[Mod(3,10^20)]
1094	4
1095	"""
1096	return int(self.lift())
1097
1098	def __long__(self):
1099	return long(self.lift())
1100
1101	def __mod__(self, right):
1102	if self.modulus() % right != 0:
1103	raise ZeroDivisionError, "Error - reduction modulo right not defined."
1104	return IntegerMod(integer_mod_ring.IntegerModRing(right), self)
1105
1106	def __pow__(IntegerMod_gmp self, right, m): # NOTE: m ignored, always use modulus of parent ring
1107	"""
1108	EXAMPLES:
1109	sage: R = Integers(10^10)
1110	sage: R(2)^1000
1111	5668069376
1112	sage: p = next_prime(11^10)
1113	sage: R = Integers(p)
1114	sage: R(9876)^(p-1)
1115	1
1116	"""
1117	cdef sage.rings.integer.Integer exp
1118	exp = sage.rings.integer_ring.Z(right)
1119	cdef IntegerMod_gmp x
1120	x = self._new_c()
1121	_sig_on
1122	mpz_powm(x.value, self.value, exp.value, self.__modulus.sageInteger.value)
1123	_sig_off
1124	return x
1125
1126	def __rshift__(IntegerMod_gmp self, int right):
1127	r"""
1128	Divide self by $2^{\text{right}}$ and take floor via bit shifting.
1129
1130	EXAMPLES:
1131	sage: e = Mod(1000001, 2^32-1)
1132	sage: e >> 5
1133	31250
1134	"""
1135	cdef IntegerMod_gmp x
1136	x = self._new_c()
1137	mpz_fdiv_q_2exp(x.value, self.value, right)
1138	return x
1139
1140	def __invert__(IntegerMod_gmp self):
1141	"""
1142	Return the multiplicative inverse of self.
1143
1144	EXAMPLES:
1145	sage: a = mod(3,10^100); type(a)
1146	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
1147	sage: ~a
1148	6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
1149	sage: ~mod(2,10^100)
1150	Traceback (most recent call last):
1151	...
1152	ZeroDivisionError: Inverse does not exist.
1153	"""
1154	if self.is_zero():
1155	raise ZeroDivisionError, "Inverse does not exist."
1156
1157	cdef IntegerMod_gmp x
1158	x = self._new_c()
1159	if (mpz_invert(x.value, self.value, self.__modulus.sageInteger.value)):
1160	return x
1161	else:
1162	raise ZeroDivisionError, "Inverse does not exist."
1163
1164	def lift(IntegerMod_gmp self):
1165	"""
1166	Lift an integer modulo $n$ to the integers.
1167
1168	EXAMPLES:
1169	sage: a = Mod(8943, 2^70); type(a)
1170	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
1171	sage: lift(a)
1172	8943
1173	sage: a.lift()
1174	8943
1175	"""
1176	cdef sage.rings.integer.Integer z
1177	z = sage.rings.integer.Integer()
1178	z.set_from_mpz(self.value)
1179	return z
1180
1181	def __float__(self):
1182	return float(self.lift())
1183
1184	def __hash__(self):
1185	"""
1186	EXAMPLES:
1187	sage: a = Mod(8943, 2^100)
1188	sage: hash(a)
1189	8943
1190	"""
1191	# return mpz_pythonhash(self.value)
1192	return hash(self.lift())
1193
1194
1195
1196	######################################################################
1197	# class IntegerMod_int
1198	######################################################################
1199
1200
1201	cdef class IntegerMod_int(IntegerMod_abstract):
1202	"""
1203	Elements of $\Z/n\Z$ for n small enough to be operated on in 32 bits
1204	AUTHORS:
1205	-- Robert Bradshaw (2006-08-24)
1206	"""
1207
1208	def __init__(self, parent, value, empty=False):
1209	"""
1210	EXAMPLES:
1211	sage: a = Mod(10,30); a
1212	10
1213	sage: loads(a.dumps()) == a
1214	True
1215	"""
1216	IntegerMod_abstract.__init__(self, parent)
1217	if empty:
1218	return
1219	cdef int_fast32_t x
1220	if PY_TYPE_CHECK(value, int):
1221	x = value
1222	self.ivalue = x % self.__modulus.int32
1223	if self.ivalue < 0:
1224	self.ivalue = self.ivalue + self.__modulus.int32
1225	return
1226	cdef sage.rings.integer.Integer z
1227	if PY_TYPE_CHECK(value, sage.rings.integer.Integer):
1228	z = value
1229	elif isinstance(value, rational.Rational):
1230	z = value % self.__modulus.sageInteger
1231	else:
1232	z = sage.rings.integer_ring.Z(value)
1233	self.set_from_mpz(z.value)
1234
1235	cdef IntegerMod_int _new_c(self, int_fast32_t value):
1236	if self.__modulus.table is not None:
1237	return self.__modulus.lookup(value)
1238	cdef IntegerMod_int x = PY_NEW(IntegerMod_int)
1239	x._parent = self._parent
1240	x.__modulus = self.__modulus
1241	x.ivalue = value
1242	return x
1243
1244	cdef void set_from_mpz(self, mpz_t value):
1245	if mpz_sgn(value) == -1 or mpz_cmp_si(value, self.__modulus.int32) >= 0:
1246	self.ivalue = mpz_fdiv_ui(value, self.__modulus.int32)
1247	else:
1248	self.ivalue = mpz_get_si(value)
1249
1250	cdef void set_from_long(self, long value):
1251	self.ivalue = value % self.__modulus.int32
1252
1253	cdef void set_from_int(IntegerMod_int self, int_fast32_t ivalue):
1254	if ivalue < 0:
1255	self.ivalue = self.__modulus.int32 + (ivalue % self.__modulus.int32)
1256	elif ivalue >= self.__modulus.int32:
1257	self.ivalue = ivalue % self.__modulus.int32
1258	else:
1259	self.ivalue = ivalue
1260
1261	cdef int_fast32_t get_int_value(IntegerMod_int self):
1262	return self.ivalue
1263
1264
1265
1266	cdef int _cmp_c_impl(left, Element right) except -2:
1267	"""
1268	EXAMPLES:
1269	sage: mod(5,13) == mod(-8,13)
1270	True
1271	sage: mod(5,13) == mod(8,13)
1272	False
1273	sage: mod(5,13) == mod(5,24)
1274	False
1275	sage: mod(0, 13) == 0
1276	True
1277	sage: mod(0, 13) == int(0)
1278	True
1279	"""
1280	if (<IntegerMod_int>left).ivalue == (<IntegerMod_int>right).ivalue:
1281	return 0
1282	elif (<IntegerMod_int>left).ivalue < (<IntegerMod_int>right).ivalue:
1283	return -1
1284	else:
1285	return 1
1286
1287	def __richcmp__(left, right, int op):
1288	return (<Element>left)._richcmp(right, op)
1289
1290
1291	def is_one(IntegerMod_int self):
1292	"""
1293	Returns \code{True} if this is $1$, otherwise \code{False}.
1294
1295	EXAMPLES:
1296	sage: mod(6,5).is_one()
1297	True
1298	sage: mod(0,5).is_one()
1299	False
1300	"""
1301	return self.ivalue == 1
1302
1303	def __nonzero__(IntegerMod_int self):
1304	"""
1305	Returns \code{True} if this is not $0$, otherwise \code{False}.
1306
1307	EXAMPLES:
1308	sage: mod(13,5).is_zero()
1309	False
1310	sage: mod(25,5).is_zero()
1311	True
1312	"""
1313	return self.ivalue != 0
1314
1315	def is_unit(IntegerMod_int self):
1316	return gcd_int(self.ivalue, self.__modulus.int32) == 1
1317
1318	def __crt(IntegerMod_int self, IntegerMod_int other):
1319	"""
1320	Use the Chinese Remainder Theorem to find an element of the
1321	integers modulo the product of the moduli that reduces to self
1322	and to other. The modulus of other must be coprime to the
1323	modulus of self.
1324	EXAMPLES:
1325	sage: a = mod(3,5)
1326	sage: b = mod(2,7)
1327	sage: a.crt(b)
1328	23
1329
1330	AUTHOR:
1331	-- Robert Bradshaw
1332	"""
1333	cdef IntegerMod_int lift
1334	cdef int_fast32_t x
1335
1336	lift = IntegerMod_int(integer_mod_ring.IntegerModRing(self.__modulus.int32 * other.__modulus.int32), None, empty=True)
1337
1338	try:
1339	x = (other.ivalue - self.ivalue % other.__modulus.int32) * mod_inverse_int(self.__modulus.int32, other.__modulus.int32)
1340	lift.set_from_int( x * self.__modulus.int32 + self.ivalue )
1341	return lift
1342	except ZeroDivisionError:
1343	raise ZeroDivisionError, "moduli must be coprime"
1344
1345
1346	def __copy__(IntegerMod_int self):
1347	return self._new_c(self.ivalue)
1348
1349	cdef ModuleElement _add_c_impl(self, ModuleElement right):
1350	"""
1351	EXAMPLES:
1352	sage: R = Integers(10)
1353	sage: R(7) + R(8)
1354	5
1355	"""
1356	cdef int_fast32_t x
1357	x = self.ivalue + (<IntegerMod_int>right).ivalue
1358	if x >= self.__modulus.int32:
1359	x = x - self.__modulus.int32
1360	return self._new_c(x)
1361
1362	cdef ModuleElement _sub_c_impl(self, ModuleElement right):
1363	"""
1364	EXAMPLES:
1365	sage: R = Integers(10)
1366	sage: R(7) - R(8)
1367	9
1368	"""
1369	cdef int_fast32_t x
1370	x = self.ivalue - (<IntegerMod_int>right).ivalue
1371	if x < 0:
1372	x = x + self.__modulus.int32
1373	return self._new_c(x)
1374
1375	cdef ModuleElement _neg_c_impl(self):
1376	"""
1377	EXAMPLES:
1378	sage: -mod(7,10)
1379	3
1380	sage: -mod(0,10)
1381	0
1382	"""
1383	if self.ivalue == 0:
1384	return self
1385	return self._new_c(self.__modulus.int32 - self.ivalue)
1386
1387	cdef RingElement _mul_c_impl(self, RingElement right):
1388	"""
1389	EXAMPLES:
1390	sage: R = Integers(10)
1391	sage: R(7) * R(8)
1392	6
1393	"""
1394	return self._new_c((self.ivalue * right.ivalue) % self.__modulus.int32)
1395
1396	cdef RingElement _div_c_impl(self, RingElement right):
1397	"""
1398	EXAMPLES:
1399	sage: R = Integers(10)
1400	sage: R(2)/3
1401	4
1402	"""
1403	if self.__modulus.inverses is not None:
1404	right_inverse = self.__modulus.inverses[(<IntegerMod_int>right).ivalue]
1405	if right_inverse is None:
1406	raise ZeroDivisionError, "Inverse does not exist."
1407	else:
1408	return self._new_c((self.ivalue * (<IntegerMod_int>right_inverse).ivalue) % self.__modulus.int32)
1409
1410	cdef int_fast32_t x
1411	x = self.ivalue * mod_inverse_int((<IntegerMod_int>right).ivalue, self.__modulus.int32)
1412	return self._new_c(x% self.__modulus.int32)
1413
1414	def __int__(IntegerMod_int self):
1415	return self.ivalue
1416
1417	def __index__(self):
1418	"""
1419	Needed so integers can be used as list indices.
1420
1421	EXAMPLES:
1422	sage: v = [1,2,3,4,5]
1423	sage: v[Mod(10,7)]
1424	4
1425	"""
1426	return self.ivalue
1427
1428	def __long__(IntegerMod_int self):
1429	return self.ivalue
1430
1431	def __mod__(IntegerMod_int self, right):
1432	right = int(right)
1433	if self.__modulus.int32 % right != 0:
1434	raise ZeroDivisionError, "reduction modulo right not defined."
1435	return integer_mod_ring.IntegerModRing(right)(self)
1436
1437	def __lshift__(IntegerMod_int self, int right):
1438	r"""
1439	Multiply self by $2^\text{right}$ very quickly via bit shifting.
1440
1441	EXAMPLES:
1442	sage: e = Mod(5, 2^10 - 1)
1443	sage: e<<5
1444	160
1445	sage: e*2^5
1446	160
1447	"""
1448	return self._new_c((self.ivalue << right) % self.__modulus.int32)
1449
1450	def __rshift__(IntegerMod_int self, int right):
1451	"""
1452	Divide self by $2^{\text{right}}$ and take floor via bit shifting.
1453
1454	EXAMPLES:
1455	sage: e = Mod(8, 2^5 - 1)
1456	sage: e >> 3
1457	1
1458	sage: int(e)/int(2^3)
1459	1
1460	"""
1461	return self._new_c(self.ivalue >> right)
1462
1463	def __pow__(IntegerMod_int self, right, m): # NOTE: m ignored, always use modulus of parent ring
1464	"""
1465	EXAMPLES:
1466	sage: R = Integers(10)
1467	sage: R(2)^10
1468	4
1469	sage: R = Integers(389)
1470	sage: R(7)^388
1471	1
1472	"""
1473	cdef sage.rings.integer.Integer exp, base
1474	exp = sage.rings.integer_ring.Z(right)
1475	cdef int_fast32_t x
1476	cdef mpz_t x_mpz
1477	if mpz_sgn(exp.value) >= 0 and mpz_cmp_si(exp.value, 100000) < 0: # TODO: test to find a good threshold
1478	x = mod_pow_int(self.ivalue, mpz_get_si(exp.value), self.__modulus.int32)
1479	else:
1480	mpz_init(x_mpz)
1481	_sig_on
1482	base = self.lift()
1483	mpz_powm(x_mpz, base.value, exp.value, self.__modulus.sageInteger.value)
1484	_sig_off
1485	x = mpz_get_si(x_mpz)
1486	mpz_clear(x_mpz)
1487	return self._new_c(x)
1488
1489
1490	def __invert__(IntegerMod_int self):
1491	"""
1492	Return the multiplicative inverse of self.
1493
1494	EXAMPLES:
1495	sage: ~mod(7,100)
1496	43
1497	"""
1498	if self.__modulus.inverses is not None:
1499	x = self.__modulus.inverses[self.ivalue]
1500	if x is None:
1501	raise ZeroDivisionError, "Inverse does not exist."
1502	else:
1503	return x
1504	else:
1505	return self._new_c(mod_inverse_int(self.ivalue, self.__modulus.int32))
1506
1507	def lift(IntegerMod_int self):
1508	"""
1509	Lift an integer modulo $n$ to the integers.
1510
1511	EXAMPLES:
1512	sage: a = Mod(8943, 2^10); type(a)
1513	<type 'sage.rings.integer_mod.IntegerMod_int'>
1514	sage: lift(a)
1515	751
1516	sage: a.lift()
1517	751
1518	"""
1519	cdef sage.rings.integer.Integer z
1520	z = sage.rings.integer.Integer()
1521	mpz_set_si(z.value, self.ivalue)
1522	return z
1523
1524	def __float__(IntegerMod_int self):
1525	return <double>self.ivalue
1526
1527	def __hash__(self):
1528	"""
1529	EXAMPLES:
1530	sage: a = Mod(89, 2^10)
1531	sage: hash(a)
1532	89
1533	"""
1534	return hash(self.ivalue)
1535
1536	cdef bint is_square_c(self) except -2:
1537	if self.ivalue <= 1:
1538	return 1
1539	moduli = self._parent.factored_order()
1540	cdef int val, e
1541	cdef int_fast32_t p
1542	if len(moduli) == 1:
1543	sage_p, e = moduli[0]
1544	p = sage_p
1545	if e == 1:
1546	return jacobi_int(self.ivalue, p) != -1
1547	elif p == 2:
1548	return self.pari().issquare() # TODO: implement directly
1549	elif self.ivalue % p == 0:
1550	val = self.lift().valuation(sage_p)
1551	return val >= e or (val % 2 == 0 and jacobi_int(self.ivalue / int(sage_p**val), p) != -1)
1552	else:
1553	return jacobi_int(self.ivalue, p) != -1
1554	else:
1555	for sage_p, e in moduli:
1556	p = sage_p
1557	if p == 2:
1558	if e > 1 and not self.pari().issquare(): # TODO: implement directly
1559	return 0
1560	elif e > 1 and self.ivalue % p == 0:
1561	val = self.lift().valuation(sage_p)
1562	if val < e and (val % 2 == 1 or jacobi_int(self.ivalue / int(sage_p**val), p) == -1):
1563	return 0
1564	elif jacobi_int(self.ivalue, p) == -1:
1565	return 0
1566	return 1
1567
1568	def sqrt(self, extend=True, all=False):
1569	cdef int_fast32_t i, n = self.__modulus.int32
1570	if n > 100:
1571	moduli = self._parent.factored_order()
1572	# Unless the modulus is tiny, test to see if we're in the really
1573	# easy case of n prime, n = 3 mod 4.
1574	if n > 100 and n % 4 == 3 and len(moduli) == 1 and moduli[0][1] == 1:
1575	if jacobi_int(self.ivalue, self.__modulus.int32) == 1:
1576	# it's a non-zero square, sqrt(a) = a^(p+1)/4
1577	i = mod_pow_int(self.ivalue, (self.__modulus.int32+1)/4, n)
1578	if i > n/2:
1579	i = n-i
1580	if all:
1581	return [self._new_c(i), self._new_c(n-i)]
1582	else:
1583	return self._new_c(i)
1584	elif self.ivalue == 0:
1585	return self
1586	elif not extend:
1587	raise ValueError, "self must be a square"
1588	# Now we use a heuristic to guess whether or not it will
1589	# be faster to just brute-force search for squares in a c loop...
1590	# TODO: more tuning?
1591	elif n <= 100 or n / (1 << len(moduli)) < 5000:
1592	if all:
1593	return [self._new_c(i) for i from 0 <= i < n if (i*i) % n == self.ivalue]
1594	else:
1595	for i from 0 <= i <= n/2:
1596	if (i*i) % n == self.ivalue:
1597	return self._new_c(i)
1598	if not extend:
1599	raise ValueError, "self must be a square"
1600	# Either it failed but extend was True, or the generic algorithm is better
1601	return IntegerMod_abstract.sqrt(self, extend=extend, all=all)
1602
1603
1604	def _balanced_abs(self):
1605	"""
1606	This function returns x or -x, whichever has a positive
1607	representative in -p/2 < x < p/2.
1608	"""
1609	if self.ivalue > self.__modulus.int32 / 2:
1610	return -self
1611	else:
1612	return self
1613
1614
1615
1616	### End of class
1617
1618
1619	cdef int_fast32_t gcd_int(int_fast32_t a, int_fast32_t b):
1620	"""
1621	Returns the gcd of a and b
1622	For use with IntegerMod_int
1623	AUTHOR:
1624	-- Robert Bradshaw
1625	"""
1626	cdef int_fast32_t tmp
1627	if a < b:
1628	tmp = b
1629	b = a
1630	a = tmp
1631	while b:
1632	tmp = b
1633	b = a % b
1634	a = tmp
1635	return a
1636
1637
1638	cdef int_fast32_t mod_inverse_int(int_fast32_t x, int_fast32_t n) except 0:
1639	"""
1640	Returns y such that xy=1 mod n
1641	For use in IntegerMod_int
1642	AUTHOR:
1643	-- Robert Bradshaw
1644	"""
1645	cdef int_fast32_t tmp, a, b, last_t, t, next_t, q
1646	a = n
1647	b = x
1648	t = 0
1649	next_t = 1
1650	while b:
1651	# a = s * n + t * x
1652	if b == 1:
1653	next_t = next_t % n
1654	if next_t < 0:
1655	next_t = next_t + n
1656	return next_t
1657	q = a / b
1658	tmp = b
1659	b = a % b
1660	a = tmp
1661	last_t = t
1662	t = next_t
1663	next_t = last_t - q * t
1664	raise ZeroDivisionError, "Inverse does not exist."
1665
1666
1667	cdef int_fast32_t mod_pow_int(int_fast32_t base, int_fast32_t exp, int_fast32_t n):
1668	"""
1669	Returns base^exp mod n
1670	For use in IntegerMod_int
1671	AUTHOR:
1672	-- Robert Bradshaw
1673	"""
1674	cdef int_fast32_t prod, pow2
1675	if exp <= 5:
1676	if exp == 0: return 1
1677	if exp == 1: return base
1678	prod = base * base % n
1679	if exp == 2: return prod
1680	if exp == 3: return (prod * base) % n
1681	if exp == 4: return (prod * prod) % n
1682
1683	pow2 = base
1684	if exp % 2: prod = base
1685	else: prod = 1
1686	exp = exp >> 1
1687	while(exp != 0):
1688	pow2 = pow2 * pow2
1689	if pow2 >= INTEGER_MOD_INT32_LIMIT: pow2 = pow2 % n
1690	if exp % 2:
1691	prod = prod * pow2
1692	if prod >= INTEGER_MOD_INT32_LIMIT: prod = prod % n
1693	exp = exp >> 1
1694
1695	if prod > n:
1696	prod = prod % n
1697	return prod
1698
1699
1700	cdef int jacobi_int(int_fast32_t a, int_fast32_t m) except -2:
1701	"""
1702	Calculates the jacobi symbol (a/n)
1703	For use in IntegerMod_int
1704	AUTHOR:
1705	-- Robert Bradshaw
1706	"""
1707	cdef int s, jacobi = 1
1708	cdef int_fast32_t b
1709
1710	a = a % m
1711
1712	while 1:
1713	if a == 0:
1714	return 0 # gcd was nontrivial
1715	elif a == 1:
1716	return jacobi
1717	s = 0
1718	while (1 << s) & a == 0:
1719	s += 1
1720	b = a >> s
1721	# Now a = 2^s * b
1722
1723	# factor out (2/m)^s term
1724	if s % 2 == 1 and (m % 8 == 3 or m % 8 == 5):
1725	jacobi = -jacobi
1726
1727	if b == 1:
1728	return jacobi
1729
1730	# quadratic reciprocity
1731	if b % 4 == 3 and m % 4 == 3:
1732	jacobi = -jacobi
1733	a = m % b
1734	m = b
1735
1736
1737	def test_gcd(a, b):
1738	return gcd_int(int(a), int(b))
1739
1740	def test_mod_inverse(a, b):
1741	return mod_inverse_int(int(a), int(b))
1742
1743
1744
1745	######################################################################
1746	# class IntegerMod_int64
1747	######################################################################
1748
1749	cdef class IntegerMod_int64(IntegerMod_abstract):
1750	"""
1751	Elements of $\Z/n\Z$ for n small enough to be operated on in 64 bits
1752	AUTHORS:
1753	-- Robert Bradshaw (2006-09-14)
1754	"""
1755
1756	def __init__(self, parent, value, empty=False):
1757	"""
1758	EXAMPLES:
1759	sage: a = Mod(10,3^10); a
1760	10
1761	sage: type(a)
1762	<type 'sage.rings.integer_mod.IntegerMod_int64'>
1763	sage: loads(a.dumps()) == a
1764	True
1765	sage: Mod(5, 2^31)
1766	5
1767	"""
1768	IntegerMod_abstract.__init__(self, parent)
1769	if empty:
1770	return
1771	cdef int_fast64_t x
1772	if PY_TYPE_CHECK(value, int):
1773	x = value
1774	self.ivalue = x % self.__modulus.int64
1775	if self.ivalue < 0:
1776	self.ivalue = self.ivalue + self.__modulus.int64
1777	return
1778	cdef sage.rings.integer.Integer z
1779	if isinstance(value, sage.rings.integer.Integer):
1780	z = value
1781	elif isinstance(value, rational.Rational):
1782	z = value % self.__modulus.sageInteger
1783	else:
1784	z = sage.rings.integer_ring.Z(value)
1785	self.set_from_mpz(z.value)
1786
1787	cdef IntegerMod_int64 _new_c(self, int_fast64_t value):
1788	cdef IntegerMod_int64 x
1789	x = PY_NEW(IntegerMod_int64)
1790	x.__modulus = self.__modulus
1791	x._parent = self._parent
1792	x.ivalue = value
1793	return x
1794
1795	cdef void set_from_mpz(self, mpz_t value):
1796	if mpz_sgn(value) == -1 or mpz_cmp_si(value, self.__modulus.int64) >= 0:
1797	self.ivalue = mpz_fdiv_ui(value, self.__modulus.int64)
1798	else:
1799	self.ivalue = mpz_get_si(value)
1800
1801	cdef void set_from_long(self, long value):
1802	self.ivalue = value % self.__modulus.int64
1803
1804	cdef void set_from_int(IntegerMod_int64 self, int_fast64_t ivalue):
1805	if ivalue < 0:
1806	self.ivalue = self.__modulus.int64 + (ivalue % self.__modulus.int64) # Is ivalue % self.__modulus.int64 actually negative?
1807	elif ivalue >= self.__modulus.int64:
1808	self.ivalue = ivalue % self.__modulus.int64
1809	else:
1810	self.ivalue = ivalue
1811
1812	cdef int_fast64_t get_int_value(IntegerMod_int64 self):
1813	return self.ivalue
1814
1815
1816	cdef int _cmp_c_impl(left, Element right) except -2:
1817	"""
1818	EXAMPLES:
1819	sage: mod(5,13^5) == mod(13^5+5,13^5)
1820	True
1821	sage: mod(5,13^5) == mod(8,13^5)
1822	False
1823	sage: mod(5,13^5) == mod(5,13)
1824	True
1825	sage: mod(0, 13^5) == 0
1826	True
1827	sage: mod(0, 13^5) == int(0)
1828	True
1829	"""
1830	if (<IntegerMod_int64>left).ivalue == (<IntegerMod_int64>right).ivalue: return 0
1831	elif (<IntegerMod_int64>left).ivalue < (<IntegerMod_int64>right).ivalue: return -1
1832	else: return 1
1833
1834	def __richcmp__(left, right, int op):
1835	return (<Element>left)._richcmp(right, op)
1836
1837
1838	def is_one(IntegerMod_int64 self):
1839	"""
1840	Returns \code{True} if this is $1$, otherwise \code{False}.
1841
1842	EXAMPLES:
1843	sage: (mod(-1,5^10)^2).is_one()
1844	True
1845	sage: mod(0,5^10).is_one()
1846	False
1847	"""
1848	return self.ivalue == 1
1849
1850	def __nonzero__(IntegerMod_int64 self):
1851	"""
1852	Returns \code{True} if this is not $0$, otherwise \code{False}.
1853
1854	EXAMPLES:
1855	sage: mod(13,5^10).is_zero()
1856	False
1857	sage: mod(5^12,5^10).is_zero()
1858	True
1859	"""
1860	return self.ivalue != 0
1861
1862	def is_unit(IntegerMod_int64 self):
1863	return gcd_int64(self.ivalue, self.__modulus.int64) == 1
1864
1865	def __crt(IntegerMod_int64 self, IntegerMod_int64 other):
1866	"""
1867	Use the Chinese Remainder Theorem to find an element of the
1868	integers modulo the product of the moduli that reduces to self
1869	and to other. The modulus of other must be coprime to the
1870	modulus of self.
1871	EXAMPLES:
1872	sage: a = mod(3,5^10)
1873	sage: b = mod(2,7)
1874	sage: a.crt(b)
1875	29296878
1876	sage: type(a.crt(b)) == type(b.crt(a)) and type(a.crt(b)) == type(mod(1, 7 * 5^10))
1877	True
1878
1879	sage: a = mod(3,10^10)
1880	sage: b = mod(2,9)
1881	sage: a.crt(b)
1882	80000000003
1883	sage: type(a.crt(b)) == type(b.crt(a)) and type(a.crt(b)) == type(mod(1, 9 * 10^10))
1884	True
1885
1886	AUTHOR:
1887	-- Robert Bradshaw
1888	"""
1889	cdef IntegerMod_int64 lift
1890	cdef int_fast64_t x
1891
1892	lift = IntegerMod_int64(integer_mod_ring.IntegerModRing(self.__modulus.int64 * other.__modulus.int64), None, empty=True)
1893
1894	try:
1895	x = (other.ivalue - self.ivalue % other.__modulus.int64) * mod_inverse_int64(self.__modulus.int64, other.__modulus.int64)
1896	lift.set_from_int( x * self.__modulus.int64 + self.ivalue )
1897	return lift
1898	except ZeroDivisionError:
1899	raise ZeroDivisionError, "moduli must be coprime"
1900
1901
1902	def __copy__(IntegerMod_int64 self):
1903	return self._new_c(self.ivalue)
1904
1905	cdef ModuleElement _add_c_impl(self, ModuleElement right):
1906	"""
1907	EXAMPLES:
1908	sage: R = Integers(10^5)
1909	sage: R(7) + R(8)
1910	15
1911	"""
1912	cdef int_fast64_t x
1913	x = self.ivalue + (<IntegerMod_int64>right).ivalue
1914	if x >= self.__modulus.int64:
1915	x = x - self.__modulus.int64
1916	return self._new_c(x)
1917
1918	cdef ModuleElement _sub_c_impl(self, ModuleElement right):
1919	"""
1920	EXAMPLES:
1921	sage: R = Integers(10^5)
1922	sage: R(7) - R(8)
1923	99999
1924	"""
1925	cdef int_fast64_t x
1926	x = self.ivalue - (<IntegerMod_int64>right).ivalue
1927	if x < 0:
1928	x = x + self.__modulus.int64
1929	return self._new_c(x)
1930
1931	cdef ModuleElement _neg_c_impl(self):
1932	"""
1933	EXAMPLES:
1934	sage: -mod(7,10^5)
1935	99993
1936	sage: -mod(0,10^6)
1937	0
1938	"""
1939	if self.ivalue == 0:
1940	return self
1941	return self._new_c(self.__modulus.int64 - self.ivalue)
1942
1943	cdef RingElement _mul_c_impl(self, RingElement right):
1944	"""
1945	EXAMPLES:
1946	sage: R = Integers(10^5)
1947	sage: R(700) * R(800)
1948	60000
1949	"""
1950	return self._new_c((self.ivalue * (<IntegerMod_int64>right).ivalue) % self.__modulus.int64)
1951
1952	cdef RingElement _div_c_impl(self, RingElement right):
1953	"""
1954	EXAMPLES:
1955	sage: R = Integers(10^5)
1956	sage: R(2)/3
1957	33334
1958	"""
1959	return self._new_c((self.ivalue * mod_inverse_int64((<IntegerMod_int64>right).ivalue,
1960	self.__modulus.int64) ) % self.__modulus.int64)
1961
1962	def __int__(IntegerMod_int64 self):
1963	return self.ivalue
1964
1965	def __long__(IntegerMod_int64 self):
1966	return self.ivalue
1967
1968	def __mod__(IntegerMod_int64 self, right):
1969	right = int(right)
1970	if self.__modulus.int64 % right != 0:
1971	raise ZeroDivisionError, "Error - reduction modulo right not defined."
1972	return integer_mod_ring.IntegerModRing(right)(self)
1973
1974	def __pow__(IntegerMod_int64 self, right, m): # NOTE: m ignored, always use modulus of parent ring
1975	"""
1976	EXAMPLES:
1977	sage: R = Integers(10)
1978	sage: R(2)^10
1979	4
1980	sage: p = next_prime(10^5)
1981	sage: R = Integers(p)
1982	sage: R(1234)^(p-1)
1983	1
1984	"""
1985	cdef sage.rings.integer.Integer exp, base
1986	exp = sage.rings.integer_ring.Z(right)
1987	cdef int_fast64_t x
1988	cdef mpz_t x_mpz
1989	if mpz_sgn(exp.value) >= 0 and mpz_cmp_si(exp.value, 100000) < 0: # TODO: test to find a good threshold
1990	x = mod_pow_int64(self.ivalue, mpz_get_si(exp.value), self.__modulus.int64)
1991	else:
1992	mpz_init(x_mpz)
1993	_sig_on
1994	base = self.lift()
1995	mpz_powm(x_mpz, base.value, exp.value, self.__modulus.sageInteger.value)
1996	_sig_off
1997	x = mpz_get_si(x_mpz)
1998	mpz_clear(x_mpz)
1999	return self._new_c(x)
2000
2001	def __lshift__(IntegerMod_int64 self, int right):
2002	r"""
2003	EXAMPLES:
2004	sage: e = Mod(5, 2^31 - 1)
2005	sage: e<<32
2006	10
2007	sage: e*2^32
2008	10
2009	"""
2010	return self._new_c((self.ivalue << right) % self.__modulus.int64)
2011
2012	def __rshift__(IntegerMod_int64 self, int right):
2013	"""
2014	Divide self by $2^{\text{right}}$ and take floor via bit shifting.
2015
2016	EXAMPLES:
2017	sage: e = Mod(8, 2^31 - 1)
2018	sage: e >> 3
2019	1
2020	sage: int(e)/int(2^3)
2021	1
2022	"""
2023	return self._new_c(self.ivalue >> right)
2024
2025	def __invert__(IntegerMod_int64 self):
2026	"""
2027	Return the multiplicative inverse of self.
2028
2029	EXAMPLES:
2030	sage: a = mod(7,2^40); type(a)
2031	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
2032	sage: ~a
2033	471219269047
2034	sage: a
2035	7
2036	"""
2037	return self._new_c(mod_inverse_int64(self.ivalue, self.__modulus.int64))
2038
2039	def lift(IntegerMod_int64 self):
2040	"""
2041	Lift an integer modulo $n$ to the integers.
2042
2043	EXAMPLES:
2044	sage: a = Mod(8943, 2^25); type(a)
2045	<type 'sage.rings.integer_mod.IntegerMod_int64'>
2046	sage: lift(a)
2047	8943
2048	sage: a.lift()
2049	8943
2050	"""
2051	cdef sage.rings.integer.Integer z
2052	z = sage.rings.integer.Integer()
2053	mpz_set_si(z.value, self.ivalue)
2054	return z
2055
2056	def __float__(IntegerMod_int64 self):
2057	"""
2058	Coerce self to a float.
2059
2060	EXAMPLES:
2061	sage: a = Mod(8943, 2^35)
2062	sage: float(a)
2063	8943.0
2064	"""
2065	return <double>self.ivalue
2066
2067	def __hash__(self):
2068	"""
2069	Compute hash of self.
2070
2071	This is a combination of the hash of the underlying integer
2072	and the modulus.
2073
2074	EXAMPLES:
2075	sage: a = Mod(8943, 2^35)
2076	sage: hash(a)
2077	8943
2078	"""
2079
2080	return hash(self.ivalue)
2081
2082	def _balanced_abs(self):
2083	"""
2084	This function returns x or -x, whichever has a positive
2085	representative in -p/2 < x < p/2.
2086	"""
2087	if self.ivalue > self.__modulus.int64 / 2:
2088	return -self
2089	else:
2090	return self
2091
2092
2093	### End of class
2094
2095
2096	cdef int_fast64_t gcd_int64(int_fast64_t a, int_fast64_t b):
2097	"""
2098	Returns the gcd of a and b
2099	For use with IntegerMod_int64
2100	AUTHOR:
2101	-- Robert Bradshaw
2102	"""
2103	cdef int_fast64_t tmp
2104	if a < b:
2105	tmp = b
2106	b = a
2107	a = tmp
2108	while b:
2109	tmp = b
2110	b = a % b
2111	a = tmp
2112	return a
2113
2114
2115	cdef int_fast64_t mod_inverse_int64(int_fast64_t x, int_fast64_t n) except 0:
2116	"""
2117	Returns y such that xy=1 mod n
2118	For use in IntegerMod_int64
2119	AUTHOR:
2120	-- Robert Bradshaw
2121	"""
2122	cdef int_fast64_t tmp, a, b, last_t, t, next_t, q
2123	a = n
2124	b = x
2125	t = 0
2126	next_t = 1
2127	while b:
2128	# a = s * n + t * x
2129	if b == 1:
2130	next_t = next_t % n
2131	if next_t < 0:
2132	next_t = next_t + n
2133	return next_t
2134	q = a / b
2135	tmp = b
2136	b = a % b
2137	a = tmp
2138	last_t = t
2139	t = next_t
2140	next_t = last_t - q * t
2141	raise ZeroDivisionError, "Inverse does not exist."
2142
2143
2144	cdef int_fast64_t mod_pow_int64(int_fast64_t base, int_fast64_t exp, int_fast64_t n):
2145	"""
2146	Returns base^exp mod n
2147	For use in IntegerMod_int64
2148	AUTHOR:
2149	-- Robert Bradshaw
2150	"""
2151	cdef int_fast64_t prod, pow2
2152	if exp <= 5:
2153	if exp == 0: return 1
2154	if exp == 1: return base
2155	prod = base * base % n
2156	if exp == 2: return prod
2157	if exp == 3: return (prod * base) % n
2158	if exp == 4: return (prod * prod) % n
2159
2160	pow2 = base
2161	if exp % 2: prod = base
2162	else: prod = 1
2163	exp = exp >> 1
2164	while(exp != 0):
2165	pow2 = pow2 * pow2
2166	if pow2 >= INTEGER_MOD_INT64_LIMIT: pow2 = pow2 % n
2167	if exp % 2:
2168	prod = prod * pow2
2169	if prod >= INTEGER_MOD_INT64_LIMIT: prod = prod % n
2170	exp = exp >> 1
2171
2172	if prod > n:
2173	prod = prod % n
2174	return prod
2175
2176
2177	cdef int jacobi_int64(int_fast64_t a, int_fast64_t m) except -2:
2178	"""
2179	Calculates the jacobi symbol (a/n)
2180	For use in IntegerMod_int64
2181	AUTHOR:
2182	-- Robert Bradshaw
2183	"""
2184	cdef int s, jacobi = 1
2185	cdef int_fast64_t b
2186
2187	a = a % m
2188
2189	while 1:
2190	if a == 0:
2191	return 0 # gcd was nontrivial
2192	elif a == 1:
2193	return jacobi
2194	s = 0
2195	while (1 << s) & a == 0:
2196	s += 1
2197	b = a >> s
2198	# Now a = 2^s * b
2199
2200	# factor out (2/m)^s term
2201	if s % 2 == 1 and (m % 8 == 3 or m % 8 == 5):
2202	jacobi = -jacobi
2203
2204	if b == 1:
2205	return jacobi
2206
2207	# quadratic reciprocity
2208	if b % 4 == 3 and m % 4 == 3:
2209	jacobi = -jacobi
2210	a = m % b
2211	m = b
2212
2213
2214	########################
2215	# Square root functions
2216	########################
2217
2218	def square_root_mod_prime_power(IntegerMod_abstract a, p, e):
2219	r"""
2220	Calculates the square root of $a$, where $a$ is an integer mod $p^e$.
2221
2222	ALGORITHM:
2223	Perform p-adically by stripping off even powers of $p$ to get
2224	a unit and lifting $\sqrt{unit} mod p$ via Newton's method.
2225
2226	AUTHOR:
2227	-- Robert Bradshaw
2228	"""
2229	if a.is_zero() or a.is_one():
2230	return a
2231
2232	if p == 2:
2233	if e == 1:
2234	return a
2235	# TODO: implement something that isn't totally idiotic.
2236	for x in a.parent():
2237	if x**2 == a:
2238	return x
2239
2240	# strip off even powers of p
2241	cdef int i, val = a.lift().valuation(p)
2242	if val % 2 == 1:
2243	raise ValueError, "self must be a square."
2244	if val > 0:
2245	unit = a._parent(a.lift() // p**val)
2246	else:
2247	unit = a
2248
2249	# find square root of unit mod p
2250	x = unit.parent()(square_root_mod_prime(mod(unit, p), p))
2251
2252	# lift p-adically using Newton iteration
2253	# this is done to higher precision than neccesary except at the last step
2254	one_half = ~(a._new_c_from_long(2))
2255	for i from 0 <= i < ceil(log(e)/log(2)) - val/2:
2256	x = (x+unit/x) * one_half
2257
2258	# multiply in powers of p (if any)
2259	if val > 0:
2260	x = p*(val // 2)
2261	return x
2262
2263	def square_root_mod_prime(IntegerMod_abstract a, p=None):
2264	r"""
2265	Calculates the square root of a, where a is an integer mod p.
2266
2267	ALGORITHM:
2268	Several cases based on residue class of p mod 16.
2269
2270	\begin{verbatim}
2271	$p$ mod 2 = 0 $\Rightarrow$ p = 2 so \sqrt{a} = a$.
2272	$p$ mod 4 = 3 $\Rightarrow \sqrt{a} = a^{(p+1)/4}$.
2273	$p$ mod 8 = 5 $\Rightarrow \sqrt{a} = \zeta i a$ where $\zeta = (2a)^{(p-5)/8}$, $i=\sqrt{-1}$.
2274	$p$ mod 16 = 9$ Similar, work in a bi-quadratic extension of $\FF_p$.
2275	$p$ mod 16 = 1$ Variant of Cipolla-Lehmer, using Lucas functions.
2276	\end{verbatim}
2277
2278	REFERENCES:
2279	Siguna M\:uller. 'On the Computation of Square Roots in Finite Fields'
2280	Designs, Codes and Cryptography, Volume 31, Issue 3 (March 2004)
2281
2282	A. Oliver L. Atkin. 'Probabilistic primality testing' (Section 4)
2283	In P. Flajolet and P. Zimmermann, editors, Analysis of Algorithms Seminar I. INRIA Research Report XXX, 1992.
2284	Summary by F. Morain.
2285	\url{http://citeseer.ist.psu.edu/atkin92probabilistic.html}
2286
2287	H. Postl. 'Fast evaluation of Dickson Polynomials'
2288	Contrib. to General Algebra, Vol. 6 (1988) pp. 223--225
2289
2290	AUTHOR:
2291	Robert Bradshaw
2292
2293	TESTS:
2294	Every case appears in the first hundred primes.
2295	sage: from sage.rings.integer_mod import square_root_mod_prime # sqrt() uses brute force for small p
2296	sage: all([square_root_mod_prime(aa)^2 == aa for p in prime_range(100) for a in Integers(p)])
2297	True
2298	"""
2299	if a.is_zero() or a.is_one():
2300	return a
2301
2302	if p is None:
2303	p = a.parent().order()
2304	if p < PyInt_GetMax():
2305	p = int(p)
2306
2307	cdef int p_mod_16 = p % 16
2308
2309	if p_mod_16 % 2 == 0: # p == 2
2310	return a
2311
2312	elif p_mod_16 % 4 == 3:
2313	return a ** ((p+1)//4)
2314
2315	elif p_mod_16 % 8 == 5:
2316	two_a = a+a
2317	zeta = two_a ** ((p-5)//8)
2318	i = two_a ** ((p-1)//4)
2319	return zetaa(i-1)
2320
2321	elif p_mod_16 == 9:
2322	s = (a+a) ** ((p-1)//4)
2323	if s.is_one():
2324	d = a._parent.quadratic_nonresidue()
2325	d2 = d*d
2326	z = (2 * d2 * a) ** ((p-9)//16)
2327	i = 2 * d2 * zz a
2328	return zda*(i-1)
2329	else:
2330	z = (a+a) ** ((p-9)//16)
2331	i = 2 * zz a
2332	return za(i-1)
2333
2334	else: # p_mod_16 == 1
2335
2336	four = a._new_c_from_long(4)
2337
2338	if a == four:
2339	return a._new_c_from_long(2)
2340
2341	if not (<IntegerMod_abstract>(a - four)).is_square_c():
2342	t = 1
2343	P = a - 2
2344	return fast_lucas((p-1) // 4, P)
2345
2346	else:
2347	t = a._new_c_from_long(2)
2348	while (<IntegerMod_abstract>(att - four)).is_square_c():
2349	t += 1
2350	P = att - 2
2351	return fast_lucas((p-1)//4, P)/t
2352
2353
2354	def fast_lucas(mm, IntegerMod_abstract P):
2355	"""
2356	Return $V_k(P, 1)$ where $V_k$ is the Lucas function
2357	defined by the recursive relation
2358
2359	$V_k(P, Q) = PV_{k-1}(P, Q) - QV_{k-2}(P, Q)$
2360
2361	with $V_0 = 2, V_1(P_Q) = P$.
2362
2363	REFERENCES:
2364	H. Postl. 'Fast evaluation of Dickson Polynomials'
2365	Contrib. to General Algebra, Vol. 6 (1988) pp. 223\202\304\354225
2366
2367	AUTHOR:
2368	Robert Bradshaw
2369
2370	TESTS:
2371	sage: from sage.rings.integer_mod import fast_lucas, slow_lucas
2372	sage: all([fast_lucas(k, a) == slow_lucas(k, a) for a in Integers(23) for k in range(13)])
2373	True
2374	"""
2375	if mm == 0:
2376	return 2
2377	elif mm == 1:
2378	return P
2379
2380	cdef sage.rings.integer.Integer m
2381	m = <sage.rings.integer.Integer>mm if PY_TYPE_CHECK(mm, sage.rings.integer.Integer) else sage.rings.integer.Integer(mm)
2382	two = P._new_c_from_long(2)
2383	d1 = P
2384	d2 = P*P - two
2385
2386	_sig_on
2387	cdef int j
2388	for j from mpz_sizeinbase(m.value, 2)-1 > j > 0:
2389	if mpz_tstbit(m.value, j):
2390	d1 = d1*d2 - P
2391	d2 = d2*d2 - two
2392	else:
2393	d2 = d1*d2 - P
2394	d1 = d1*d1 - two
2395	_sig_off
2396	if mpz_odd_p(m.value):
2397	return d1*d2 - P
2398	else:
2399	return d1*d1 - two
2400
2401	def slow_lucas(k, P, Q=1):
2402	"""
2403	Lucas function defined using the definition, for consitancy testing.
2404	"""
2405	if k == 0:
2406	return 2
2407	elif k == 1:
2408	return P
2409	else:
2410	return Pslow_lucas(k-1, P, Q) - Qslow_lucas(k-2, P, Q)
2411
2412
2413	############# Homomorphisms ###############
2414
2415	cdef class IntegerMod_hom(Morphism):
2416	cdef IntegerMod_abstract zero
2417	cdef NativeIntStruct modulus
2418	def __init__(self, parent):
2419	Morphism.__init__(self, parent)
2420	self.zero = self._codomain(0)
2421	self.modulus = self._codomain._pyx_order
2422	cdef Element _call_c(self, Element x):
2423	return IntegerMod(self.codomain, x)
2424
2425	cdef class IntegerMod_to_IntegerMod(IntegerMod_hom):
2426	"""
2427	Very fast IntegerMod to IntegerMod homomorphism.
2428
2429	sage: from sage.rings.integer_mod import IntegerMod_to_IntegerMod
2430	sage: R = Integers(2^8); S = Integers(2^3)
2431	sage: hom = IntegerMod_to_IntegerMod(HomsetWithBase(R, S))
2432	sage: hom(R(-1))
2433	7
2434	sage: R = Integers(2^18); S = Integers(2^3)
2435	sage: type(R(1))
2436	<type 'sage.rings.integer_mod.IntegerMod_int64'>
2437	sage: hom = IntegerMod_to_IntegerMod(HomsetWithBase(R, S))
2438	sage: hom(R(-1))
2439	7
2440	sage: R = Integers(2^118); S = Integers(2^3)
2441	sage: type(R(1))
2442	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
2443	sage: hom = IntegerMod_to_IntegerMod(HomsetWithBase(R, S))
2444	sage: hom(R(-1))
2445	7
2446	sage: R = Integers(3^118); S = Integers(3^109)
2447	sage: type(S(1))
2448	<type 'sage.rings.integer_mod.IntegerMod_gmp'>
2449	sage: hom = IntegerMod_to_IntegerMod(HomsetWithBase(R, S))
2450	sage: hom(R(-1)) + 1
2451	0
2452	"""
2453	cdef Element _call_c(self, Element x):
2454	cdef IntegerMod_abstract a
2455	if PY_TYPE_CHECK(x, IntegerMod_int):
2456	return (<IntegerMod_int>self.zero)._new_c((<IntegerMod_int>x).ivalue % self.modulus.int32)
2457	elif PY_TYPE_CHECK(x, IntegerMod_int64):
2458	return self.zero._new_c_from_long((<IntegerMod_int64>x).ivalue % self.modulus.int64)
2459	else: # PY_TYPE_CHECK(x, IntegerMod_gmp)
2460	a = self.zero._new_c_from_long(0)
2461	a.set_from_mpz((<IntegerMod_gmp>x).value)
2462	return a
2463
2464	cdef class Integer_to_IntegerMod(IntegerMod_hom):
2465	"""
2466	Fast $\Z \rightarrow \Z/n\Z$ morphism.
2467
2468	sage: from sage.rings.integer_mod import Integer_to_IntegerMod
2469	sage: R = Integers(5)
2470	sage: hom = Integer_to_IntegerMod(HomsetWithBase(ZZ, R))
2471	sage: hom(-1)
2472	4
2473	sage: R = Integers(2^20)
2474	sage: hom = Integer_to_IntegerMod(HomsetWithBase(ZZ, R))
2475	sage: hom(5)
2476	5
2477	sage: R = Integers(5^120)
2478	sage: hom = Integer_to_IntegerMod(HomsetWithBase(ZZ, R))
2479	sage: hom(2)
2480	2
2481	"""
2482	cdef Element _call_c(self, Element x):
2483	cdef IntegerMod_abstract a
2484	cdef Py_ssize_t res
2485	if self.modulus.table is not None:
2486	res = x % self.modulus.int64
2487	a = self.modulus.lookup(res)
2488	if a._parent is not self._codomain:
2489	a._parent = self._codomain
2490	# print (<Element>a)._parent, " is not ", parent
2491	return a
2492	else:
2493	a = self.zero._new_c_from_long(0)
2494	a.set_from_mpz((<Integer>x).value)
2495	return a
2496
2497	cdef class Int_to_IntegerMode(IntegerMod_hom):
2498	cdef Element _call_c(self, Element x):
2499	cdef IntegerMod_abstract a
2500	cdef Py_ssize_t res
2501	if self.modulus.table is not None:
2502	res = x
2503	res %= self.modulus.int64
2504	a = self.modulus.lookup(res)
2505	if a._parent is not self._codomain:
2506	a._parent = self._codomain
2507	# print (<Element>a)._parent, " is not ", parent
2508	return a
2509	else:
2510	return self.zero._new_c_from_long(x)

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

Download in other formats: