Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/rings/integer.pyx
4036 views
1
r"""
2
Elements of the ring `\ZZ` of integers
3
4
AUTHORS:
5
6
- William Stein (2005): initial version
7
8
- Gonzalo Tornaria (2006-03-02): vastly improved python/GMP
9
conversion; hashing
10
11
- Didier Deshommes (2006-03-06): numerous examples
12
and docstrings
13
14
- William Stein (2006-03-31): changes to reflect GMP bug fixes
15
16
- William Stein (2006-04-14): added GMP factorial method (since it's
17
now very fast).
18
19
- David Harvey (2006-09-15): added nth_root, exact_log
20
21
- David Harvey (2006-09-16): attempt to optimise Integer constructor
22
23
- Rishikesh (2007-02-25): changed quo_rem so that the rem is positive
24
25
- David Harvey, Martin Albrecht, Robert Bradshaw (2007-03-01):
26
optimized Integer constructor and pool
27
28
- Pablo De Napoli (2007-04-01): multiplicative_order should return
29
+infinity for non zero numbers
30
31
- Robert Bradshaw (2007-04-12): is_perfect_power, Jacobi symbol (with
32
Kronecker extension). Convert some methods to use GMP directly
33
rather than PARI, Integer(), PY_NEW(Integer)
34
35
- David Roe (2007-03-21): sped up valuation and is_square, added
36
val_unit, is_power, is_power_of and divide_knowing_divisible_by
37
38
- Robert Bradshaw (2008-03-26): gamma function, multifactorials
39
40
- Robert Bradshaw (2008-10-02): bounded squarefree part
41
42
- David Loeffler (2011-01-15): fixed bug #10625 (inverse_mod should accept an ideal as argument)
43
44
- Vincent Delecroix (2010-12-28): added unicode in Integer.__init__
45
46
EXAMPLES:
47
48
Add 2 integers::
49
50
sage: a = Integer(3) ; b = Integer(4)
51
sage: a + b == 7
52
True
53
54
Add an integer and a real number::
55
56
sage: a + 4.0
57
7.00000000000000
58
59
Add an integer and a rational number::
60
61
sage: a + Rational(2)/5
62
17/5
63
64
Add an integer and a complex number::
65
66
sage: b = ComplexField().0 + 1.5
67
sage: loads((a+b).dumps()) == a+b
68
True
69
70
sage: z = 32
71
sage: -z
72
-32
73
sage: z = 0; -z
74
0
75
sage: z = -0; -z
76
0
77
sage: z = -1; -z
78
1
79
80
Multiplication::
81
82
sage: a = Integer(3) ; b = Integer(4)
83
sage: a * b == 12
84
True
85
sage: loads((a * 4.0).dumps()) == a*b
86
True
87
sage: a * Rational(2)/5
88
6/5
89
90
::
91
92
sage: list([2,3]) * 4
93
[2, 3, 2, 3, 2, 3, 2, 3]
94
95
::
96
97
sage: 'sage'*Integer(3)
98
'sagesagesage'
99
100
COERCIONS:
101
102
Returns version of this integer in the multi-precision floating
103
real field R::
104
105
sage: n = 9390823
106
sage: RR = RealField(200)
107
sage: RR(n)
108
9.3908230000000000000000000000000000000000000000000000000000e6
109
110
"""
111
#*****************************************************************************
112
# Copyright (C) 2004,2006 William Stein <[email protected]>
113
# Copyright (C) 2006 Gonzalo Tornaria <[email protected]>
114
# Copyright (C) 2006 Didier Deshommes <[email protected]>
115
# Copyright (C) 2007 David Harvey <[email protected]>
116
# Copyright (C) 2007 Martin Albrecht <[email protected]>
117
# Copyright (C) 2007,2008 Robert Bradshaw <[email protected]>
118
# Copyright (C) 2007 David Roe <[email protected]>
119
#
120
# Distributed under the terms of the GNU General Public License (GPL)
121
# as published by the Free Software Foundation; either version 2 of
122
# the License, or (at your option) any later version.
123
# http://www.gnu.org/licenses/
124
#*****************************************************************************
125
126
doc="""
127
Integers
128
"""
129
130
import cython
131
132
import operator
133
134
import sys
135
136
include "../ext/gmp.pxi"
137
include "../ext/interrupt.pxi" # ctrl-c interrupt block support
138
include "../ext/stdsage.pxi"
139
include "../ext/python_list.pxi"
140
include "../ext/python_number.pxi"
141
include "../ext/python_int.pxi"
142
include "../structure/coerce.pxi" # for parent_c
143
include "../libs/pari/decl.pxi"
144
145
cdef extern from "limits.h":
146
long LONG_MAX
147
148
cdef extern from "math.h":
149
double sqrt_double "sqrt"(double)
150
cdef double log_c "log" (double)
151
cdef double ceil_c "ceil" (double)
152
int isnan(double)
153
154
cdef extern from "mpz_pylong.h":
155
cdef mpz_get_pylong(mpz_t src)
156
cdef mpz_get_pyintlong(mpz_t src)
157
cdef int mpz_set_pylong(mpz_t dst, src) except -1
158
cdef long mpz_pythonhash(mpz_t src)
159
160
cdef extern from "mpz_longlong.h":
161
cdef void mpz_set_longlong(mpz_t dst, long long src)
162
cdef void mpz_set_ulonglong(mpz_t dst, unsigned long long src)
163
164
cdef extern from "convert.h":
165
cdef void t_INT_to_ZZ( mpz_t value, long *g )
166
167
from sage.libs.pari.gen cimport gen as pari_gen, PariInstance
168
from sage.libs.flint.long_extras cimport *
169
170
import sage.rings.infinity
171
import sage.libs.pari.all
172
173
from sage.structure.element import canonical_coercion
174
175
cdef object numpy_long_interface = {'typestr': '=i4' if sizeof(long) == 4 else '=i8' }
176
cdef object numpy_int64_interface = {'typestr': '=i8'}
177
cdef object numpy_object_interface = {'typestr': '|O'}
178
179
cdef mpz_t mpz_tmp
180
mpz_init(mpz_tmp)
181
182
cdef public int set_mpz(Integer self, mpz_t value):
183
mpz_set(self.value, value)
184
185
cdef set_from_Integer(Integer self, Integer other):
186
mpz_set(self.value, other.value)
187
188
cdef set_from_int(Integer self, int other):
189
mpz_set_si(self.value, other)
190
191
cdef public mpz_t* get_value(Integer self):
192
return &self.value
193
194
195
def _test_mpz_set_longlong(long long v):
196
"""
197
TESTS::
198
199
sage: from sage.rings.integer import _test_mpz_set_longlong
200
sage: _test_mpz_set_longlong(1)
201
1
202
sage: _test_mpz_set_longlong(-1)
203
-1
204
sage: _test_mpz_set_longlong(100)
205
100
206
sage: _test_mpz_set_longlong(100000000000)
207
100000000000
208
sage: _test_mpz_set_longlong(2^32) == 2^32
209
True
210
sage: _test_mpz_set_longlong(-2^32) == -2^32
211
True
212
sage: all([_test_mpz_set_longlong(2^n) == 2^n for n in range(63)])
213
True
214
"""
215
cdef Integer z = PY_NEW(Integer)
216
mpz_set_longlong(z.value, v)
217
return z
218
219
cdef _digits_naive(mpz_t v,l,int offset,Integer base,digits):
220
"""
221
This method fills in digit entries in the list, l, using the most
222
basic digit algorithm -- repeat division by base.
223
224
INPUT:
225
226
- ``v`` - the value whose digits we want to put into the list
227
228
- ``l`` - the list to file
229
230
- ``offset`` - offset from the beginning of the list that we want
231
to fill at
232
233
- ``base`` -- the base to which we finding digits
234
235
- ``digits`` - a python sequence type with objects to use for digits
236
note that python negative index semantics are relied upon
237
238
AUTHORS:
239
240
- Joel B. Mohler (2009-01-16)
241
"""
242
cdef mpz_t mpz_value
243
cdef mpz_t mpz_res # used on one side of the 'if'
244
cdef Integer z # used on the other side of the 'if'
245
246
mpz_init(mpz_value)
247
mpz_set(mpz_value, v)
248
249
# we aim to avoid sage Integer creation if possible
250
if digits is None:
251
while mpz_cmp_si(mpz_value,0):
252
z = PY_NEW(Integer)
253
mpz_tdiv_qr(mpz_value, z.value, mpz_value, base.value)
254
l[offset] = z
255
offset += 1
256
else:
257
mpz_init(mpz_res)
258
while mpz_cmp_si(mpz_value,0):
259
mpz_tdiv_qr(mpz_value, mpz_res, mpz_value, base.value)
260
l[offset] = digits[mpz_get_si(mpz_res)]
261
offset += 1
262
mpz_clear(mpz_res)
263
264
mpz_clear(mpz_value)
265
266
cdef _digits_internal(mpz_t v,l,int offset,int power_index,power_list,digits):
267
"""
268
INPUT:
269
270
- ``v`` - the value whose digits we want to put into the list
271
272
- ``l`` - the list to file
273
274
- ``offset`` - offset from the beginning of the list that we want
275
to fill at
276
277
- ``power_index`` - a measure of size to fill and index to
278
power_list we're filling 1 << (power_index+1) digits
279
280
- ``power_list`` - a list of powers of the base, precomputed in
281
method digits digits - a python sequence type with objects to
282
use for digits note that python negative index semantics are
283
relied upon
284
285
AUTHORS:
286
287
- Joel B. Mohler (2008-03-13)
288
"""
289
cdef mpz_t mpz_res
290
cdef mpz_t mpz_quot
291
cdef Integer temp
292
cdef int v_int
293
if power_index < 5:
294
# It turns out that simple repeated division is very fast for
295
# relatively few digits. I don't think this is a real algorithmic
296
# statement, it's an annoyance introduced by memory allocation.
297
# I think that manual memory management with mpn_* would make the
298
# divide & conquer approach even faster, but the code would be much
299
# more complicated.
300
_digits_naive(v,l,offset,power_list[0],digits)
301
else:
302
mpz_init(mpz_quot)
303
mpz_init(mpz_res)
304
temp = power_list[power_index]
305
mpz_tdiv_qr(mpz_quot, mpz_res, v, temp.value)
306
if mpz_sgn(mpz_res) != 0:
307
_digits_internal(mpz_res,l,offset,power_index-1,power_list,digits)
308
if mpz_sgn(mpz_quot) != 0:
309
_digits_internal(mpz_quot,l,offset+(1<<power_index),power_index-1,power_list,digits)
310
mpz_clear(mpz_quot)
311
mpz_clear(mpz_res)
312
313
arith = None
314
cdef void late_import():
315
global arith
316
if arith is None:
317
import sage.rings.arith
318
arith = sage.rings.arith
319
320
MAX_UNSIGNED_LONG = 2 * sys.maxint
321
322
from sage.structure.sage_object cimport SageObject
323
from sage.structure.element cimport EuclideanDomainElement, ModuleElement, Element
324
from sage.structure.element import bin_op
325
from sage.structure.coerce_exceptions import CoercionException
326
327
import integer_ring
328
the_integer_ring = integer_ring.ZZ
329
330
initialized = False
331
cdef set_zero_one_elements():
332
global the_integer_ring, initialized
333
if initialized: return
334
the_integer_ring._zero_element = Integer(0)
335
the_integer_ring._one_element = Integer(1)
336
init_mpz_globals()
337
initialized = True
338
set_zero_one_elements()
339
340
cdef Integer zero = the_integer_ring._zero_element
341
cdef Integer one = the_integer_ring._one_element
342
343
# The documentation for the ispseudoprime() function in the PARI
344
# manual states that its result is always prime up to this 10^13.
345
cdef mpz_t PARI_PSEUDOPRIME_LIMIT
346
mpz_init(PARI_PSEUDOPRIME_LIMIT)
347
mpz_set_str(PARI_PSEUDOPRIME_LIMIT, "10000000000000", 10)
348
349
def is_Integer(x):
350
"""
351
Return true if x is of the Sage integer type.
352
353
EXAMPLES::
354
355
sage: from sage.rings.integer import is_Integer
356
sage: is_Integer(2)
357
True
358
sage: is_Integer(2/1)
359
False
360
sage: is_Integer(int(2))
361
False
362
sage: is_Integer(long(2))
363
False
364
sage: is_Integer('5')
365
False
366
"""
367
return PY_TYPE_CHECK(x, Integer)
368
369
cdef inline Integer as_Integer(x):
370
if PY_TYPE_CHECK(x, Integer):
371
return <Integer>x
372
else:
373
return Integer(x)
374
375
cdef class IntegerWrapper(Integer):
376
r"""
377
Rationale for the ``IntegerWrapper`` class:
378
379
With ``Integers``, the allocation/deallocation function slots are
380
hijacked with custom functions that stick already allocated
381
``Integers`` (with initialized ``parent`` and ``mpz_t`` fields)
382
into a pool on "deallocation" and then pull them out whenever a
383
new one is needed. Because ``Integers`` are so common, this is
384
actually a significant savings. However , this does cause issues
385
with subclassing a Python class directly from ``Integer`` (but
386
that's ok for a Cython class).
387
388
As a workaround, one can instead derive a class from the
389
intermediate class ``IntegerWrapper``, which sets statically its
390
alloc/dealloc methods to the *original* ``Integer`` alloc/dealloc
391
methods, before they are swapped manually for the custom ones.
392
393
The constructor of ``IntegerWrapper`` further allows for
394
specifying an alternative parent to ``IntegerRing()``.
395
"""
396
397
def __init__(self, x=None, unsigned int base=0, parent=None):
398
"""
399
We illustrate how to create integers with parents different
400
from ``IntegerRing()``::
401
402
sage: from sage.rings.integer import IntegerWrapper
403
404
sage: n = IntegerWrapper(3, parent=Primes()) # indirect doctest
405
sage: n
406
3
407
sage: n.parent()
408
Set of all prime numbers: 2, 3, 5, 7, ...
409
410
Pickling seems to work now (as of #10314)
411
412
sage: nn = loads(dumps(n))
413
sage: nn
414
3
415
sage: nn.parent()
416
Integer Ring
417
418
sage: TestSuite(n).run()
419
"""
420
if parent is not None:
421
Element.__init__(self, parent=parent)
422
Integer.__init__(self, x, base=base)
423
424
cdef class Integer(sage.structure.element.EuclideanDomainElement):
425
r"""
426
The ``Integer`` class represents arbitrary precision
427
integers. It derives from the ``Element`` class, so
428
integers can be used as ring elements anywhere in Sage.
429
430
Integer() interprets numbers and strings that begin with 0 as octal
431
numbers, and numbers and strings that begin with 0x as hexadecimal
432
numbers.
433
434
The class ``Integer`` is implemented in Cython, as a
435
wrapper of the GMP ``mpz_t`` integer type.
436
437
EXAMPLES::
438
439
sage: Integer(010)
440
8
441
sage: Integer(0x10)
442
16
443
sage: Integer(10)
444
10
445
sage: Integer('0x12')
446
18
447
sage: Integer('012')
448
10
449
450
Conversion from PARI::
451
452
sage: Integer(pari('-10380104371593008048799446356441519384'))
453
-10380104371593008048799446356441519384
454
sage: Integer(pari('Pol([-3])'))
455
-3
456
"""
457
458
def __cinit__(self):
459
mpz_init(self.value)
460
self._parent = <SageObject>the_integer_ring
461
462
def __pyxdoc__init__(self):
463
"""
464
You can create an integer from an int, long, string literal, or
465
integer modulo N.
466
467
EXAMPLES::
468
469
sage: Integer(495)
470
495
471
sage: Integer('495949209809328523')
472
495949209809328523
473
sage: Integer(Mod(3,7))
474
3
475
sage: 2^3
476
8
477
"""
478
479
def __init__(self, x=None, base=0):
480
"""
481
EXAMPLES::
482
483
sage: a = long(-901824309821093821093812093810928309183091832091)
484
sage: b = ZZ(a); b
485
-901824309821093821093812093810928309183091832091
486
sage: ZZ(b)
487
-901824309821093821093812093810928309183091832091
488
sage: ZZ('-901824309821093821093812093810928309183091832091')
489
-901824309821093821093812093810928309183091832091
490
sage: ZZ(int(-93820984323))
491
-93820984323
492
sage: ZZ(ZZ(-901824309821093821093812093810928309183091832091))
493
-901824309821093821093812093810928309183091832091
494
sage: ZZ(QQ(-901824309821093821093812093810928309183091832091))
495
-901824309821093821093812093810928309183091832091
496
sage: ZZ(RR(2.0)^80)
497
1208925819614629174706176
498
sage: ZZ(QQbar(sqrt(28-10*sqrt(3)) + sqrt(3)))
499
5
500
sage: ZZ(AA(32).nth_root(5))
501
2
502
sage: ZZ(pari('Mod(-3,7)'))
503
4
504
sage: ZZ('sage')
505
Traceback (most recent call last):
506
...
507
TypeError: unable to convert x (=sage) to an integer
508
sage: Integer('zz',36).str(36)
509
'zz'
510
sage: ZZ('0x3b').str(16)
511
'3b'
512
sage: ZZ( ZZ(5).digits(3) , 3)
513
5
514
sage: import numpy
515
sage: ZZ(numpy.int64(7^7))
516
823543
517
sage: ZZ(numpy.ubyte(-7))
518
249
519
sage: ZZ(True)
520
1
521
sage: ZZ(False)
522
0
523
sage: ZZ(1==0)
524
0
525
sage: ZZ('+10')
526
10
527
528
::
529
530
sage: k = GF(2)
531
sage: ZZ( (k(0),k(1)), 2)
532
2
533
534
::
535
536
sage: ZZ(float(2.0))
537
2
538
sage: ZZ(float(1.0/0.0))
539
Traceback (most recent call last):
540
...
541
OverflowError: cannot convert float infinity to integer
542
sage: ZZ(float(0.0/0.0))
543
Traceback (most recent call last):
544
...
545
ValueError: cannot convert float NaN to integer
546
547
::
548
549
sage: class MyInt(int):
550
... pass
551
sage: class MyLong(long):
552
... pass
553
sage: class MyFloat(float):
554
... pass
555
sage: ZZ(MyInt(3))
556
3
557
sage: ZZ(MyLong(4))
558
4
559
sage: ZZ(MyFloat(5))
560
5
561
562
::
563
564
sage: Integer(u'0')
565
0
566
sage: Integer(u'0X2AEEF')
567
175855
568
569
Test conversion from PARI (#11685)::
570
571
sage: ZZ(pari(-3))
572
-3
573
sage: ZZ(pari("-3.0"))
574
-3
575
sage: ZZ(pari("-3.5"))
576
Traceback (most recent call last):
577
...
578
TypeError: Attempt to coerce non-integral real number to an Integer
579
sage: ZZ(pari("1e100"))
580
Traceback (most recent call last):
581
...
582
PariError: precision too low (10)
583
sage: ZZ(pari("10^50"))
584
100000000000000000000000000000000000000000000000000
585
sage: ZZ(pari("Pol(3)"))
586
3
587
sage: ZZ(GF(3^20,'t')(1))
588
1
589
sage: ZZ(pari(GF(3^20,'t')(1)))
590
1
591
sage: x = polygen(QQ)
592
sage: K.<a> = NumberField(x^2+3)
593
sage: ZZ(a^2)
594
-3
595
sage: ZZ(pari(a)^2)
596
-3
597
sage: ZZ(pari("Mod(x, x^3+x+1)")) # Note error message refers to lifted element
598
Traceback (most recent call last):
599
...
600
TypeError: Unable to coerce PARI x to an Integer
601
602
Test coercion of p-adic with negative valuation::
603
604
sage: ZZ(pari(Qp(11)(11^-7)))
605
Traceback (most recent call last):
606
...
607
TypeError: Cannot convert p-adic with negative valuation to an integer
608
609
Test converting a list with a very large base::
610
611
sage: a=ZZ(randint(0,2^128-1))
612
sage: L = a.digits(2^64)
613
sage: a == sum([x * 2^(64*i) for i,x in enumerate(L)])
614
True
615
sage: a == ZZ(L,base=2^64)
616
True
617
"""
618
# TODO: All the code below should somehow be in an external
619
# cdef'd function. Then e.g., if a matrix or vector or
620
# polynomial is getting filled by mpz_t's, it can use the
621
# rules below to do the fill construction of mpz_t's, but
622
# without the overhead of creating any Python objects at all.
623
# The cdef's function should be of the form
624
# mpz_init_set_sage(mpz_t y, object x)
625
# Then this function becomes the one liner:
626
# mpz_init_set_sage(self.value, x)
627
628
cdef Integer tmp
629
cdef char* xs
630
cdef int paritype
631
cdef Py_ssize_t j
632
cdef object otmp
633
cdef unsigned int ibase
634
635
cdef Element lift
636
637
if x is None:
638
if mpz_sgn(self.value) != 0:
639
mpz_set_si(self.value, 0)
640
641
else:
642
# First do all the type-check versions; these are fast.
643
644
if PY_TYPE_CHECK(x, Integer):
645
set_from_Integer(self, <Integer>x)
646
647
elif PY_TYPE_CHECK(x, bool):
648
mpz_set_si(self.value, PyInt_AS_LONG(x))
649
650
elif PyInt_Check(x):
651
mpz_set_si(self.value, PyInt_AS_LONG(x))
652
653
elif PyLong_Check(x):
654
mpz_set_pylong(self.value, x)
655
656
elif PyFloat_Check(x):
657
n = long(x)
658
if n == x:
659
mpz_set_pylong(self.value, n)
660
else:
661
raise TypeError, "Cannot convert non-integral float to integer"
662
663
elif PY_TYPE_CHECK(x, pari_gen):
664
# Simplify and lift until we get an integer
665
while typ((<pari_gen>x).g) != t_INT:
666
x = x.simplify()
667
paritype = typ((<pari_gen>x).g)
668
if paritype == t_INT:
669
break
670
elif paritype == t_REAL:
671
# Check that the fractional part is zero
672
if not x.frac().gequal0():
673
raise TypeError, "Attempt to coerce non-integral real number to an Integer"
674
# floor yields an integer
675
x = x.floor()
676
break
677
elif paritype == t_PADIC:
678
if x._valp() < 0:
679
raise TypeError("Cannot convert p-adic with negative valuation to an integer")
680
# Lifting a PADIC yields an integer
681
x = x.lift()
682
break
683
elif paritype == t_INTMOD:
684
# Lifting an INTMOD yields an integer
685
x = x.lift()
686
break
687
elif paritype == t_POLMOD:
688
x = x.lift()
689
else:
690
raise TypeError, "Unable to coerce PARI %s to an Integer"%x
691
692
# Now we have a true PARI integer, convert it to Sage
693
t_INT_to_ZZ(self.value, (<pari_gen>x).g)
694
695
elif PyString_Check(x) or PY_TYPE_CHECK(x,unicode):
696
if base < 0 or base > 36:
697
raise ValueError, "`base` (=%s) must be between 2 and 36."%base
698
ibase = base
699
xs = x
700
if xs[0] == c'+':
701
xs += 1
702
if mpz_set_str(self.value, xs, ibase) != 0:
703
raise TypeError, "unable to convert x (=%s) to an integer"%x
704
705
elif PyObject_HasAttrString(x, "_integer_"):
706
# TODO: Note that PyObject_GetAttrString returns NULL if
707
# the attribute was not found. If we could test for this,
708
# we could skip the double lookup. Unfortunately Cython doesn't
709
# seem to let us do this; it flags an error if the function
710
# returns NULL, because it can't construct an "object" object
711
# out of the NULL pointer. This really sucks. Perhaps we could
712
# make the function prototype have return type void*, but
713
# then how do we make Cython handle the reference counting?
714
set_from_Integer(self, (<object> PyObject_GetAttrString(x, "_integer_"))(the_integer_ring))
715
716
elif (PY_TYPE_CHECK(x, list) or PY_TYPE_CHECK(x, tuple)) and base > 1:
717
b = the_integer_ring(base)
718
if b == 2: # we use a faster method
719
for j from 0 <= j < len(x):
720
otmp = x[j]
721
if not PY_TYPE_CHECK(otmp, Integer):
722
# should probably also have fast code for Python ints...
723
otmp = Integer(otmp)
724
if mpz_cmp_si((<Integer>otmp).value, 1) == 0:
725
mpz_setbit(self.value, j)
726
elif mpz_sgn((<Integer>otmp).value) != 0:
727
# one of the entries was something other than 0 or 1.
728
break
729
else:
730
return
731
tmp = the_integer_ring(0)
732
for i in range(len(x)):
733
tmp += the_integer_ring(x[i])*b**i
734
mpz_set(self.value, tmp.value)
735
736
else:
737
import numpy
738
if isinstance(x, numpy.integer):
739
mpz_set_pylong(self.value, x.__long__())
740
return
741
742
elif PY_TYPE_CHECK(x, Element):
743
try:
744
lift = x.lift()
745
if lift._parent != (<Element>x)._parent:
746
tmp = the_integer_ring(lift)
747
mpz_swap(tmp.value, self.value)
748
return
749
except AttributeError:
750
pass
751
752
raise TypeError, "unable to coerce %s to an integer" % type(x)
753
754
def __reduce__(self):
755
"""
756
This is used when pickling integers.
757
758
EXAMPLES::
759
760
sage: n = 5
761
sage: t = n.__reduce__(); t
762
(<built-in function make_integer>, ('5',))
763
sage: t[0](*t[1])
764
5
765
sage: loads(dumps(n)) == n
766
True
767
"""
768
# This single line below took me HOURS to figure out.
769
# It is the *trick* needed to pickle Cython extension types.
770
# The trick is that you must put a pure Python function
771
# as the first argument, and that function must return
772
# the result of unpickling with the argument in the second
773
# tuple as input. All kinds of problems happen
774
# if we don't do this.
775
return sage.rings.integer.make_integer, (self.str(32),)
776
777
cdef _reduce_set(self, s):
778
"""
779
Set this integer from a string in base 32.
780
781
.. note::
782
783
Integers are supposed to be immutable, so you should not
784
use this function.
785
"""
786
mpz_set_str(self.value, s, 32)
787
788
def __index__(self):
789
"""
790
Needed so integers can be used as list indices.
791
792
EXAMPLES::
793
794
sage: v = [1,2,3,4,5]
795
sage: v[Integer(3)]
796
4
797
sage: v[Integer(2):Integer(4)]
798
[3, 4]
799
"""
800
return mpz_get_pyintlong(self.value)
801
802
def _im_gens_(self, codomain, im_gens):
803
"""
804
Return the image of self under the map that sends the generators of
805
the parent to im_gens. Since ZZ maps canonically in the category
806
of rings, this is just the natural coercion.
807
808
EXAMPLES::
809
810
sage: n = -10
811
sage: R = GF(17)
812
sage: n._im_gens_(R, [R(1)])
813
7
814
"""
815
return codomain._coerce_(self)
816
817
cdef _xor(Integer self, Integer other):
818
cdef Integer x
819
x = PY_NEW(Integer)
820
mpz_xor(x.value, self.value, other.value)
821
return x
822
823
def __xor__(x, y):
824
"""
825
Compute the exclusive or of x and y.
826
827
EXAMPLES::
828
829
sage: n = ZZ(2); m = ZZ(3)
830
sage: n.__xor__(m)
831
1
832
"""
833
if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
834
return (<Integer>x)._xor(y)
835
return bin_op(x, y, operator.xor)
836
837
def __richcmp__(left, right, int op):
838
"""
839
cmp for integers
840
841
EXAMPLES::
842
843
sage: 2 < 3
844
True
845
sage: 2 > 3
846
False
847
sage: 2 == 3
848
False
849
sage: 3 > 2
850
True
851
sage: 3 < 2
852
False
853
854
sage: 1000000000000000000000000000000000000000000000000000.0r==1000000000000000000000000000000000000000000000000000
855
False
856
sage: 1000000000000000000000000000000000000000000000000000.1r==1000000000000000000000000000000000000000000000000000
857
False
858
859
Canonical coercions are used but non-canonical ones are not.
860
861
::
862
863
sage: 4 == 4/1
864
True
865
sage: 4 == '4'
866
False
867
868
TESTS::
869
870
sage: 3 < 4r
871
True
872
sage: 3r < 4
873
True
874
sage: 3 >= 4r
875
False
876
sage: 4r <= 3
877
False
878
sage: 12345678901234567890123456789r == 12345678901234567890123456789
879
True
880
sage: 12345678901234567890123456788 < 12345678901234567890123456789r
881
True
882
sage: 2 < 2.7r
883
True
884
sage: 4 < 3.1r
885
False
886
sage: -1r < 1
887
True
888
sage: -1.5r < 3
889
True
890
sage: Ilist = [-2,-1,0,1,2,12345678901234567890123456788]
891
sage: ilist = [-4r,-1r,0r,1r,2r,5r]
892
sage: llist = [-12345678901234567890123456788r, 12345678901234567890123456788r, 12345678901234567890123456900r]
893
sage: flist = [-21.8r, -1.2r, -.000005r, 0.0r, .999999r, 1000000000000.0r]
894
sage: all([(a < b) == (RR(a) < RR(b)) for (a, b) in zip(Ilist, ilist)])
895
True
896
sage: all([(a > b) == (RR(a) > RR(b)) for (a, b) in zip(Ilist, ilist)])
897
True
898
sage: all([(a == b) == (RR(a) == RR(b)) for (a, b) in zip(Ilist, ilist)])
899
True
900
sage: all([(a <= b) == (RR(a) <= RR(b)) for (a, b) in zip(Ilist, ilist)])
901
True
902
sage: all([(a >= b) == (RR(a) >= RR(b)) for (a, b) in zip(Ilist, ilist)])
903
True
904
sage: all([(a != b) == (RR(a) != RR(b)) for (a, b) in zip(Ilist, ilist)])
905
True
906
sage: all([(a < b) == (RR(a) < RR(b)) for (a, b) in zip(Ilist, llist)])
907
True
908
sage: all([(a > b) == (RR(a) > RR(b)) for (a, b) in zip(Ilist, llist)])
909
True
910
sage: all([(a < b) == (RR(a) < RR(b)) for (a, b) in zip(Ilist, flist)])
911
True
912
sage: all([(a > b) == (RR(a) > RR(b)) for (a, b) in zip(Ilist, flist)])
913
True
914
915
Verify that trac 12149 was fixed (and the fix is consistent
916
with Python ints)::
917
918
sage: a = int(1); b = 1; n = float('nan')
919
sage: a == n
920
False
921
sage: a == n, b == n
922
(False, False)
923
sage: a != n, b != n, n != b
924
(True, True, True)
925
sage: a < n, b < n, n > b
926
(False, False, False)
927
sage: a > n, b > n, n < b
928
(False, False, False)
929
sage: a <= n, b <= n, n >= b
930
(False, False, False)
931
sage: a >= n, b >= n, n <= b
932
(False, False, False)
933
"""
934
cdef int c
935
cdef double d
936
if PY_TYPE_CHECK(left, Integer):
937
if PY_TYPE_CHECK(right, Integer):
938
c = mpz_cmp((<Integer>left).value, (<Integer>right).value)
939
elif PyInt_CheckExact(right):
940
c = mpz_cmp_si((<Integer>left).value, PyInt_AS_LONG(right))
941
elif PyLong_CheckExact(right):
942
mpz_set_pylong(mpz_tmp, right)
943
c = mpz_cmp((<Integer>left).value, mpz_tmp)
944
elif PyFloat_CheckExact(right):
945
d = PyFloat_AsDouble(right)
946
if isnan(d): return op == 3
947
c = mpz_cmp_d((<Integer>left).value, d)
948
else:
949
return (<sage.structure.element.Element>left)._richcmp(right, op)
950
951
else: # right is an Integer
952
if PyInt_CheckExact(left):
953
c = -mpz_cmp_si((<Integer>right).value, PyInt_AS_LONG(left))
954
elif PyLong_CheckExact(left):
955
mpz_set_pylong(mpz_tmp, left)
956
c = mpz_cmp(mpz_tmp, (<Integer>right).value)
957
elif PyFloat_CheckExact(left):
958
d = PyFloat_AsDouble(left)
959
if isnan(d): return op == 3
960
c = -mpz_cmp_d((<Integer>right).value, d)
961
else:
962
return (<sage.structure.element.Element>left)._richcmp(right, op)
963
return (<sage.structure.element.Element>left)._rich_to_bool(op, c)
964
965
cdef int _cmp_c_impl(left, sage.structure.element.Element right) except -2:
966
cdef int i
967
i = mpz_cmp((<Integer>left).value, (<Integer>right).value)
968
if i < 0: return -1
969
elif i == 0: return 0
970
else: return 1
971
972
def __copy__(self):
973
"""
974
Return a copy of the integer.
975
976
EXAMPLES::
977
978
sage: n = 2
979
sage: copy(n)
980
2
981
sage: copy(n) is n
982
False
983
"""
984
cdef Integer z
985
z = PY_NEW(Integer)
986
set_mpz(z,self.value)
987
return z
988
989
def list(self):
990
"""
991
Return a list with this integer in it, to be compatible with the
992
method for number fields.
993
994
EXAMPLES::
995
996
sage: m = 5
997
sage: m.list()
998
[5]
999
"""
1000
return [ self ]
1001
1002
def __dealloc__(self):
1003
mpz_clear(self.value)
1004
1005
def __repr__(self):
1006
"""
1007
Return string representation of this integer.
1008
1009
EXAMPLES::
1010
1011
sage: n = -5; n.__repr__()
1012
'-5'
1013
"""
1014
return self.str()
1015
1016
def _latex_(self):
1017
"""
1018
Return latex representation of this integer. This is just the
1019
underlying string representation and nothing more. This is called
1020
by the latex function.
1021
1022
EXAMPLES::
1023
1024
sage: n = -5; n._latex_()
1025
'-5'
1026
sage: latex(n)
1027
-5
1028
"""
1029
return self.str()
1030
1031
def _sympy_(self):
1032
"""
1033
Convert Sage Integer() to SymPy Integer.
1034
1035
EXAMPLES::
1036
1037
sage: n = 5; n._sympy_()
1038
5
1039
sage: n = -5; n._sympy_()
1040
-5
1041
"""
1042
import sympy
1043
return sympy.sympify(int(self))
1044
1045
def _mathml_(self):
1046
"""
1047
Return mathml representation of this integer.
1048
1049
EXAMPLES::
1050
1051
sage: mathml(-45)
1052
<mn>-45</mn>
1053
sage: (-45)._mathml_()
1054
'<mn>-45</mn>'
1055
"""
1056
return '<mn>%s</mn>'%self
1057
1058
def str(self, int base=10):
1059
r"""
1060
Return the string representation of ``self`` in the
1061
given base.
1062
1063
EXAMPLES::
1064
1065
sage: Integer(2^10).str(2)
1066
'10000000000'
1067
sage: Integer(2^10).str(17)
1068
'394'
1069
1070
::
1071
1072
sage: two=Integer(2)
1073
sage: two.str(1)
1074
Traceback (most recent call last):
1075
...
1076
ValueError: base (=1) must be between 2 and 36
1077
1078
::
1079
1080
sage: two.str(37)
1081
Traceback (most recent call last):
1082
...
1083
ValueError: base (=37) must be between 2 and 36
1084
1085
::
1086
1087
sage: big = 10^5000000
1088
sage: s = big.str() # long time (> 20 seconds)
1089
sage: len(s) # long time (depends on above defn of s)
1090
5000001
1091
sage: s[:10] # long time (depends on above defn of s)
1092
'1000000000'
1093
"""
1094
if base < 2 or base > 36:
1095
raise ValueError, "base (=%s) must be between 2 and 36"%base
1096
cdef size_t n
1097
cdef char *s
1098
n = mpz_sizeinbase(self.value, base) + 2
1099
s = <char *>PyMem_Malloc(n)
1100
if s == NULL:
1101
raise MemoryError, "Unable to allocate enough memory for the string representation of an integer."
1102
sig_on()
1103
mpz_get_str(s, base, self.value)
1104
sig_off()
1105
k = <object> PyString_FromString(s)
1106
PyMem_Free(s)
1107
return k
1108
1109
def __format__(self, *args, **kwargs):
1110
"""
1111
Returns a string representation using Python's Format protocol.
1112
Valid format descriptions are exactly those for Python integers.
1113
1114
EXAMPLES::
1115
1116
sage: "{0:#x}; {0:#b}; {0:+05d}".format(ZZ(17))
1117
'0x11; 0b10001; +0017'
1118
1119
"""
1120
return int(self).__format__(*args,**kwargs)
1121
1122
def ordinal_str(self):
1123
"""
1124
Returns a string representation of the ordinal associated to self.
1125
1126
EXAMPLES::
1127
1128
sage: [ZZ(n).ordinal_str() for n in range(25)]
1129
['0th',
1130
'1st',
1131
'2nd',
1132
'3rd',
1133
'4th',
1134
...
1135
'10th',
1136
'11th',
1137
'12th',
1138
'13th',
1139
'14th',
1140
...
1141
'20th',
1142
'21st',
1143
'22nd',
1144
'23rd',
1145
'24th']
1146
1147
sage: ZZ(1001).ordinal_str()
1148
'1001st'
1149
1150
sage: ZZ(113).ordinal_str()
1151
'113th'
1152
sage: ZZ(112).ordinal_str()
1153
'112th'
1154
sage: ZZ(111).ordinal_str()
1155
'111th'
1156
1157
"""
1158
if self<0:
1159
raise ValueError, "Negative integers are not ordinals."
1160
n = self.abs()
1161
if ((n%100)!=11 and n%10==1):
1162
th = 'st'
1163
elif ((n%100)!=12 and n%10==2):
1164
th = 'nd'
1165
elif ((n%100)!=13 and n%10==3):
1166
th = 'rd'
1167
else:
1168
th = 'th'
1169
return n.str()+th
1170
1171
def __hex__(self):
1172
r"""
1173
Return the hexadecimal digits of self in lower case.
1174
1175
.. note::
1176
1177
'0x' is *not* prepended to the result like is done by the
1178
corresponding Python function on int or long. This is for
1179
efficiency sake--adding and stripping the string wastes
1180
time; since this function is used for conversions from
1181
integers to other C-library structures, it is important
1182
that it be fast.
1183
1184
EXAMPLES::
1185
1186
sage: print hex(Integer(15))
1187
f
1188
sage: print hex(Integer(16))
1189
10
1190
sage: print hex(Integer(16938402384092843092843098243))
1191
36bb1e3929d1a8fe2802f083
1192
sage: print hex(long(16938402384092843092843098243))
1193
0x36bb1e3929d1a8fe2802f083L
1194
"""
1195
return self.str(16)
1196
1197
def __oct__(self):
1198
r"""
1199
Return the digits of self in base 8.
1200
1201
.. note::
1202
1203
'0' is *not* prepended to the result like is done by the
1204
corresponding Python function on int or long. This is for
1205
efficiency sake--adding and stripping the string wastes
1206
time; since this function is used for conversions from
1207
integers to other C-library structures, it is important
1208
that it be fast.
1209
1210
EXAMPLES::
1211
1212
sage: print oct(Integer(800))
1213
1440
1214
sage: print oct(Integer(8))
1215
10
1216
sage: print oct(Integer(-50))
1217
-62
1218
sage: print oct(Integer(-899))
1219
-1603
1220
sage: print oct(Integer(16938402384092843092843098243))
1221
15535436162247215217705000570203
1222
1223
Behavior of Sage integers vs. Python integers::
1224
1225
sage: oct(Integer(10))
1226
'12'
1227
sage: oct(int(10))
1228
'012'
1229
sage: oct(Integer(-23))
1230
'-27'
1231
sage: oct(int(-23))
1232
'-027'
1233
"""
1234
return self.str(8)
1235
1236
def binary(self):
1237
"""
1238
Return the binary digits of self as a string.
1239
1240
EXAMPLES::
1241
1242
sage: print Integer(15).binary()
1243
1111
1244
sage: print Integer(16).binary()
1245
10000
1246
sage: print Integer(16938402384092843092843098243).binary()
1247
1101101011101100011110001110010010100111010001101010001111111000101000000000101111000010000011
1248
"""
1249
return self.str(2)
1250
1251
def bits(self):
1252
"""
1253
Return the bits in self as a list, least significant first.
1254
1255
EXAMPLES::
1256
1257
sage: 500.bits()
1258
[0, 0, 1, 0, 1, 1, 1, 1, 1]
1259
sage: 11.bits()
1260
[1, 1, 0, 1]
1261
"""
1262
return self.digits(base=2)
1263
1264
def nbits(self):
1265
"""
1266
Return the number of bits in self.
1267
1268
EXAMPLES::
1269
1270
sage: 500.nbits()
1271
9
1272
sage: 5.nbits()
1273
3
1274
sage: 0.nbits() == len(0.bits()) == 0.ndigits(base=2)
1275
True
1276
sage: 12345.nbits() == len(12345.binary())
1277
True
1278
"""
1279
# mpz_sizeinbase(0,2) always returns 1
1280
if mpz_cmp_si(self.value,0) == 0:
1281
return int(0)
1282
else:
1283
return int(mpz_sizeinbase(self.value, 2))
1284
1285
def trailing_zero_bits(self):
1286
"""
1287
Return the number of trailing zero bits in self, i.e.
1288
the exponent of the largest power of 2 dividing self.
1289
1290
EXAMPLES::
1291
1292
sage: 11.trailing_zero_bits()
1293
0
1294
sage: (-11).trailing_zero_bits()
1295
0
1296
sage: (11<<5).trailing_zero_bits()
1297
5
1298
sage: (-11<<5).trailing_zero_bits()
1299
5
1300
sage: 0.trailing_zero_bits()
1301
0
1302
1303
"""
1304
if mpz_sgn(self.value) == 0:
1305
return int(0)
1306
return int(mpz_scan1(self.value, 0))
1307
1308
def digits(self, base=10, digits=None, padto=0):
1309
r"""
1310
Return a list of digits for ``self`` in the given base in little
1311
endian order.
1312
1313
The returned value is unspecified if self is a negative number
1314
and the digits are given.
1315
1316
INPUT:
1317
1318
- ``base`` - integer (default: 10)
1319
1320
- ``digits`` - optional indexable object as source for
1321
the digits
1322
1323
- ``padto`` - the minimal length of the returned list,
1324
sufficient number of zeros are added to make the list minimum that
1325
length (default: 0)
1326
1327
As a shorthand for ``digits(2)``, you can use :meth:`.bits`.
1328
1329
Also see :meth:`ndigits`.
1330
1331
EXAMPLES::
1332
1333
sage: 17.digits()
1334
[7, 1]
1335
sage: 5.digits(base=2, digits=["zero","one"])
1336
['one', 'zero', 'one']
1337
sage: 5.digits(3)
1338
[2, 1]
1339
sage: 0.digits(base=10) # 0 has 0 digits
1340
[]
1341
sage: 0.digits(base=2) # 0 has 0 digits
1342
[]
1343
sage: 10.digits(16,'0123456789abcdef')
1344
['a']
1345
sage: 0.digits(16,'0123456789abcdef')
1346
[]
1347
sage: 0.digits(16,'0123456789abcdef',padto=1)
1348
['0']
1349
sage: 123.digits(base=10,padto=5)
1350
[3, 2, 1, 0, 0]
1351
sage: 123.digits(base=2,padto=3) # padto is the minimal length
1352
[1, 1, 0, 1, 1, 1, 1]
1353
sage: 123.digits(base=2,padto=10,digits=(1,-1))
1354
[-1, -1, 1, -1, -1, -1, -1, 1, 1, 1]
1355
sage: a=9939082340; a.digits(10)
1356
[0, 4, 3, 2, 8, 0, 9, 3, 9, 9]
1357
sage: a.digits(512)
1358
[100, 302, 26, 74]
1359
sage: (-12).digits(10)
1360
[-2, -1]
1361
sage: (-12).digits(2)
1362
[0, 0, -1, -1]
1363
1364
We support large bases.
1365
1366
::
1367
1368
sage: n=2^6000
1369
sage: n.digits(2^3000)
1370
[0, 0, 1]
1371
1372
::
1373
1374
sage: base=3; n=25
1375
sage: l=n.digits(base)
1376
sage: # the next relationship should hold for all n,base
1377
sage: sum(base^i*l[i] for i in range(len(l)))==n
1378
True
1379
sage: base=3; n=-30; l=n.digits(base); sum(base^i*l[i] for i in range(len(l)))==n
1380
True
1381
1382
The inverse of this method -- constructing an integer from a
1383
list of digits and a base -- can be done using the above method
1384
or by simply using :class:`ZZ()
1385
<sage.rings.integer_ring.IntegerRing_class>` with a base::
1386
1387
sage: x = 123; ZZ(x.digits(), 10)
1388
123
1389
sage: x == ZZ(x.digits(6), 6)
1390
True
1391
sage: x == ZZ(x.digits(25), 25)
1392
True
1393
1394
Using :func:`sum` and :func:`enumerate` to do the same thing is
1395
slightly faster in many cases (and
1396
:func:`~sage.misc.misc_c.balanced_sum` may be faster yet). Of
1397
course it gives the same result::
1398
1399
sage: base = 4
1400
sage: sum(digit * base^i for i, digit in enumerate(x.digits(base))) == ZZ(x.digits(base), base)
1401
True
1402
1403
Note: In some cases it is faster to give a digits collection. This
1404
would be particularly true for computing the digits of a series of
1405
small numbers. In these cases, the code is careful to allocate as
1406
few python objects as reasonably possible.
1407
1408
::
1409
1410
sage: digits = range(15)
1411
sage: l=[ZZ(i).digits(15,digits) for i in range(100)]
1412
sage: l[16]
1413
[1, 1]
1414
1415
This function is comparable to ``str`` for speed.
1416
1417
::
1418
1419
sage: n=3^100000
1420
sage: n.digits(base=10)[-1] # slightly slower than str
1421
1
1422
sage: n=10^10000
1423
sage: n.digits(base=10)[-1] # slightly faster than str
1424
1
1425
1426
AUTHORS:
1427
1428
- Joel B. Mohler (2008-03-02): significantly rewrote this entire function
1429
"""
1430
cdef Integer _base
1431
cdef Integer self_abs = self
1432
cdef int power_index = 0
1433
cdef list power_list
1434
cdef list l
1435
cdef int i
1436
cdef size_t s
1437
1438
if PY_TYPE_CHECK(base, Integer):
1439
_base = <Integer>base
1440
else:
1441
_base = Integer(base)
1442
1443
if mpz_cmp_si(_base.value,2) < 0:
1444
raise ValueError, "base must be >= 2"
1445
1446
if mpz_sgn(self.value) < 0:
1447
self_abs = -self
1448
1449
cdef bint do_sig_on
1450
if mpz_sgn(self.value) == 0:
1451
l = [the_integer_ring._zero_element if digits is None else digits[0]]*padto
1452
elif mpz_cmp_si(_base.value,2) == 0:
1453
s = mpz_sizeinbase(self.value, 2)
1454
if digits:
1455
o = digits[1]
1456
z = digits[0]
1457
else:
1458
if mpz_sgn(self.value) == 1:
1459
o = the_integer_ring._one_element
1460
else:
1461
o = -the_integer_ring._one_element
1462
z = the_integer_ring._zero_element
1463
l = [z]*(s if s >= padto else padto)
1464
for i from 0<= i < s:
1465
# mpz_tstbit seems to return 0 for the high-order bit of
1466
# negative numbers?!
1467
if mpz_tstbit(self_abs.value,i):
1468
l[i] = o
1469
else:
1470
s = mpz_sizeinbase(self.value, 2)
1471
do_sig_on = (s > 256)
1472
if do_sig_on: sig_on()
1473
1474
# We use a divide and conquer approach (suggested by the prior
1475
# author, malb?, of the digits method) here: for base b, compute
1476
# b^2, b^4, b^8, ... (repeated squaring) until you get larger
1477
# than your number; then compute (n // b^256, n % b^256)
1478
# (if b^512 > number) to split the number in half and recurse
1479
1480
# Pre-computing the exact number of digits up-front is actually
1481
# faster (especially for large values of self) than trimming off
1482
# trailing zeros after the fact. It also seems that it would
1483
# avoid duplicating the list in memory with a list-slice.
1484
z = the_integer_ring._zero_element if digits is None else digits[0]
1485
s = self_abs.exact_log(_base)
1486
l = [z]*(s+1 if s+1 >= padto else padto)
1487
1488
# set up digits for optimal access once we get inside the worker
1489
# functions
1490
if not digits is None:
1491
# list objects have fastest access in the innermost loop
1492
if not PyList_CheckExact(digits):
1493
digits = [digits[i] for i in range(_base)]
1494
elif mpz_cmp_ui(_base.value,s) < 0 and mpz_cmp_ui(_base.value,10000):
1495
# We can get a speed boost by pre-allocating digit values in
1496
# big cases.
1497
# We do this we have more digits than the base and the base
1498
# is not too extremely large (currently, "extremely" means
1499
# larger than 10000 -- that's very arbitrary.)
1500
if mpz_sgn(self.value) > 0:
1501
digits = [Integer(i) for i in range(_base)]
1502
else:
1503
# All the digits will be negated in the recursive function.
1504
# we'll just compensate for python index semantics
1505
digits = [Integer(i) for i in range(-_base,0)]
1506
digits[0] = the_integer_ring._zero_element
1507
1508
if s < 40:
1509
_digits_naive(self.value,l,0,_base,digits)
1510
else:
1511
# count the bits of s
1512
i = 0
1513
while s != 0:
1514
s >>= 1
1515
i += 1
1516
1517
power_list = [_base]*i
1518
for power_index from 1 <= power_index < i:
1519
power_list[power_index] = power_list[power_index-1]**2
1520
1521
# Note that it may appear that the recursive calls to
1522
# _digit_internal would be assigning list elements i in l for
1523
# anywhere from 0<=i<(1<<power_index). However, this is not
1524
# the case due to the optimization of skipping assigns
1525
# assigning zero.
1526
_digits_internal(self.value,l,0,i-1,power_list,digits)
1527
1528
if do_sig_on: sig_off()
1529
1530
# padding should be taken care of with-in the function
1531
# all we need to do is return
1532
return l
1533
1534
def ndigits(self, base=10):
1535
"""
1536
Return the number of digits of self expressed in the given base.
1537
1538
INPUT:
1539
1540
- ``base`` - integer (default: 10)
1541
1542
EXAMPLES::
1543
1544
sage: n = 52
1545
sage: n.ndigits()
1546
2
1547
sage: n = -10003
1548
sage: n.ndigits()
1549
5
1550
sage: n = 15
1551
sage: n.ndigits(2)
1552
4
1553
sage: n = 1000**1000000+1
1554
sage: n.ndigits()
1555
3000001
1556
sage: n = 1000**1000000-1
1557
sage: n.ndigits()
1558
3000000
1559
sage: n = 10**10000000-10**9999990
1560
sage: n.ndigits()
1561
10000000
1562
"""
1563
cdef Integer temp
1564
1565
if mpz_sgn(self.value) == 0:
1566
temp = PY_NEW(Integer)
1567
mpz_set_ui(temp.value, 0)
1568
return temp
1569
1570
if mpz_sgn(self.value) > 0:
1571
temp = self.exact_log(base)
1572
mpz_add_ui(temp.value, temp.value, 1)
1573
return temp
1574
else:
1575
return self.abs().exact_log(base) + 1
1576
1577
cdef void set_from_mpz(Integer self, mpz_t value):
1578
mpz_set(self.value, value)
1579
1580
cdef mpz_t* get_value(Integer self):
1581
return &self.value
1582
1583
cdef void _to_ZZ(self, ZZ_c *z):
1584
sig_on()
1585
mpz_to_ZZ(z, &self.value)
1586
sig_off()
1587
1588
cpdef ModuleElement _add_(self, ModuleElement right):
1589
"""
1590
Integer addition.
1591
1592
TESTS::
1593
1594
sage: Integer(32) + Integer(23)
1595
55
1596
sage: sum(Integer(i) for i in [1..100])
1597
5050
1598
sage: a = ZZ.random_element(10^50000)
1599
sage: b = ZZ.random_element(10^50000)
1600
sage: a+b == b+a
1601
True
1602
"""
1603
# self and right are guaranteed to be Integers
1604
cdef Integer x = <Integer>PY_NEW(Integer)
1605
mpz_add(x.value, self.value, (<Integer>right).value)
1606
return x
1607
1608
cpdef ModuleElement _iadd_(self, ModuleElement right):
1609
# self and right are guaranteed to be Integers, self safe to mutate
1610
mpz_add(self.value, self.value, (<Integer>right).value)
1611
return self
1612
1613
cdef RingElement _add_long(self, long n):
1614
"""
1615
Fast path for adding a C long.
1616
1617
TESTS::
1618
sage: int(10) + Integer(100)
1619
110
1620
sage: Integer(100) + int(10)
1621
110
1622
sage: Integer(10^100) + int(10)
1623
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010
1624
1625
Also called for subtraction::
1626
1627
sage: Integer(100) - int(10)
1628
90
1629
sage: Integer(10^100) - int(10)
1630
9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999990
1631
1632
Make sure it works when -<long>n would overflow::
1633
1634
sage: most_neg_long = int(-sys.maxint - 1)
1635
sage: type(most_neg_long), type(-most_neg_long)
1636
(<type 'int'>, <type 'long'>)
1637
sage: 0 + most_neg_long == most_neg_long
1638
True
1639
sage: 0 - most_neg_long == -most_neg_long
1640
True
1641
"""
1642
cdef Integer x = <Integer>PY_NEW(Integer)
1643
if n > 0:
1644
mpz_add_ui(x.value, self.value, n)
1645
else:
1646
# Note that 0-<unsigned long>n is always -n as an unsigned
1647
# long (whereas -n may overflow).
1648
mpz_sub_ui(x.value, self.value, 0 - <unsigned long>n)
1649
return x
1650
1651
cpdef ModuleElement _sub_(self, ModuleElement right):
1652
"""
1653
Integer subtraction.
1654
1655
TESTS::
1656
1657
sage: Integer(32) - Integer(23)
1658
9
1659
sage: Integer(10^100) - Integer(1)
1660
9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
1661
sage: Integer(1) - Integer(10^100)
1662
-9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
1663
sage: a = ZZ.random_element(10^50000)
1664
sage: b = ZZ.random_element(10^50000)
1665
sage: a-b == -(b-a) == a + -b
1666
True
1667
"""
1668
# self and right are guaranteed to be Integers
1669
cdef Integer x = <Integer>PY_NEW(Integer)
1670
mpz_sub(x.value, self.value, (<Integer>right).value)
1671
return x
1672
1673
cpdef ModuleElement _isub_(self, ModuleElement right):
1674
mpz_sub(self.value, self.value, (<Integer>right).value)
1675
return self
1676
1677
def __neg__(self):
1678
"""
1679
TESTS::
1680
1681
sage: a = Integer(3)
1682
sage: -a
1683
-3
1684
sage: a = Integer(3^100); a
1685
515377520732011331036461129765621272702107522001
1686
sage: -a
1687
-515377520732011331036461129765621272702107522001
1688
"""
1689
cdef Integer x = <Integer>PY_NEW(Integer)
1690
mpz_neg(x.value, self.value)
1691
return x
1692
1693
cpdef ModuleElement _neg_(self):
1694
cdef Integer x = <Integer>PY_NEW(Integer)
1695
mpz_neg(x.value, self.value)
1696
return x
1697
1698
cpdef _act_on_(self, s, bint self_on_left):
1699
"""
1700
EXAMPLES::
1701
1702
sage: 8 * [0] #indirect doctest
1703
[0, 0, 0, 0, 0, 0, 0, 0]
1704
sage: 'hi' * 8
1705
'hihihihihihihihi'
1706
"""
1707
if isinstance(s, (list, tuple, basestring)):
1708
if mpz_fits_slong_p(self.value):
1709
return s * mpz_get_si(self.value)
1710
else:
1711
return s * int(self) # will raise the appropriate exception
1712
1713
cdef ModuleElement _mul_long(self, long n):
1714
"""
1715
Fast path for multiplying a C long.
1716
1717
TESTS::
1718
1719
sage: Integer(25) * int(4)
1720
100
1721
sage: int(4) * Integer(25)
1722
100
1723
sage: Integer(10^100) * int(4)
1724
40000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
1725
"""
1726
cdef Integer x = <Integer>PY_NEW(Integer)
1727
if mpz_size(self.value) > 100000:
1728
sig_on()
1729
mpz_mul_si(x.value, self.value, n)
1730
sig_off()
1731
else:
1732
mpz_mul_si(x.value, self.value, n)
1733
return x
1734
1735
cpdef RingElement _mul_(self, RingElement right):
1736
"""
1737
Integer multiplication.
1738
1739
sage: Integer(25) * Integer(4)
1740
100
1741
sage: Integer(5^100) * Integer(2^100)
1742
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
1743
sage: a = ZZ.random_element(10^50000)
1744
sage: b = ZZ.random_element(10^50000)
1745
sage: a*b == b*a
1746
True
1747
"""
1748
# self and right are guaranteed to be Integers
1749
cdef Integer x = <Integer>PY_NEW(Integer)
1750
if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
1751
# We only use the signal handler (to enable ctrl-c out) when the
1752
# product might take a while to compute
1753
sig_on()
1754
mpz_mul(x.value, self.value, (<Integer>right).value)
1755
sig_off()
1756
else:
1757
mpz_mul(x.value, self.value, (<Integer>right).value)
1758
return x
1759
1760
cpdef RingElement _imul_(self, RingElement right):
1761
if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
1762
# We only use the signal handler (to enable ctrl-c out) when the
1763
# product might take a while to compute
1764
sig_on()
1765
mpz_mul(self.value, self.value, (<Integer>right).value)
1766
sig_off()
1767
else:
1768
mpz_mul(self.value, self.value, (<Integer>right).value)
1769
return self
1770
1771
cpdef RingElement _div_(self, RingElement right):
1772
r"""
1773
Computes `\frac{a}{b}`
1774
1775
EXAMPLES::
1776
1777
sage: a = Integer(3) ; b = Integer(4)
1778
sage: a / b == Rational(3) / 4
1779
True
1780
sage: Integer(32) / Integer(32)
1781
1
1782
"""
1783
# This is vastly faster than doing it here, since here
1784
# we can't cimport rationals.
1785
return the_integer_ring._div(self, right)
1786
1787
def __floordiv__(x, y):
1788
r"""
1789
Computes the whole part of `\frac{x}{y}`.
1790
1791
EXAMPLES::
1792
1793
sage: a = Integer(321) ; b = Integer(10)
1794
sage: a // b
1795
32
1796
sage: z = Integer(-231)
1797
sage: z // 2
1798
-116
1799
sage: z = Integer(231)
1800
sage: z // 2
1801
115
1802
sage: z // -2
1803
-116
1804
sage: z // 0
1805
Traceback (most recent call last):
1806
...
1807
ZeroDivisionError: Integer division by zero
1808
sage: 101 // int(5)
1809
20
1810
sage: 100 // int(-3)
1811
-34
1812
1813
TESTS::
1814
1815
sage: signs = [(11,5), (11,-5), (-11,5), (-11,-5)]
1816
sage: control = [int(a) // int(b) for a, b in signs]
1817
sage: [a // b for a,b in signs] == control
1818
True
1819
sage: [a // int(b) for a,b in signs] == control
1820
True
1821
sage: [int(a) // b for a,b in signs] == control
1822
True
1823
"""
1824
cdef Integer z = <Integer>PY_NEW(Integer)
1825
cdef long yy, res
1826
if PY_TYPE(x) is PY_TYPE(y):
1827
if not mpz_sgn((<Integer>y).value):
1828
raise ZeroDivisionError, "Integer division by zero"
1829
if mpz_size((<Integer>x).value) > 100000:
1830
sig_on()
1831
mpz_fdiv_q(z.value, (<Integer>x).value, (<Integer>y).value)
1832
sig_off()
1833
else:
1834
mpz_fdiv_q(z.value, (<Integer>x).value, (<Integer>y).value)
1835
return z
1836
1837
elif PyInt_CheckExact(y):
1838
yy = PyInt_AS_LONG(y)
1839
if yy > 0:
1840
mpz_fdiv_q_ui(z.value, (<Integer>x).value, yy)
1841
elif yy == 0:
1842
raise ZeroDivisionError, "Integer division by zero"
1843
else:
1844
res = mpz_fdiv_q_ui(z.value, (<Integer>x).value, -yy)
1845
mpz_neg(z.value, z.value)
1846
if res:
1847
mpz_sub_ui(z.value, z.value, 1)
1848
return z
1849
1850
else:
1851
return bin_op(x, y, operator.floordiv)
1852
1853
def __pow__(self, n, modulus):
1854
r"""
1855
Computes `\text{self}^n`
1856
1857
EXAMPLES::
1858
1859
sage: 2^-6
1860
1/64
1861
sage: 2^6
1862
64
1863
sage: 2^0
1864
1
1865
sage: 2^-0
1866
1
1867
sage: (-1)^(1/3)
1868
(-1)^(1/3)
1869
1870
For consistency with Python and MPFR, 0^0 is defined to be 1 in
1871
Sage::
1872
1873
sage: 0^0
1874
1
1875
1876
The base need not be an integer (it can be a builtin Python type).
1877
1878
::
1879
1880
sage: int(2)^10
1881
1024
1882
sage: float(2.5)^10
1883
9536.7431640625
1884
sage: 'sage'^3
1885
'sagesagesage'
1886
1887
The exponent must fit in a long unless the base is -1, 0, or 1.
1888
1889
::
1890
1891
sage: x = 2^100000000000000000000000
1892
Traceback (most recent call last):
1893
...
1894
RuntimeError: exponent must be at most 2147483647 # 32-bit
1895
RuntimeError: exponent must be at most 9223372036854775807 # 64-bit
1896
sage: (-1)^100000000000000000000000
1897
1
1898
1899
We raise 2 to various interesting exponents::
1900
1901
sage: 2^x # symbolic x
1902
2^x
1903
sage: 2^1.5 # real number
1904
2.82842712474619
1905
sage: 2^float(1.5) # python float
1906
2.8284271247461903
1907
sage: 2^I # complex number
1908
2^I
1909
sage: f = 2^(sin(x)-cos(x)); f
1910
2^(sin(x) - cos(x))
1911
sage: f(x=3)
1912
2^(sin(3) - cos(3))
1913
1914
A symbolic sum::
1915
1916
sage: x,y,z = var('x,y,z')
1917
sage: 2^(x+y+z)
1918
2^(x + y + z)
1919
sage: 2^(1/2)
1920
sqrt(2)
1921
sage: 2^(-1/2)
1922
1/2*sqrt(2)
1923
1924
TESTS::
1925
1926
sage: complex(0,1)^2
1927
(-1+0j)
1928
sage: R.<t> = QQ[]
1929
sage: 2^t
1930
Traceback (most recent call last):
1931
...
1932
TypeError: non-integral exponents not supported
1933
sage: int(3)^-3
1934
1/27
1935
sage: type(int(3)^2)
1936
<type 'sage.rings.integer.Integer'>
1937
sage: type(int(3)^int(2))
1938
<type 'int'>
1939
"""
1940
if modulus is not None:
1941
#raise RuntimeError, "__pow__ dummy argument ignored"
1942
from sage.rings.finite_rings.integer_mod import Mod
1943
return Mod(self, modulus) ** n
1944
1945
if not PY_TYPE_CHECK(self, Integer):
1946
if isinstance(self, str):
1947
return self * n
1948
if not PY_TYPE_CHECK(self, int):
1949
return self ** int(n)
1950
else:
1951
self = Integer(self) #convert from int to Integer
1952
cdef Integer _self = <Integer>self
1953
cdef long nn
1954
1955
try:
1956
nn = PyNumber_Index(n)
1957
except TypeError:
1958
s = parent_c(n)(self)
1959
return s**n
1960
1961
except OverflowError:
1962
if mpz_cmp_si(_self.value, 1) == 0:
1963
return self
1964
elif mpz_cmp_si(_self.value, 0) == 0:
1965
return self
1966
elif mpz_cmp_si(_self.value, -1) == 0:
1967
return self if n % 2 else -self
1968
raise RuntimeError, "exponent must be at most %s" % sys.maxint
1969
1970
if nn == 0:
1971
return one
1972
1973
cdef Integer x = PY_NEW(Integer)
1974
1975
sig_on()
1976
mpz_pow_ui(x.value, (<Integer>self).value, nn if nn > 0 else -nn)
1977
sig_off()
1978
1979
if nn < 0:
1980
return ~x
1981
else:
1982
return x
1983
1984
def nth_root(self, int n, bint truncate_mode=0):
1985
r"""
1986
Returns the (possibly truncated) n'th root of self.
1987
1988
INPUT:
1989
1990
- ``n`` - integer >= 1 (must fit in C int type).
1991
1992
- ``truncate_mode`` - boolean, whether to allow truncation if
1993
self is not an n'th power.
1994
1995
OUTPUT:
1996
1997
If truncate_mode is 0 (default), then returns the exact n'th root
1998
if self is an n'th power, or raises a ValueError if it is not.
1999
2000
If truncate_mode is 1, then if either n is odd or self is
2001
positive, returns a pair (root, exact_flag) where root is the
2002
truncated nth root (rounded towards zero) and exact_flag is a
2003
boolean indicating whether the root extraction was exact;
2004
otherwise raises a ValueError.
2005
2006
AUTHORS:
2007
2008
- David Harvey (2006-09-15)
2009
- Interface changed by John Cremona (2009-04-04)
2010
2011
EXAMPLES::
2012
2013
sage: Integer(125).nth_root(3)
2014
5
2015
sage: Integer(124).nth_root(3)
2016
Traceback (most recent call last):
2017
...
2018
ValueError: 124 is not a 3rd power
2019
sage: Integer(124).nth_root(3, truncate_mode=1)
2020
(4, False)
2021
sage: Integer(125).nth_root(3, truncate_mode=1)
2022
(5, True)
2023
sage: Integer(126).nth_root(3, truncate_mode=1)
2024
(5, False)
2025
2026
::
2027
2028
sage: Integer(-125).nth_root(3)
2029
-5
2030
sage: Integer(-125).nth_root(3,truncate_mode=1)
2031
(-5, True)
2032
sage: Integer(-124).nth_root(3,truncate_mode=1)
2033
(-4, False)
2034
sage: Integer(-126).nth_root(3,truncate_mode=1)
2035
(-5, False)
2036
2037
::
2038
2039
sage: Integer(125).nth_root(2, True)
2040
(11, False)
2041
sage: Integer(125).nth_root(3, True)
2042
(5, True)
2043
2044
::
2045
2046
sage: Integer(125).nth_root(-5)
2047
Traceback (most recent call last):
2048
...
2049
ValueError: n (=-5) must be positive
2050
2051
::
2052
2053
sage: Integer(-25).nth_root(2)
2054
Traceback (most recent call last):
2055
...
2056
ValueError: cannot take even root of negative number
2057
2058
::
2059
2060
sage: a=9
2061
sage: a.nth_root(3)
2062
Traceback (most recent call last):
2063
...
2064
ValueError: 9 is not a 3rd power
2065
2066
sage: a.nth_root(22)
2067
Traceback (most recent call last):
2068
...
2069
ValueError: 9 is not a 22nd power
2070
2071
sage: ZZ(2^20).nth_root(21)
2072
Traceback (most recent call last):
2073
...
2074
ValueError: 1048576 is not a 21st power
2075
2076
sage: ZZ(2^20).nth_root(21, truncate_mode=1)
2077
(1, False)
2078
2079
"""
2080
if n < 1:
2081
raise ValueError, "n (=%s) must be positive" % n
2082
if (mpz_sgn(self.value) < 0) and not (n & 1):
2083
raise ValueError, "cannot take even root of negative number"
2084
cdef Integer x
2085
cdef bint is_exact
2086
x = PY_NEW(Integer)
2087
sig_on()
2088
is_exact = mpz_root(x.value, self.value, n)
2089
sig_off()
2090
2091
if truncate_mode:
2092
return x, is_exact
2093
else:
2094
if is_exact:
2095
return x
2096
else:
2097
raise ValueError, "%s is not a %s power"%(self,integer_ring.ZZ(n).ordinal_str())
2098
2099
cpdef size_t _exact_log_log2_iter(self,Integer m):
2100
"""
2101
This is only for internal use only. You should expect it to crash
2102
and burn for negative or other malformed input. In particular, if
2103
the base `2 \leq m < 4` the log2 approximation of m is 1 and certain
2104
input causes endless loops. Along these lines, it is clear that
2105
this function is most useful for m with a relatively large number
2106
of bits.
2107
2108
For ``small`` values (which I'll leave quite ambiguous), this function
2109
is a fast path for exact log computations. Any integer division with
2110
such input tends to dominate the runtime. Thus we avoid division
2111
entirely in this function.
2112
2113
AUTHOR::
2114
2115
- Joel B. Mohler (2009-04-10)
2116
2117
EXAMPLES::
2118
2119
sage: Integer(125)._exact_log_log2_iter(4)
2120
3
2121
sage: Integer(5^150)._exact_log_log2_iter(5)
2122
150
2123
"""
2124
cdef size_t n_log2
2125
cdef size_t m_log2
2126
cdef size_t l_min
2127
cdef size_t l_max
2128
cdef size_t l
2129
cdef Integer result
2130
cdef mpz_t accum
2131
cdef mpz_t temp_exp
2132
2133
if mpz_cmp_si(m.value,4) < 0:
2134
raise ValueError, "This is undefined or possibly non-convergent with this algorithm."
2135
2136
n_log2=mpz_sizeinbase(self.value,2)-1
2137
m_log2=mpz_sizeinbase(m.value,2)-1
2138
l_min=n_log2/(m_log2+1)
2139
l_max=n_log2/m_log2
2140
if l_min != l_max:
2141
sig_on()
2142
mpz_init(accum)
2143
mpz_init(temp_exp)
2144
mpz_set_ui(accum,1)
2145
l = 0
2146
while l_min != l_max:
2147
# print "self=...",m,l_min,l_max
2148
if l_min + 1 == l_max:
2149
mpz_pow_ui(temp_exp,m.value,l_min+1-l)
2150
# This might over-shoot and make accum > self, but
2151
# we'll know that it's only over by a factor of m^1.
2152
mpz_mul(accum,accum,temp_exp)
2153
if mpz_cmp(self.value,accum) >= 0:
2154
l_min += 1
2155
break
2156
mpz_pow_ui(temp_exp,m.value,l_min-l)
2157
mpz_mul(accum,accum,temp_exp)
2158
l = l_min
2159
2160
# Let x=n_log2-(mpz_sizeinbase(accum,2)-1) and y=m_log2.
2161
# Now, with x>0 and y>0, we have the following observation.
2162
# If floor((x-1)/(y+1))=0, then x-1<y+1 which implies that
2163
# x/y<1+2/y.
2164
# So long as y>=2, this means that floor(x/y)<=1. This shows
2165
# that this iteration is forced to converge for input m >= 4.
2166
# If m=3, we can find input so that floor((x-1)/(y+1))=0 and
2167
# floor(x/y)=2 which results in non-convergence.
2168
2169
# We need the additional '-1' in the l_min computation
2170
# because mpz_sizeinbase(accum,2)-1 is smaller than the
2171
# true log_2(accum)
2172
l_min=l+(n_log2-(mpz_sizeinbase(accum,2)-1)-1)/(m_log2+1)
2173
l_max=l+(n_log2-(mpz_sizeinbase(accum,2)-1))/m_log2
2174
mpz_clear(temp_exp)
2175
mpz_clear(accum)
2176
sig_off()
2177
return l_min
2178
2179
cpdef size_t _exact_log_mpfi_log(self,m):
2180
"""
2181
This is only for internal use only. You should expect it to crash
2182
and burn for negative or other malformed input.
2183
2184
I avoid using this function until the input is large. The overhead
2185
associated with computing the floating point log entirely dominates
2186
the runtime for small values. Note that this is most definitely not
2187
an artifact of format conversion. Tricks with log2 approximations
2188
and using exact integer arithmetic are much better for small input.
2189
2190
AUTHOR::
2191
2192
- Joel B. Mohler (2009-04-10)
2193
2194
EXAMPLES::
2195
2196
sage: Integer(125)._exact_log_mpfi_log(3)
2197
4
2198
sage: Integer(5^150)._exact_log_mpfi_log(5)
2199
150
2200
"""
2201
cdef int i
2202
cdef list pow_2_things
2203
cdef int pow_2
2204
cdef size_t upper,lower,middle
2205
2206
import real_mpfi
2207
R=real_mpfi.RIF
2208
2209
rif_self = R(self)
2210
2211
sig_on()
2212
rif_m = R(m)
2213
rif_log = rif_self.log()/rif_m.log()
2214
# upper is *greater* than the answer
2215
try:
2216
upper = rif_log.upper().ceiling()
2217
except:
2218
# ceiling is probably Infinity
2219
# I'm not sure what to do now
2220
upper = 0
2221
lower = rif_log.lower().floor()
2222
# since the log function is monotonic increasing, lower
2223
# and upper bracket our desired answer
2224
2225
# if upper - lower == 1: "we are done"
2226
if upper - lower == 2:
2227
# You could test it by checking rif_m**(lower+1), but I think
2228
# that's a waste of time since it won't be conclusive.
2229
# We must test with exact integer arithmetic which takes all
2230
# the bits of self into account.
2231
sig_off()
2232
if self >= m**(lower+1):
2233
return lower + 1
2234
else:
2235
return lower
2236
elif upper - lower > 2:
2237
# this case would only happen in cases with extremely large 'self'
2238
rif_m = R(m)
2239
min_power = rif_m**lower
2240
middle = upper-lower
2241
pow_2 = 0
2242
while middle != 0:
2243
middle >>= 1
2244
pow_2 += 1
2245
# if middle was an exact power of 2, adjust down
2246
if (1 << (pow_2-1)) == upper-lower:
2247
pow_2 -= 1
2248
#print upper, lower, pow_2
2249
pow_2_things = [rif_m]*pow_2
2250
for i from 1<=i<pow_2:
2251
pow_2_things[i] = pow_2_things[i-1]**2
2252
for i from pow_2>i>=0:
2253
middle = lower + int(2)**i
2254
#print "Upper: %i; Lower: %i; Middle: %i" % (upper,lower, middle)
2255
exp = min_power*pow_2_things[i]
2256
if exp > rif_self:
2257
upper = middle
2258
elif exp < rif_self:
2259
lower = middle
2260
min_power = exp
2261
else:
2262
sig_off()
2263
if m**middle <= self:
2264
return middle
2265
else:
2266
return lower
2267
sig_off()
2268
2269
if upper == 0:
2270
raise ValueError, "The input for exact_log is too large and support is not implemented."
2271
2272
return lower
2273
2274
def exact_log(self, m):
2275
r"""
2276
Returns the largest integer `k` such that `m^k \leq \text{self}`,
2277
i.e., the floor of `\log_m(\text{self})`.
2278
2279
This is guaranteed to return the correct answer even when the usual
2280
log function doesn't have sufficient precision.
2281
2282
INPUT:
2283
2284
- ``m`` - integer >= 2
2285
2286
AUTHORS:
2287
2288
- David Harvey (2006-09-15)
2289
- Joel B. Mohler (2009-04-08) -- rewrote this to handle small cases
2290
and/or easy cases up to 100x faster..
2291
2292
EXAMPLES::
2293
2294
sage: Integer(125).exact_log(5)
2295
3
2296
sage: Integer(124).exact_log(5)
2297
2
2298
sage: Integer(126).exact_log(5)
2299
3
2300
sage: Integer(3).exact_log(5)
2301
0
2302
sage: Integer(1).exact_log(5)
2303
0
2304
sage: Integer(178^1700).exact_log(178)
2305
1700
2306
sage: Integer(178^1700-1).exact_log(178)
2307
1699
2308
sage: Integer(178^1700+1).exact_log(178)
2309
1700
2310
sage: # we need to exercise the large base code path too
2311
sage: Integer(1780^1700-1).exact_log(1780)
2312
1699
2313
2314
sage: # The following are very very fast.
2315
sage: # Note that for base m a perfect power of 2, we get the exact log by counting bits.
2316
sage: n=2983579823750185701375109835; m=32
2317
sage: n.exact_log(m)
2318
18
2319
sage: # The next is a favorite of mine. The log2 approximate is exact and immediately provable.
2320
sage: n=90153710570912709517902579010793251709257901270941709247901209742124;m=213509721309572
2321
sage: n.exact_log(m)
2322
4
2323
2324
::
2325
2326
sage: x = 3^100000
2327
sage: RR(log(RR(x), 3))
2328
100000.000000000
2329
sage: RR(log(RR(x + 100000), 3))
2330
100000.000000000
2331
2332
::
2333
2334
sage: x.exact_log(3)
2335
100000
2336
sage: (x+1).exact_log(3)
2337
100000
2338
sage: (x-1).exact_log(3)
2339
99999
2340
2341
::
2342
2343
sage: x.exact_log(2.5)
2344
Traceback (most recent call last):
2345
...
2346
TypeError: Attempt to coerce non-integral RealNumber to Integer
2347
"""
2348
cdef Integer _m
2349
cdef Integer result
2350
cdef size_t n_log2
2351
cdef size_t m_log2
2352
cdef size_t guess # this will contain the final answer
2353
cdef bint guess_filled = 0 # this variable is only used in one branch below
2354
cdef mpz_t z
2355
if PY_TYPE_CHECK(m,Integer):
2356
_m=<Integer>m
2357
else:
2358
_m=<Integer>Integer(m)
2359
2360
if mpz_sgn(self.value) <= 0 or mpz_sgn(_m.value) <= 0:
2361
raise ValueError, "both self and m must be positive"
2362
if mpz_cmp_si(_m.value,2) < 0:
2363
raise ValueError, "m must be at least 2"
2364
2365
n_log2=mpz_sizeinbase(self.value,2)-1
2366
m_log2=mpz_sizeinbase(_m.value,2)-1
2367
if mpz_divisible_2exp_p(_m.value,m_log2):
2368
# Here, m is a power of 2 and the correct answer is found
2369
# by a log 2 approximation.
2370
guess = n_log2/m_log2 # truncating division
2371
elif n_log2/(m_log2+1) == n_log2/m_log2:
2372
# In this case, we have an upper bound and lower bound which
2373
# give the same answer, thus, the correct answer.
2374
guess = n_log2/m_log2
2375
elif m_log2 < 8: # i.e. m<256
2376
# if the base m is at most 256, we can use mpz_sizeinbase
2377
# to get the following guess which is either the exact
2378
# log, or 1+ the exact log
2379
guess = mpz_sizeinbase(self.value, mpz_get_si(_m.value)) - 1
2380
2381
# we've already excluded the case when m is an exact power of 2
2382
2383
if n_log2/m_log2 > 8000:
2384
# If we have a very large number of digits, it can be a nice
2385
# shortcut to test the guess using interval arithmetic.
2386
# (suggested by David Harvey and Carl Witty)
2387
# "for randomly distributed integers, the chance of this
2388
# interval-based comparison failing is absurdly low"
2389
import real_mpfi
2390
approx_compare = real_mpfi.RIF(m)**guess
2391
if self > approx_compare:
2392
guess_filled = 1
2393
elif self < approx_compare:
2394
guess_filled = 1
2395
guess = guess - 1
2396
if not guess_filled:
2397
# At this point, either
2398
# 1) self is close enough to a perfect power of m that we
2399
# need an exact comparison, or
2400
# 2) the numbers are small enough that converting to the
2401
# interval field is more work than the exact comparison.
2402
compare = _m**guess
2403
if self < compare:
2404
guess = guess - 1
2405
elif n_log2 < 5000:
2406
# for input with small exact log, it's very fast to work in exact
2407
# integer arithmetic starting from log2 approximations
2408
guess = self._exact_log_log2_iter(_m)
2409
else:
2410
# finally, we are out of easy cases this subroutine uses interval
2411
# arithmetic to guess and check the exact log.
2412
guess = self._exact_log_mpfi_log(_m)
2413
2414
result = PY_NEW(Integer)
2415
mpz_set_ui(result.value,guess)
2416
return result
2417
2418
def log(self, m=None, prec=None):
2419
r"""
2420
Returns symbolic log by default, unless the logarithm is exact (for
2421
an integer base). When precision is given, the RealField
2422
approximation to that bit precision is used.
2423
2424
This function is provided primarily so that Sage integers may be
2425
treated in the same manner as real numbers when convenient. Direct
2426
use of exact_log is probably best for arithmetic log computation.
2427
2428
INPUT:
2429
2430
- ``m`` - default: natural log base e
2431
2432
- ``prec`` - integer (default: None): if None, returns
2433
symbolic, else to given bits of precision as in RealField
2434
2435
EXAMPLES::
2436
2437
sage: Integer(124).log(5)
2438
log(124)/log(5)
2439
sage: Integer(124).log(5,100)
2440
2.9950093311241087454822446806
2441
sage: Integer(125).log(5)
2442
3
2443
sage: Integer(125).log(5,prec=53)
2444
3.00000000000000
2445
sage: log(Integer(125))
2446
log(125)
2447
2448
For extremely large numbers, this works::
2449
2450
sage: x = 3^100000
2451
sage: log(x,3)
2452
100000
2453
2454
With the new Pynac symbolic backend, log(x) also
2455
works in a reasonable amount of time for this x::
2456
2457
sage: x = 3^100000
2458
sage: log(x)
2459
log(1334971414230...5522000001)
2460
2461
But approximations are probably more useful in this
2462
case, and work to as high a precision as we desire::
2463
2464
sage: x.log(3,53) # default precision for RealField
2465
100000.000000000
2466
sage: (x+1).log(3,53)
2467
100000.000000000
2468
sage: (x+1).log(3,1000)
2469
100000.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
2470
2471
We can use non-integer bases, with default e::
2472
2473
sage: x.log(2.5,prec=53)
2474
119897.784671579
2475
2476
We also get logarithms of negative integers, via the
2477
symbolic ring, using the branch from `-pi` to `pi`::
2478
2479
sage: log(-1)
2480
I*pi
2481
2482
The logarithm of zero is done likewise::
2483
2484
sage: log(0)
2485
-Infinity
2486
"""
2487
if mpz_sgn(self.value) <= 0:
2488
from sage.symbolic.all import SR
2489
return SR(self).log()
2490
if m <= 0 and m != None:
2491
raise ValueError, "m must be positive"
2492
if prec:
2493
from sage.rings.real_mpfr import RealField
2494
if m is None:
2495
return RealField(prec)(self).log()
2496
return RealField(prec)(self).log(m)
2497
if type(m)==Integer and type(self)==Integer and m**(self.exact_log(m))==self:
2498
return self.exact_log(m)
2499
2500
from sage.symbolic.all import SR
2501
from sage.functions.log import function_log
2502
if m is None:
2503
return function_log(self,dont_call_method_on_arg=True)
2504
return function_log(self,dont_call_method_on_arg=True)/\
2505
function_log(m,dont_call_method_on_arg=True)
2506
2507
def exp(self, prec=None):
2508
r"""
2509
Returns the exponential function of self as a real number.
2510
2511
This function is provided only so that Sage integers may be treated
2512
in the same manner as real numbers when convenient.
2513
2514
INPUT:
2515
2516
2517
- ``prec`` - integer (default: None): if None, returns
2518
symbolic, else to given bits of precision as in RealField
2519
2520
2521
EXAMPLES::
2522
2523
sage: Integer(8).exp()
2524
e^8
2525
sage: Integer(8).exp(prec=100)
2526
2980.9579870417282747435920995
2527
sage: exp(Integer(8))
2528
e^8
2529
2530
For even fairly large numbers, this may not be useful.
2531
2532
::
2533
2534
sage: y=Integer(145^145)
2535
sage: y.exp()
2536
e^25024207011349079210459585279553675697932183658421565260323592409432707306554163224876110094014450895759296242775250476115682350821522931225499163750010280453185147546962559031653355159703678703793369785727108337766011928747055351280379806937944746847277089168867282654496776717056860661614337004721164703369140625
2537
sage: y.exp(prec=53) # default RealField precision
2538
+infinity
2539
"""
2540
from sage.functions.all import exp
2541
res = exp(self, dont_call_method_on_arg=True)
2542
if prec:
2543
return res.n(prec=prec)
2544
return res
2545
2546
def prime_to_m_part(self, m):
2547
"""
2548
Returns the prime-to-m part of self, i.e., the largest divisor of
2549
self that is coprime to m.
2550
2551
INPUT:
2552
2553
- ``m`` - Integer
2554
2555
OUTPUT: Integer
2556
2557
EXAMPLES::
2558
2559
sage: z = 43434
2560
sage: z.prime_to_m_part(20)
2561
21717
2562
"""
2563
late_import()
2564
return sage.rings.arith.prime_to_m_part(self, m)
2565
2566
def prime_divisors(self):
2567
"""
2568
The prime divisors of self, sorted in increasing order. If n is
2569
negative, we do *not* include -1 among the prime divisors, since
2570
-1 is not a prime number.
2571
2572
EXAMPLES::
2573
2574
sage: a = 1; a.prime_divisors()
2575
[]
2576
sage: a = 100; a.prime_divisors()
2577
[2, 5]
2578
sage: a = -100; a.prime_divisors()
2579
[2, 5]
2580
sage: a = 2004; a.prime_divisors()
2581
[2, 3, 167]
2582
"""
2583
late_import()
2584
return sage.rings.arith.prime_divisors(self)
2585
2586
prime_factors = prime_divisors
2587
2588
def divisors(self):
2589
"""
2590
Returns a list of all positive integer divisors of the integer
2591
self.
2592
2593
EXAMPLES::
2594
2595
sage: a = -3; a.divisors()
2596
[1, 3]
2597
sage: a = 6; a.divisors()
2598
[1, 2, 3, 6]
2599
sage: a = 28; a.divisors()
2600
[1, 2, 4, 7, 14, 28]
2601
sage: a = 2^5; a.divisors()
2602
[1, 2, 4, 8, 16, 32]
2603
sage: a = 100; a.divisors()
2604
[1, 2, 4, 5, 10, 20, 25, 50, 100]
2605
sage: a = 1; a.divisors()
2606
[1]
2607
sage: a = 0; a.divisors()
2608
Traceback (most recent call last):
2609
...
2610
ValueError: n must be nonzero
2611
sage: a = 2^3 * 3^2 * 17; a.divisors()
2612
[1, 2, 3, 4, 6, 8, 9, 12, 17, 18, 24, 34, 36, 51, 68, 72, 102, 136, 153, 204, 306, 408, 612, 1224]
2613
sage: a = odd_part(factorial(31))
2614
sage: v = a.divisors(); len(v)
2615
172800
2616
sage: prod(e+1 for p,e in factor(a))
2617
172800
2618
sage: all([t.divides(a) for t in v])
2619
True
2620
2621
.. note::
2622
2623
If one first computes all the divisors and then sorts it,
2624
the sorting step can easily dominate the runtime. Note,
2625
however, that (non-negative) multiplication on the left
2626
preserves relative order. One can leverage this fact to
2627
keep the list in order as one computes it using a process
2628
similar to that of the merge sort when adding new elements.
2629
"""
2630
cdef list all, prev, sorted
2631
cdef long tip, top
2632
cdef long i, j, e, ee
2633
cdef Integer apn, p, pn, z, all_tip
2634
2635
if not self:
2636
raise ValueError, "n must be nonzero"
2637
f = self.factor()
2638
2639
# All of the declarations below are for optimizing the word-sized
2640
# case. Operations are performed in c as far as possible without
2641
# overflow before moving to python objects.
2642
cdef long long p_c, pn_c, apn_c
2643
cdef long all_len, sorted_len, prev_len
2644
cdef long long* ptr
2645
cdef long long* empty_c
2646
cdef long long* swap_tmp
2647
cdef long long* all_c
2648
cdef long long* sorted_c
2649
cdef long long* prev_c
2650
2651
cdef long divisor_count = 1
2652
for p,e in f: divisor_count *= (1+e)
2653
ptr = <long long*>sage_malloc(sizeof(long long) * 3 * divisor_count)
2654
if ptr == NULL:
2655
raise MemoryError
2656
all_c = ptr
2657
sorted_c = ptr + divisor_count
2658
prev_c = ptr + (2*divisor_count)
2659
2660
# These are used to keep track of whether or not we are able to
2661
# perform the operations in machine words. A factor of two safety
2662
# margin is added to cover any floating-point rounding issues.
2663
cdef bint fits_c = True
2664
cdef double cur_max = 1
2665
cdef double fits_max = 2.0**(sizeof(long long)*8-2)
2666
2667
sorted_c[0] = 1
2668
sorted_len = 1
2669
2670
for p, e in f:
2671
2672
cur_max *= (<double>p)**e
2673
if fits_c and cur_max > fits_max:
2674
sorted = []
2675
for i from 0 <= i < sorted_len:
2676
z = <Integer>PY_NEW(Integer)
2677
mpz_set_longlong(z.value, sorted_c[i])
2678
sorted.append(z)
2679
sage_free(ptr)
2680
fits_c = False
2681
2682
# The two cases below are essentially the same algorithm, one
2683
# operating on Integers in Python lists, the other on long longs.
2684
if fits_c:
2685
2686
pn_c = p_c = p
2687
2688
swap_tmp = sorted_c
2689
sorted_c = prev_c
2690
prev_c = swap_tmp
2691
prev_len = sorted_len
2692
sorted_len = 0
2693
2694
tip = 0
2695
prev_c[prev_len] = prev_c[prev_len-1] * pn_c
2696
for i from 0 <= i < prev_len:
2697
apn_c = prev_c[i] * pn_c
2698
while prev_c[tip] < apn_c:
2699
sorted_c[sorted_len] = prev_c[tip]
2700
sorted_len += 1
2701
tip += 1
2702
sorted_c[sorted_len] = apn_c
2703
sorted_len += 1
2704
2705
for ee in range(1, e):
2706
2707
swap_tmp = all_c
2708
all_c = sorted_c
2709
sorted_c = swap_tmp
2710
all_len = sorted_len
2711
sorted_len = 0
2712
2713
pn_c *= p_c
2714
tip = 0
2715
all_c[all_len] = prev_c[prev_len-1] * pn_c
2716
for i from 0 <= i < prev_len:
2717
apn_c = prev_c[i] * pn_c
2718
while all_c[tip] < apn_c:
2719
sorted_c[sorted_len] = all_c[tip]
2720
sorted_len += 1
2721
tip += 1
2722
sorted_c[sorted_len] = apn_c
2723
sorted_len += 1
2724
2725
else:
2726
prev = sorted
2727
pn = <Integer>PY_NEW(Integer)
2728
mpz_set_ui(pn.value, 1)
2729
for ee in range(e):
2730
all = sorted
2731
sorted = []
2732
tip = 0
2733
top = len(all)
2734
mpz_mul(pn.value, pn.value, p.value) # pn *= p
2735
for a in prev:
2736
# apn = a*pn
2737
apn = <Integer>PY_NEW(Integer)
2738
mpz_mul(apn.value, (<Integer>a).value, pn.value)
2739
while tip < top:
2740
all_tip = <Integer>all[tip]
2741
if mpz_cmp(all_tip.value, apn.value) > 0:
2742
break
2743
sorted.append(all_tip)
2744
tip += 1
2745
sorted.append(apn)
2746
2747
if fits_c:
2748
# all the data is in sorted_c
2749
sorted = []
2750
for i from 0 <= i < sorted_len:
2751
z = <Integer>PY_NEW(Integer)
2752
mpz_set_longlong(z.value, sorted_c[i])
2753
sorted.append(z)
2754
sage_free(ptr)
2755
2756
return sorted
2757
2758
2759
def __pos__(self):
2760
"""
2761
EXAMPLES::
2762
2763
sage: z=43434
2764
sage: z.__pos__()
2765
43434
2766
"""
2767
return self
2768
2769
def __abs__(self):
2770
"""
2771
Computes `|self|`
2772
2773
EXAMPLES::
2774
2775
sage: z = -1
2776
sage: abs(z)
2777
1
2778
sage: abs(z) == abs(1)
2779
True
2780
"""
2781
cdef Integer x = PY_NEW(Integer)
2782
mpz_abs(x.value, self.value)
2783
return x
2784
2785
def __mod__(x, y):
2786
r"""
2787
Returns x modulo y.
2788
2789
EXAMPLES::
2790
2791
sage: z = 43
2792
sage: z % 2
2793
1
2794
sage: z % 0
2795
Traceback (most recent call last):
2796
...
2797
ZeroDivisionError: Integer modulo by zero
2798
sage: -5 % 7
2799
2
2800
sage: -5 % -7
2801
-5
2802
sage: 5 % -7
2803
-2
2804
2805
TESTS::
2806
2807
sage: signs = [(11,5), (11,-5), (-11,5), (-11,-5)]
2808
sage: control = [int(a) % int(b) for a, b in signs]
2809
sage: [a % b for a,b in signs] == control
2810
True
2811
sage: [a % int(b) for a,b in signs] == control
2812
True
2813
sage: [int(a) % b for a,b in signs] == control
2814
True
2815
2816
This example caused trouble in trac #6083::
2817
2818
sage: a = next_prime(2**31)
2819
sage: b = Integers(a)(100)
2820
sage: a % b
2821
59
2822
"""
2823
cdef Integer z = PY_NEW(Integer)
2824
cdef long yy, res
2825
2826
# first case: Integer % Integer
2827
if PY_TYPE(x) is PY_TYPE(y):
2828
if not mpz_sgn((<Integer>y).value):
2829
raise ZeroDivisionError, "Integer modulo by zero"
2830
if mpz_size((<Integer>x).value) > 100000:
2831
sig_on()
2832
mpz_fdiv_r(z.value, (<Integer>x).value, (<Integer>y).value)
2833
sig_off()
2834
else:
2835
mpz_fdiv_r(z.value, (<Integer>x).value, (<Integer>y).value)
2836
return z
2837
2838
# next: Integer % python int
2839
elif PyInt_CheckExact(y):
2840
yy = PyInt_AS_LONG(y)
2841
if yy > 0:
2842
mpz_fdiv_r_ui(z.value, (<Integer>x).value, yy)
2843
elif yy == 0:
2844
raise ZeroDivisionError, "Integer modulo by zero"
2845
else:
2846
res = mpz_fdiv_r_ui(z.value, (<Integer>x).value, -yy)
2847
if res:
2848
mpz_sub_ui(z.value, z.value, -yy)
2849
return z
2850
2851
# all other cases
2852
else:
2853
try:
2854
# we explicitly try coercing both to ZZ here to
2855
# avoid infinite loops in some cases (such as
2856
# Integers and Integers(n)), see trac #6083
2857
x = integer(x)
2858
y = integer(y)
2859
return x % y
2860
except ValueError:
2861
return bin_op(x, y, operator.mod)
2862
2863
def quo_rem(Integer self, other):
2864
"""
2865
Returns the quotient and the remainder of self divided by other.
2866
Note that the remainder returned is always either zero or of the
2867
same sign as other.
2868
2869
INPUT:
2870
2871
- ``other`` - the divisor
2872
2873
OUTPUT:
2874
2875
- ``q`` - the quotient of self/other
2876
2877
- ``r`` - the remainder of self/other
2878
2879
EXAMPLES::
2880
2881
sage: z = Integer(231)
2882
sage: z.quo_rem(2)
2883
(115, 1)
2884
sage: z.quo_rem(-2)
2885
(-116, -1)
2886
sage: z.quo_rem(0)
2887
Traceback (most recent call last):
2888
...
2889
ZeroDivisionError: Integer division by zero
2890
2891
sage: a = ZZ.random_element(10**50)
2892
sage: b = ZZ.random_element(10**15)
2893
sage: q, r = a.quo_rem(b)
2894
sage: q*b + r == a
2895
True
2896
2897
sage: 3.quo_rem(ZZ['x'].0)
2898
(0, 3)
2899
2900
TESTS:
2901
2902
The divisor can be rational as well, although the remainder
2903
will always be zero (trac #7965)::
2904
2905
sage: 5.quo_rem(QQ(2))
2906
(5/2, 0)
2907
sage: 5.quo_rem(2/3)
2908
(15/2, 0)
2909
2910
"""
2911
cdef Integer q = PY_NEW(Integer)
2912
cdef Integer r = PY_NEW(Integer)
2913
cdef long d, res
2914
2915
if PyInt_CheckExact(other):
2916
d = PyInt_AS_LONG(other)
2917
if d > 0:
2918
mpz_fdiv_qr_ui(q.value, r.value, self.value, d)
2919
elif d == 0:
2920
raise ZeroDivisionError, "Integer division by zero"
2921
else:
2922
res = mpz_fdiv_qr_ui(q.value, r.value, self.value, -d)
2923
mpz_neg(q.value, q.value)
2924
if res:
2925
mpz_sub_ui(q.value, q.value, 1)
2926
mpz_sub_ui(r.value, r.value, -d)
2927
2928
elif PY_TYPE_CHECK_EXACT(other, Integer):
2929
if mpz_sgn((<Integer>other).value) == 0:
2930
raise ZeroDivisionError, "Integer division by zero"
2931
if mpz_size((<Integer>x).value) > 100000:
2932
sig_on()
2933
mpz_fdiv_qr(q.value, r.value, self.value, (<Integer>other).value)
2934
sig_off()
2935
else:
2936
mpz_fdiv_qr(q.value, r.value, self.value, (<Integer>other).value)
2937
2938
else:
2939
left, right = canonical_coercion(self, other)
2940
return left.quo_rem(right)
2941
2942
return q, r
2943
2944
def powermod(self, exp, mod):
2945
"""
2946
Compute self\*\*exp modulo mod.
2947
2948
EXAMPLES::
2949
2950
sage: z = 2
2951
sage: z.powermod(31,31)
2952
2
2953
sage: z.powermod(0,31)
2954
1
2955
sage: z.powermod(-31,31) == 2^-31 % 31
2956
True
2957
2958
As expected, the following is invalid::
2959
2960
sage: z.powermod(31,0)
2961
Traceback (most recent call last):
2962
...
2963
ZeroDivisionError: cannot raise to a power modulo 0
2964
"""
2965
cdef Integer x, _exp, _mod
2966
_exp = Integer(exp); _mod = Integer(mod)
2967
if mpz_cmp_si(_mod.value,0) == 0:
2968
raise ZeroDivisionError, "cannot raise to a power modulo 0"
2969
2970
x = PY_NEW(Integer)
2971
2972
sig_on()
2973
mpz_powm(x.value, self.value, _exp.value, _mod.value)
2974
sig_off()
2975
2976
return x
2977
2978
def rational_reconstruction(self, Integer m):
2979
"""
2980
Return the rational reconstruction of this integer modulo m, i.e.,
2981
the unique (if it exists) rational number that reduces to self
2982
modulo m and whose numerator and denominator is bounded by
2983
sqrt(m/2).
2984
2985
EXAMPLES::
2986
2987
sage: (3/7)%100
2988
29
2989
sage: (29).rational_reconstruction(100)
2990
3/7
2991
2992
TEST:
2993
2994
Check that ticket #9345 is fixed::
2995
2996
sage: ZZ(1).rational_reconstruction(0)
2997
Traceback (most recent call last):
2998
...
2999
ZeroDivisionError: The modulus cannot be zero
3000
sage: m = ZZ.random_element(-10^6,10^6)
3001
sage: m.rational_reconstruction(0)
3002
Traceback (most recent call last):
3003
...
3004
ZeroDivisionError: The modulus cannot be zero
3005
"""
3006
import rational
3007
return rational.pyrex_rational_reconstruction(self, m)
3008
3009
def powermodm_ui(self, exp, mod):
3010
r"""
3011
Computes self\*\*exp modulo mod, where exp is an unsigned long
3012
integer.
3013
3014
EXAMPLES::
3015
3016
sage: z = 32
3017
sage: z.powermodm_ui(2, 4)
3018
0
3019
sage: z.powermodm_ui(2, 14)
3020
2
3021
sage: z.powermodm_ui(2^32-2, 14)
3022
2
3023
sage: z.powermodm_ui(2^32-1, 14)
3024
Traceback (most recent call last): # 32-bit
3025
... # 32-bit
3026
OverflowError: exp (=4294967295) must be <= 4294967294 # 32-bit
3027
8 # 64-bit
3028
sage: z.powermodm_ui(2^65, 14)
3029
Traceback (most recent call last):
3030
...
3031
OverflowError: exp (=36893488147419103232) must be <= 4294967294 # 32-bit
3032
OverflowError: exp (=36893488147419103232) must be <= 18446744073709551614 # 64-bit
3033
"""
3034
if exp < 0:
3035
raise ValueError, "exp (=%s) must be nonnegative"%exp
3036
elif exp > MAX_UNSIGNED_LONG:
3037
raise OverflowError, "exp (=%s) must be <= %s"%(exp, MAX_UNSIGNED_LONG)
3038
cdef Integer x, _mod
3039
_mod = Integer(mod)
3040
x = PY_NEW(Integer)
3041
3042
sig_on()
3043
mpz_powm_ui(x.value, self.value, exp, _mod.value)
3044
sig_off()
3045
3046
return x
3047
3048
def __int__(self):
3049
"""
3050
Return the Python int (or long) corresponding to this Sage
3051
integer.
3052
3053
EXAMPLES::
3054
3055
sage: n = 920938; n
3056
920938
3057
sage: int(n)
3058
920938
3059
sage: type(n.__int__())
3060
<type 'int'>
3061
sage: n = 99028390823409823904823098490238409823490820938; n
3062
99028390823409823904823098490238409823490820938
3063
sage: type(n.__int__())
3064
<type 'long'>
3065
"""
3066
# TODO -- this crashes on sage.math, since it is evidently written incorrectly.
3067
return mpz_get_pyintlong(self.value)
3068
#return int(mpz_get_pylong(self.value))
3069
3070
def __long__(self):
3071
"""
3072
Return long integer corresponding to this Sage integer.
3073
3074
EXAMPLES::
3075
3076
sage: n = 9023408290348092849023849820934820938490234290; n
3077
9023408290348092849023849820934820938490234290
3078
sage: long(n)
3079
9023408290348092849023849820934820938490234290L
3080
sage: n = 920938; n
3081
920938
3082
sage: long(n)
3083
920938L
3084
sage: n.__long__()
3085
920938L
3086
"""
3087
return mpz_get_pylong(self.value)
3088
3089
def __float__(self):
3090
"""
3091
Return double precision floating point representation of this
3092
integer.
3093
3094
EXAMPLES::
3095
3096
sage: n = Integer(17); float(n)
3097
17.0
3098
sage: n = Integer(902834098234908209348209834092834098); float(n)
3099
9.028340982349081e+35
3100
sage: n = Integer(-57); float(n)
3101
-57.0
3102
sage: n.__float__()
3103
-57.0
3104
sage: type(n.__float__())
3105
<type 'float'>
3106
"""
3107
return mpz_get_d(self.value)
3108
3109
def _rpy_(self):
3110
"""
3111
Returns int(self) so that rpy can convert self into an object it
3112
knows how to work with.
3113
3114
EXAMPLES::
3115
3116
sage: n = 100
3117
sage: n._rpy_()
3118
100
3119
sage: type(n._rpy_())
3120
<type 'int'>
3121
"""
3122
return self.__int__()
3123
3124
def __hash__(self):
3125
"""
3126
Return the hash of this integer.
3127
3128
This agrees with the Python hash of the corresponding Python int or
3129
long.
3130
3131
EXAMPLES::
3132
3133
sage: n = -920384; n.__hash__()
3134
-920384
3135
sage: hash(int(n))
3136
-920384
3137
sage: n = -920390823904823094890238490238484; n.__hash__()
3138
-873977844 # 32-bit
3139
6874330978542788722 # 64-bit
3140
sage: hash(long(n))
3141
-873977844 # 32-bit
3142
6874330978542788722 # 64-bit
3143
3144
TESTS::
3145
3146
sage: hash(0)
3147
0
3148
sage: hash(-1)
3149
-2
3150
sage: n = 2^31 + 2^63 + 2^95 + 2^127 + 2^128*(2^32-2)
3151
sage: hash(n) == hash(long(n))
3152
True
3153
sage: hash(n-1) == hash(long(n-1))
3154
True
3155
sage: hash(-n) == hash(long(-n))
3156
True
3157
sage: hash(1-n) == hash(long(1-n))
3158
True
3159
sage: n = 2^63 + 2^127 + 2^191 + 2^255 + 2^256*(2^64-2)
3160
sage: hash(n) == hash(long(n))
3161
True
3162
sage: hash(n-1) == hash(long(n-1))
3163
True
3164
sage: hash(-n) == hash(long(-n))
3165
True
3166
sage: hash(1-n) == hash(long(1-n))
3167
True
3168
3169
These tests come from Trac #4957:
3170
3171
sage: n = 2^31 + 2^13
3172
sage: hash(n)
3173
-2147475456 # 32-bit
3174
2147491840 # 64-bit
3175
sage: hash(n) == hash(int(n))
3176
True
3177
sage: n = 2^63 + 2^13
3178
sage: hash(n)
3179
-2147475456 # 32-bit
3180
-9223372036854767616 # 64-bit
3181
sage: hash(n) == hash(int(n))
3182
True
3183
"""
3184
return mpz_pythonhash(self.value)
3185
3186
cdef hash_c(self):
3187
"""
3188
A C version of the __hash__ function.
3189
"""
3190
return mpz_pythonhash(self.value)
3191
3192
def trial_division(self, long bound=LONG_MAX):
3193
"""
3194
Return smallest prime divisor of self up to bound,
3195
or abs(self) if no such divisor is found.
3196
3197
INPUT:
3198
3199
- ``bound`` -- a positive integer that fits in a C signed long
3200
3201
OUTPUT:
3202
3203
- a positive integer
3204
3205
EXAMPLES::
3206
3207
sage: n = next_prime(10^6)*next_prime(10^7); n.trial_division()
3208
1000003
3209
sage: (-n).trial_division()
3210
1000003
3211
sage: n.trial_division(bound=100)
3212
10000049000057
3213
sage: n.trial_division(bound=-10)
3214
Traceback (most recent call last):
3215
...
3216
ValueError: bound must be positive
3217
sage: n.trial_division(bound=0)
3218
Traceback (most recent call last):
3219
...
3220
ValueError: bound must be positive
3221
sage: ZZ(0).trial_division()
3222
Traceback (most recent call last):
3223
...
3224
ValueError: self must be nonzero
3225
3226
sage: n = next_prime(10^5) * next_prime(10^40); n.trial_division()
3227
100003
3228
sage: n.trial_division(bound=10^4)
3229
1000030000000000000000000000000000000012100363
3230
sage: (-n).trial_division(bound=10^4)
3231
1000030000000000000000000000000000000012100363
3232
sage: (-n).trial_division()
3233
100003
3234
sage: n = 2 * next_prime(10^40); n.trial_division()
3235
2
3236
sage: n = 3 * next_prime(10^40); n.trial_division()
3237
3
3238
sage: n = 5 * next_prime(10^40); n.trial_division()
3239
5
3240
sage: n = 2 * next_prime(10^4); n.trial_division()
3241
2
3242
sage: n = 3 * next_prime(10^4); n.trial_division()
3243
3
3244
sage: n = 5 * next_prime(10^4); n.trial_division()
3245
5
3246
"""
3247
if bound <= 0:
3248
raise ValueError, "bound must be positive"
3249
if mpz_sgn(self.value) == 0:
3250
raise ValueError, "self must be nonzero"
3251
cdef unsigned long n, m=7, i=1, limit, dif[8]
3252
dif[0]=6;dif[1]=4;dif[2]=2;dif[3]=4;dif[4]=2;dif[5]=4;dif[6]=6;dif[7]=2
3253
cdef Integer x = PY_NEW(Integer)
3254
if mpz_fits_ulong_p(self.value):
3255
n = mpz_get_ui(self.value) # ignores the sign automatically
3256
if n == 1: return one
3257
if n%2==0:
3258
mpz_set_ui(x.value,2); return x
3259
if n%3==0:
3260
mpz_set_ui(x.value,3); return x
3261
if n%5==0:
3262
mpz_set_ui(x.value,5); return x
3263
3264
limit = <unsigned long> sqrt_double(<double> n)
3265
if bound < limit: limit = bound
3266
# Algorithm: only trial divide by numbers that
3267
# are congruent to 1,7,11,13,17,19,23,29 mod 30=2*3*5.
3268
while m <= limit:
3269
if n%m == 0:
3270
mpz_set_ui(x.value, m); return x
3271
m += dif[i%8]
3272
i += 1
3273
mpz_abs(x.value, self.value)
3274
return x
3275
else:
3276
# self is big -- it doesn't fit in unsigned long.
3277
if mpz_even_p(self.value):
3278
mpz_set_ui(x.value,2); return x
3279
if mpz_divisible_ui_p(self.value,3):
3280
mpz_set_ui(x.value,3); return x
3281
if mpz_divisible_ui_p(self.value,5):
3282
mpz_set_ui(x.value,5); return x
3283
3284
# x.value = floor(sqrt(self.value))
3285
sig_on()
3286
mpz_abs(x.value, self.value)
3287
mpz_sqrt(x.value, x.value)
3288
if mpz_cmp_si(x.value, bound) < 0:
3289
limit = mpz_get_ui(x.value)
3290
else:
3291
limit = bound
3292
while m <= limit:
3293
if mpz_divisible_ui_p(self.value, m):
3294
mpz_set_ui(x.value, m)
3295
sig_off()
3296
return x
3297
m += dif[i%8]
3298
i += 1
3299
mpz_abs(x.value, self.value)
3300
sig_off()
3301
return x
3302
3303
def factor(self, algorithm='pari', proof=None, limit=None, int_=False,
3304
verbose=0):
3305
"""
3306
Return the prime factorization of this integer as a
3307
formal Factorization object.
3308
3309
INPUT:
3310
3311
- ``algorithm`` - string
3312
3313
- ``'pari'`` - (default) use the PARI library
3314
3315
- ``'kash'`` - use the KASH computer algebra system (requires
3316
the optional kash package)
3317
3318
- ``proof`` - bool (default: True) whether or not to
3319
prove primality of each factor (only applicable for PARI).
3320
3321
- ``limit`` - int or None (default: None) if limit is
3322
given it must fit in a signed int, and the factorization is done
3323
using trial division and primes up to limit.
3324
3325
OUTPUT:
3326
3327
- a Factorization object containing the prime factors and
3328
their multiplicities
3329
3330
EXAMPLES::
3331
3332
sage: n = 2^100 - 1; n.factor()
3333
3 * 5^3 * 11 * 31 * 41 * 101 * 251 * 601 * 1801 * 4051 * 8101 * 268501
3334
3335
This factorization can be converted into a list of pairs `(p,
3336
e)`, where `p` is prime and `e` is a positive integer. Each
3337
pair can also be accessed directly by its index (ordered by
3338
increasing size of the prime)::
3339
3340
sage: f = 60.factor()
3341
sage: list(f)
3342
[(2, 2), (3, 1), (5, 1)]
3343
sage: f[2]
3344
(5, 1)
3345
3346
Similarly, the factorization can be converted to a dictionary
3347
so the exponent can be extracted for each prime::
3348
3349
sage: f = (3^6).factor()
3350
sage: dict(f)
3351
{3: 6}
3352
sage: dict(f)[3]
3353
6
3354
3355
We use proof=False, which doesn't prove correctness of the primes
3356
that appear in the factorization::
3357
3358
sage: n = 920384092842390423848290348203948092384082349082
3359
sage: n.factor(proof=False)
3360
2 * 11 * 1531 * 4402903 * 10023679 * 619162955472170540533894518173
3361
sage: n.factor(proof=True)
3362
2 * 11 * 1531 * 4402903 * 10023679 * 619162955472170540533894518173
3363
3364
We factor using trial division only::
3365
3366
sage: n.factor(limit=1000)
3367
2 * 11 * 41835640583745019265831379463815822381094652231
3368
3369
TESTS::
3370
3371
sage: n.factor(algorithm='foobar')
3372
Traceback (most recent call last):
3373
...
3374
ValueError: Algorithm is not known
3375
"""
3376
from sage.structure.factorization import Factorization
3377
from sage.structure.factorization_integer import IntegerFactorization
3378
3379
cdef Integer n, p, unit
3380
cdef int i
3381
cdef factor_t f
3382
3383
if mpz_sgn(self.value) == 0:
3384
raise ArithmeticError, "Prime factorization of 0 not defined."
3385
3386
if mpz_sgn(self.value) > 0:
3387
n = self
3388
unit = ONE
3389
else:
3390
n = PY_NEW(Integer)
3391
unit = PY_NEW(Integer)
3392
mpz_neg(n.value, self.value)
3393
mpz_set_si(unit.value, -1)
3394
3395
if mpz_cmpabs_ui(n.value, 1) == 0:
3396
return IntegerFactorization([], unit=unit, unsafe=True,
3397
sort=False, simplify=False)
3398
3399
if limit is not None:
3400
from sage.rings.factorint import factor_trial_division
3401
return factor_trial_division(self, limit)
3402
3403
if mpz_fits_slong_p(n.value):
3404
if proof is None:
3405
from sage.structure.proof.proof import get_flag
3406
proof = get_flag(proof, "arithmetic")
3407
z_factor(&f, mpz_get_ui(n.value), proof)
3408
F = [(Integer(f.p[i]), int(f.exp[i])) for i from 0 <= i < f.num]
3409
F.sort()
3410
return IntegerFactorization(F, unit=unit, unsafe=True,
3411
sort=False, simplify=False)
3412
3413
if mpz_sizeinbase(n.value, 2) < 40:
3414
from sage.rings.factorint import factor_trial_division
3415
return factor_trial_division(self)
3416
3417
if algorithm == 'pari':
3418
from sage.rings.factorint import factor_using_pari
3419
F = factor_using_pari(n, int_=int_, debug_level=verbose, proof=proof)
3420
F.sort()
3421
return IntegerFactorization(F, unit=unit, unsafe=True,
3422
sort=False, simplify=False)
3423
elif algorithm in ['kash', 'magma']:
3424
if algorithm == 'kash':
3425
from sage.interfaces.all import kash as I
3426
else:
3427
from sage.interfaces.all import magma as I
3428
F = I.eval('Factorization(%s)'%n)
3429
i = F.rfind(']') + 1
3430
F = F[:i]
3431
F = F.replace("<","(").replace(">",")")
3432
F = eval(F)
3433
if not int_:
3434
F = [(Integer(a), int(b)) for a,b in F]
3435
return Factorization(F, unit)
3436
else:
3437
raise ValueError, "Algorithm is not known"
3438
3439
def support(self):
3440
"""
3441
Return a sorted list of the primes dividing this integer.
3442
3443
OUTPUT: The sorted list of primes appearing in the factorization of
3444
this rational with positive exponent.
3445
3446
EXAMPLES::
3447
3448
sage: factorial(10).support()
3449
[2, 3, 5, 7]
3450
sage: (-999).support()
3451
[3, 37]
3452
3453
Trying to find the support of 0 gives an arithmetic error::
3454
3455
sage: 0.support()
3456
Traceback (most recent call last):
3457
...
3458
ArithmeticError: Support of 0 not defined.
3459
"""
3460
if self.is_zero():
3461
raise ArithmeticError, "Support of 0 not defined."
3462
return sage.rings.arith.prime_factors(self)
3463
3464
def coprime_integers(self, m):
3465
"""
3466
Return the positive integers `< m` that are coprime to
3467
self.
3468
3469
EXAMPLES::
3470
3471
sage: n = 8
3472
sage: n.coprime_integers(8)
3473
[1, 3, 5, 7]
3474
sage: n.coprime_integers(11)
3475
[1, 3, 5, 7, 9]
3476
sage: n = 5; n.coprime_integers(10)
3477
[1, 2, 3, 4, 6, 7, 8, 9]
3478
sage: n.coprime_integers(5)
3479
[1, 2, 3, 4]
3480
sage: n = 99; n.coprime_integers(99)
3481
[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]
3482
3483
AUTHORS:
3484
3485
- Naqi Jaffery (2006-01-24): examples
3486
3487
ALGORITHM: Naive - compute lots of GCD's. If this isn't good enough
3488
for you, please code something better and submit a patch.
3489
"""
3490
# TODO -- make VASTLY faster
3491
v = []
3492
for n in range(1,m):
3493
if self.gcd(n) == 1:
3494
v.append(Integer(n))
3495
return v
3496
3497
def divides(self, n):
3498
"""
3499
Return True if self divides n.
3500
3501
EXAMPLES::
3502
3503
sage: Z = IntegerRing()
3504
sage: Z(5).divides(Z(10))
3505
True
3506
sage: Z(0).divides(Z(5))
3507
False
3508
sage: Z(10).divides(Z(5))
3509
False
3510
"""
3511
cdef bint t
3512
cdef Integer _n
3513
_n = Integer(n)
3514
if mpz_sgn(self.value) == 0:
3515
return mpz_sgn(_n.value) == 0
3516
sig_on()
3517
t = mpz_divisible_p(_n.value, self.value)
3518
sig_off()
3519
return t
3520
3521
cpdef RingElement _valuation(Integer self, Integer p):
3522
r"""
3523
Return the p-adic valuation of self.
3524
3525
We do not require that p be prime, but it must be at least 2. For
3526
more documentation see ``valuation``
3527
3528
AUTHORS:
3529
3530
- David Roe (3/31/07)
3531
"""
3532
if mpz_sgn(self.value) == 0:
3533
return sage.rings.infinity.infinity
3534
if mpz_cmp_ui(p.value, 2) < 0:
3535
raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
3536
3537
cdef Integer v = PY_NEW(Integer)
3538
cdef mpz_t u
3539
mpz_init(u)
3540
sig_on()
3541
mpz_set_ui(v.value, mpz_remove(u, self.value, p.value))
3542
sig_off()
3543
mpz_clear(u)
3544
return v
3545
3546
cdef object _val_unit(Integer self, Integer p):
3547
r"""
3548
Returns a pair: the p-adic valuation of self, and the p-adic unit
3549
of self.
3550
3551
We do not require the p be prime, but it must be at least 2. For
3552
more documentation see ``val_unit``
3553
3554
AUTHORS:
3555
3556
- David Roe (2007-03-31)
3557
"""
3558
cdef Integer v, u
3559
if mpz_cmp_ui(p.value, 2) < 0:
3560
raise ValueError, "You can only compute the valuation with respect to a integer larger than 1."
3561
if self == 0:
3562
u = ONE
3563
Py_INCREF(ONE)
3564
return (sage.rings.infinity.infinity, u)
3565
v = PY_NEW(Integer)
3566
u = PY_NEW(Integer)
3567
sig_on()
3568
mpz_set_ui(v.value, mpz_remove(u.value, self.value, p.value))
3569
sig_off()
3570
return (v, u)
3571
3572
def valuation(self, p):
3573
"""
3574
Return the p-adic valuation of self.
3575
3576
INPUT:
3577
3578
- ``p`` - an integer at least 2.
3579
3580
EXAMPLE::
3581
3582
sage: n = 60
3583
sage: n.valuation(2)
3584
2
3585
sage: n.valuation(3)
3586
1
3587
sage: n.valuation(7)
3588
0
3589
sage: n.valuation(1)
3590
Traceback (most recent call last):
3591
...
3592
ValueError: You can only compute the valuation with respect to a integer larger than 1.
3593
3594
We do not require that p is a prime::
3595
3596
sage: (2^11).valuation(4)
3597
5
3598
"""
3599
return self._valuation(Integer(p))
3600
3601
# Alias for valuation
3602
ord = valuation
3603
3604
def val_unit(self, p):
3605
r"""
3606
Returns a pair: the p-adic valuation of self, and the p-adic unit
3607
of self.
3608
3609
INPUT:
3610
3611
- ``p`` - an integer at least 2.
3612
3613
OUTPUT:
3614
3615
- ``v_p(self)`` - the p-adic valuation of ``self``
3616
3617
- ``u_p(self)`` - ``self`` / `p^{v_p(\mathrm{self})}`
3618
3619
EXAMPLE::
3620
3621
sage: n = 60
3622
sage: n.val_unit(2)
3623
(2, 15)
3624
sage: n.val_unit(3)
3625
(1, 20)
3626
sage: n.val_unit(7)
3627
(0, 60)
3628
sage: (2^11).val_unit(4)
3629
(5, 2)
3630
sage: 0.val_unit(2)
3631
(+Infinity, 1)
3632
"""
3633
return self._val_unit(Integer(p))
3634
3635
def odd_part(self):
3636
r"""
3637
The odd part of the integer `n`. This is `n / 2^v`,
3638
where `v = \mathrm{valuation}(n,2)`.
3639
3640
IMPLEMENTATION:
3641
3642
Currently returns 0 when self is 0. This behaviour is fairly arbitrary,
3643
and in Sage 4.6 this special case was not handled at all, eventually
3644
propagating a TypeError. The caller should not rely on the behaviour
3645
in case self is 0.
3646
3647
EXAMPLES::
3648
3649
sage: odd_part(5)
3650
5
3651
sage: odd_part(4)
3652
1
3653
sage: odd_part(factorial(31))
3654
122529844256906551386796875
3655
"""
3656
cdef Integer odd
3657
cdef unsigned long bits
3658
3659
if mpz_cmpabs_ui(self.value, 1) <= 0:
3660
return self
3661
3662
odd = PY_NEW(Integer)
3663
bits = mpz_scan1(self.value, 0)
3664
mpz_tdiv_q_2exp(odd.value, self.value, bits)
3665
return odd
3666
3667
cdef Integer _divide_knowing_divisible_by(Integer self, Integer right):
3668
r"""
3669
Returns the integer self / right when self is divisible by right.
3670
3671
If self is not divisible by right, the return value is undefined,
3672
and may not even be close to self/right. For more documentation see
3673
``divide_knowing_divisible_by``
3674
3675
AUTHORS:
3676
3677
- David Roe (2007-03-31)
3678
"""
3679
if mpz_cmp_ui(right.value, 0) == 0:
3680
raise ZeroDivisionError
3681
cdef Integer x
3682
x = PY_NEW(Integer)
3683
if mpz_size(self.value) + mpz_size((<Integer>right).value) > 100000:
3684
# Only use the signal handler (to enable ctrl-c out) when the
3685
# quotient might take a while to compute
3686
sig_on()
3687
mpz_divexact(x.value, self.value, right.value)
3688
sig_off()
3689
else:
3690
mpz_divexact(x.value, self.value, right.value)
3691
return x
3692
3693
def divide_knowing_divisible_by(self, right):
3694
r"""
3695
Returns the integer self / right when self is divisible by right.
3696
3697
If self is not divisible by right, the return value is undefined,
3698
and may not even be close to self/right for multi-word integers.
3699
3700
EXAMPLES::
3701
3702
sage: a = 8; b = 4
3703
sage: a.divide_knowing_divisible_by(b)
3704
2
3705
sage: (100000).divide_knowing_divisible_by(25)
3706
4000
3707
sage: (100000).divide_knowing_divisible_by(26) # close (random)
3708
3846
3709
3710
However, often it's way off.
3711
3712
::
3713
3714
sage: a = 2^70; a
3715
1180591620717411303424
3716
sage: a // 11 # floor divide
3717
107326510974310118493
3718
sage: a.divide_knowing_divisible_by(11) # way off and possibly random
3719
43215361478743422388970455040
3720
"""
3721
return self._divide_knowing_divisible_by(right)
3722
3723
def _lcm(self, Integer n):
3724
"""
3725
Returns the least common multiple of self and `n`.
3726
3727
EXAMPLES::
3728
3729
sage: n = 60
3730
sage: n._lcm(150)
3731
300
3732
"""
3733
cdef Integer z = PY_NEW(Integer)
3734
sig_on()
3735
mpz_lcm(z.value, self.value, n.value)
3736
sig_off()
3737
return z
3738
3739
def denominator(self):
3740
"""
3741
Return the denominator of this integer.
3742
3743
EXAMPLES::
3744
3745
sage: x = 5
3746
sage: x.denominator()
3747
1
3748
sage: x = 0
3749
sage: x.denominator()
3750
1
3751
"""
3752
return ONE
3753
3754
def numerator(self):
3755
"""
3756
Return the numerator of this integer.
3757
3758
EXAMPLES::
3759
3760
sage: x = 5
3761
sage: x.numerator()
3762
5
3763
3764
::
3765
3766
sage: x = 0
3767
sage: x.numerator()
3768
0
3769
"""
3770
return self
3771
3772
def factorial(self):
3773
r"""
3774
Return the factorial `n! = 1 \cdot 2 \cdot 3 \cdots n`.
3775
3776
If the input does not fit in an ``unsigned long int`` a symbolic
3777
expression is returned.
3778
3779
EXAMPLES::
3780
3781
sage: for n in srange(7):
3782
... print n, n.factorial()
3783
0 1
3784
1 1
3785
2 2
3786
3 6
3787
4 24
3788
5 120
3789
6 720
3790
sage: 234234209384023842034.factorial()
3791
factorial(234234209384023842034)
3792
"""
3793
if mpz_sgn(self.value) < 0:
3794
raise ValueError, "factorial -- self = (%s) must be nonnegative"%self
3795
3796
if not mpz_fits_uint_p(self.value):
3797
from sage.functions.all import factorial
3798
return factorial(self, hold=True)
3799
3800
cdef Integer z = PY_NEW(Integer)
3801
3802
sig_on()
3803
mpz_fac_ui(z.value, mpz_get_ui(self.value))
3804
sig_off()
3805
3806
return z
3807
3808
@cython.cdivision(True)
3809
def multifactorial(self, int k):
3810
r"""
3811
Computes the k-th factorial `n!^{(k)}` of self. For k=1
3812
this is the standard factorial, and for k greater than one it is
3813
the product of every k-th terms down from self to k. The recursive
3814
definition is used to extend this function to the negative
3815
integers.
3816
3817
EXAMPLES::
3818
3819
sage: 5.multifactorial(1)
3820
120
3821
sage: 5.multifactorial(2)
3822
15
3823
sage: 23.multifactorial(2)
3824
316234143225
3825
sage: prod([1..23, step=2])
3826
316234143225
3827
sage: (-29).multifactorial(7)
3828
1/2640
3829
"""
3830
if k <= 0:
3831
raise ValueError, "multifactorial only defined for positive values of k"
3832
3833
if not mpz_fits_sint_p(self.value):
3834
raise ValueError, "multifactorial not implemented for n >= 2^32.\nThis is probably OK, since the answer would have billions of digits."
3835
3836
cdef int n = mpz_get_si(self.value)
3837
3838
# base case
3839
if 0 < n < k:
3840
return ONE
3841
3842
# easy to calculate
3843
elif n % k == 0:
3844
factorial = Integer(n/k).factorial()
3845
if k == 2:
3846
return factorial << (n/k)
3847
else:
3848
return factorial * Integer(k)**(n/k)
3849
3850
# negative base case
3851
elif -k < n < 0:
3852
return ONE / (self+k)
3853
3854
# reflection case
3855
elif n < -k:
3856
if (n/k) % 2:
3857
sign = -ONE
3858
else:
3859
sign = ONE
3860
return sign / Integer(-k-n).multifactorial(k)
3861
3862
# compute the actual product, optimizing the number of large
3863
# multiplications
3864
cdef int i,j
3865
3866
# we need (at most) log_2(#factors) concurrent sub-products
3867
cdef int prod_count = <int>ceil_c(log_c(n/k+1)/log_c(2))
3868
cdef mpz_t* sub_prods = <mpz_t*>sage_malloc(prod_count * sizeof(mpz_t))
3869
if sub_prods == NULL:
3870
raise MemoryError
3871
for i from 0 <= i < prod_count:
3872
mpz_init(sub_prods[i])
3873
3874
sig_on()
3875
3876
cdef residue = n % k
3877
cdef int tip = 0
3878
for i from 1 <= i <= n//k:
3879
mpz_set_ui(sub_prods[tip], k*i + residue)
3880
# for the i-th terms we use the bits of i to calculate how many
3881
# times we need to multiply "up" the stack of sub-products
3882
for j from 0 <= j < 32:
3883
if i & (1 << j):
3884
break
3885
tip -= 1
3886
mpz_mul(sub_prods[tip], sub_prods[tip], sub_prods[tip+1])
3887
tip += 1
3888
cdef int last = tip-1
3889
for tip from last > tip >= 0:
3890
mpz_mul(sub_prods[tip], sub_prods[tip], sub_prods[tip+1])
3891
3892
sig_off()
3893
3894
cdef Integer z = PY_NEW(Integer)
3895
mpz_swap(z.value, sub_prods[0])
3896
3897
for i from 0 <= i < prod_count:
3898
mpz_clear(sub_prods[i])
3899
sage_free(sub_prods)
3900
3901
return z
3902
3903
def gamma(self):
3904
r"""
3905
The gamma function on integers is the factorial function (shifted by
3906
one) on positive integers, and `\pm \infty` on non-positive integers.
3907
3908
EXAMPLES::
3909
3910
sage: gamma(5)
3911
24
3912
sage: gamma(0)
3913
Infinity
3914
sage: gamma(-1)
3915
Infinity
3916
sage: gamma(-2^150)
3917
Infinity
3918
"""
3919
if mpz_sgn(self.value) > 0:
3920
return (self-ONE).factorial()
3921
else:
3922
from sage.rings.infinity import unsigned_infinity
3923
return unsigned_infinity
3924
3925
def floor(self):
3926
"""
3927
Return the floor of self, which is just self since self is an
3928
integer.
3929
3930
EXAMPLES::
3931
3932
sage: n = 6
3933
sage: n.floor()
3934
6
3935
"""
3936
return self
3937
3938
def ceil(self):
3939
"""
3940
Return the ceiling of self, which is self since self is an
3941
integer.
3942
3943
EXAMPLES::
3944
3945
sage: n = 6
3946
sage: n.ceil()
3947
6
3948
"""
3949
return self
3950
3951
def real(self):
3952
"""
3953
Returns the real part of self, which is self.
3954
3955
EXAMPLES::
3956
3957
sage: Integer(-4).real()
3958
-4
3959
"""
3960
return self
3961
3962
def imag(self):
3963
"""
3964
Returns the imaginary part of self, which is zero.
3965
3966
EXAMPLES::
3967
3968
sage: Integer(9).imag()
3969
0
3970
"""
3971
return zero
3972
3973
def is_one(self):
3974
r"""
3975
Returns ``True`` if the integer is `1`, otherwise ``False``.
3976
3977
EXAMPLES::
3978
3979
sage: Integer(1).is_one()
3980
True
3981
sage: Integer(0).is_one()
3982
False
3983
"""
3984
return mpz_cmp_si(self.value, 1) == 0
3985
3986
def __nonzero__(self):
3987
r"""
3988
Returns ``True`` if the integer is not `0`, otherwise ``False``.
3989
3990
EXAMPLES::
3991
3992
sage: Integer(1).is_zero()
3993
False
3994
sage: Integer(0).is_zero()
3995
True
3996
"""
3997
return mpz_sgn(self.value) != 0
3998
3999
def is_integral(self):
4000
"""
4001
Return ``True`` since integers are integral, i.e.,
4002
satisfy a monic polynomial with integer coefficients.
4003
4004
EXAMPLES::
4005
4006
sage: Integer(3).is_integral()
4007
True
4008
"""
4009
return True
4010
4011
def is_unit(self):
4012
r"""
4013
Returns ``true`` if this integer is a unit, i.e., 1 or `-1`.
4014
4015
EXAMPLES::
4016
4017
sage: for n in srange(-2,3):
4018
... print n, n.is_unit()
4019
-2 False
4020
-1 True
4021
0 False
4022
1 True
4023
2 False
4024
"""
4025
return mpz_cmpabs_ui(self.value, 1) == 0
4026
4027
def is_square(self):
4028
r"""
4029
Returns ``True`` if self is a perfect square.
4030
4031
EXAMPLES::
4032
4033
sage: Integer(4).is_square()
4034
True
4035
sage: Integer(41).is_square()
4036
False
4037
"""
4038
return mpz_perfect_square_p(self.value)
4039
4040
def is_power(self):
4041
r"""
4042
Returns ``True`` if self is a perfect power, i.e. if there exist
4043
integers `a` and `b`, `b > 1` with `self = a^b`.
4044
4045
EXAMPLES::
4046
4047
sage: Integer(-27).is_power()
4048
True
4049
sage: Integer(12).is_power()
4050
False
4051
"""
4052
return mpz_perfect_power_p(self.value)
4053
4054
cdef bint _is_power_of(Integer self, Integer n):
4055
r"""
4056
Returns a non-zero int if there is an integer b with
4057
`\mathtt{self} = n^b`.
4058
4059
For more documentation see ``is_power_of``.
4060
4061
AUTHORS:
4062
4063
- David Roe (2007-03-31)
4064
"""
4065
cdef int a
4066
cdef unsigned long b, c
4067
cdef mpz_t u, sabs, nabs
4068
a = mpz_cmp_ui(n.value, 2)
4069
if a <= 0: # n <= 2
4070
if a == 0: # n == 2
4071
if mpz_popcount(self.value) == 1: #number of bits set in self == 1
4072
return 1
4073
else:
4074
return 0
4075
a = mpz_cmp_si(n.value, -2)
4076
if a >= 0: # -2 <= n < 2:
4077
a = mpz_get_si(n.value)
4078
if a == 1: # n == 1
4079
if mpz_cmp_ui(self.value, 1) == 0: # Only 1 is a power of 1
4080
return 1
4081
else:
4082
return 0
4083
elif a == 0: # n == 0
4084
if mpz_cmp_ui(self.value, 0) == 0 or mpz_cmp_ui(self.value, 1) == 0: # 0^0 = 1, 0^x = 0
4085
return 1
4086
else:
4087
return 0
4088
elif a == -1: # n == -1
4089
if mpz_cmp_ui(self.value, 1) == 0 or mpz_cmp_si(self.value, -1) == 0: # 1 and -1 are powers of -1
4090
return 1
4091
else:
4092
return 0
4093
elif a == -2: # n == -2
4094
mpz_init(sabs)
4095
mpz_abs(sabs, self.value)
4096
if mpz_popcount(sabs) == 1: # number of bits set in |self| == 1
4097
b = mpz_scan1(sabs, 0) % 2 # b == 1 if |self| is an odd power of 2, 0 if |self| is an even power
4098
mpz_clear(sabs)
4099
if (b == 1 and mpz_cmp_ui(self.value, 0) < 0) or (b == 0 and mpz_cmp_ui(self.value, 0) > 0):
4100
# An odd power of -2 is negative, an even power must be positive.
4101
return 1
4102
else: # number of bits set in |self| is not 1, so self cannot be a power of -2
4103
return 0
4104
else: # |self| is not a power of 2, so self cannot be a power of -2
4105
return 0
4106
else: # n < -2
4107
mpz_init(nabs)
4108
mpz_neg(nabs, n.value)
4109
if mpz_popcount(nabs) == 1: # |n| = 2^k for k >= 2. We special case this for speed
4110
mpz_init(sabs)
4111
mpz_abs(sabs, self.value)
4112
if mpz_popcount(sabs) == 1: # |self| = 2^L for some L >= 0.
4113
b = mpz_scan1(sabs, 0) # the bit that self is set at
4114
c = mpz_scan1(nabs, 0) # the bit that n is set at
4115
# Having obtained b and c, we're done with nabs and sabs (on this branch anyway)
4116
mpz_clear(nabs)
4117
mpz_clear(sabs)
4118
if b % c == 0: # Now we know that |self| is a power of |n|
4119
b = (b // c) % 2 # Whether b // c is even or odd determines whether (-2^c)^(b // c) is positive or negative
4120
a = mpz_cmp_ui(self.value, 0)
4121
if b == 0 and a > 0 or b == 1 and a < 0:
4122
# These two cases are that b // c is even and self positive, or b // c is odd and self negative
4123
return 1
4124
else: # The sign of self is wrong
4125
return 0
4126
else: # Since |self| is not a power of |n|, self cannot be a power of n
4127
return 0
4128
else: # self is not a power of 2, and thus cannot be a power of n, which is a power of 2.
4129
mpz_clear(nabs)
4130
mpz_clear(sabs)
4131
return 0
4132
else: # |n| is not a power of 2, so we use mpz_remove
4133
mpz_init(u)
4134
sig_on()
4135
b = mpz_remove(u, self.value, nabs)
4136
sig_off()
4137
# Having obtained b and u, we're done with nabs
4138
mpz_clear(nabs)
4139
if mpz_cmp_ui(u, 1) == 0: # self is a power of |n|
4140
mpz_clear(u)
4141
if b % 2 == 0: # an even power of |n|, and since self > 0, this means that self is a power of n
4142
return 1
4143
else:
4144
return 0
4145
elif mpz_cmp_si(u, -1) == 0: # -self is a power of |n|
4146
mpz_clear(u)
4147
if b % 2 == 1: # an odd power of |n|, and thus self is a power of n
4148
return 1
4149
else:
4150
return 0
4151
else: # |self| is not a power of |n|, so self cannot be a power of n
4152
mpz_clear(u)
4153
return 0
4154
elif mpz_popcount(n.value) == 1: # n > 2 and in fact n = 2^k for k >= 2
4155
if mpz_popcount(self.value) == 1: # since n is a power of 2, so must self be.
4156
if mpz_scan1(self.value, 0) % mpz_scan1(n.value, 0) == 0: # log_2(self) is divisible by log_2(n)
4157
return 1
4158
else:
4159
return 0
4160
else: # self is not a power of 2, and thus not a power of n
4161
return 0
4162
else: # n > 2, but not a power of 2, so we use mpz_remove
4163
mpz_init(u)
4164
sig_on()
4165
mpz_remove(u, self.value, n.value)
4166
sig_off()
4167
a = mpz_cmp_ui(u, 1)
4168
mpz_clear(u)
4169
if a == 0:
4170
return 1
4171
else:
4172
return 0
4173
4174
def is_power_of(Integer self, n):
4175
r"""
4176
Returns ``True`` if there is an integer b with
4177
`\mathtt{self} = n^b`.
4178
4179
EXAMPLES::
4180
4181
sage: Integer(64).is_power_of(4)
4182
True
4183
sage: Integer(64).is_power_of(16)
4184
False
4185
4186
TESTS::
4187
4188
sage: Integer(-64).is_power_of(-4)
4189
True
4190
sage: Integer(-32).is_power_of(-2)
4191
True
4192
sage: Integer(1).is_power_of(1)
4193
True
4194
sage: Integer(-1).is_power_of(-1)
4195
True
4196
sage: Integer(0).is_power_of(1)
4197
False
4198
sage: Integer(0).is_power_of(0)
4199
True
4200
sage: Integer(1).is_power_of(0)
4201
True
4202
sage: Integer(1).is_power_of(8)
4203
True
4204
sage: Integer(-8).is_power_of(2)
4205
False
4206
sage: Integer(-81).is_power_of(-3)
4207
False
4208
4209
.. note::
4210
4211
For large integers self, is_power_of() is faster than
4212
is_power(). The following examples gives some indication of
4213
how much faster.
4214
4215
::
4216
4217
sage: b = lcm(range(1,10000))
4218
sage: b.exact_log(2)
4219
14446
4220
sage: t=cputime()
4221
sage: for a in range(2, 1000): k = b.is_power()
4222
sage: cputime(t) # random
4223
0.53203299999999976
4224
sage: t=cputime()
4225
sage: for a in range(2, 1000): k = b.is_power_of(2)
4226
sage: cputime(t) # random
4227
0.0
4228
sage: t=cputime()
4229
sage: for a in range(2, 1000): k = b.is_power_of(3)
4230
sage: cputime(t) # random
4231
0.032002000000000308
4232
4233
::
4234
4235
sage: b = lcm(range(1, 1000))
4236
sage: b.exact_log(2)
4237
1437
4238
sage: t=cputime()
4239
sage: for a in range(2, 10000): k = b.is_power() # note that we change the range from the example above
4240
sage: cputime(t) # random
4241
0.17201100000000036
4242
sage: t=cputime(); TWO=int(2)
4243
sage: for a in range(2, 10000): k = b.is_power_of(TWO)
4244
sage: cputime(t) # random
4245
0.0040000000000000036
4246
sage: t=cputime()
4247
sage: for a in range(2, 10000): k = b.is_power_of(3)
4248
sage: cputime(t) # random
4249
0.040003000000000011
4250
sage: t=cputime()
4251
sage: for a in range(2, 10000): k = b.is_power_of(a)
4252
sage: cputime(t) # random
4253
0.02800199999999986
4254
"""
4255
if not PY_TYPE_CHECK(n, Integer):
4256
n = Integer(n)
4257
return self._is_power_of(n)
4258
4259
def is_prime_power(self, flag=0):
4260
r"""
4261
Returns True if `x` is a prime power, and False otherwise.
4262
4263
INPUT:
4264
4265
- ``flag`` (for primality testing) - int
4266
- ``0`` (default): use a combination of algorithms.
4267
- ``1``: certify primality using the Pocklington-Lehmer test.
4268
- ``2``: certify primality using the APRCL test.
4269
4270
EXAMPLES::
4271
4272
sage: (-10).is_prime_power()
4273
False
4274
sage: (10).is_prime_power()
4275
False
4276
sage: (64).is_prime_power()
4277
True
4278
sage: (3^10000).is_prime_power()
4279
True
4280
sage: (10000).is_prime_power(flag=1)
4281
False
4282
4283
We check that \#4777 is fixed::
4284
4285
sage: n = 150607571^14
4286
sage: n.is_prime_power()
4287
True
4288
"""
4289
if self.is_zero():
4290
return False
4291
elif self.is_one():
4292
return True
4293
elif mpz_sgn(self.value) < 0:
4294
return False
4295
if self.is_prime():
4296
return True
4297
if not self.is_perfect_power():
4298
return False
4299
k, g = self._pari_().ispower()
4300
if not k:
4301
raise RuntimeError, "Inconsistent results between GMP and PARI"
4302
return g.isprime(flag=flag)
4303
4304
def is_prime(self, proof=None):
4305
r"""
4306
Returns ``True`` if ``self`` is prime.
4307
4308
INPUT:
4309
4310
- ``proof`` -- If False, use a strong pseudo-primality test.
4311
If True, use a provable primality test. If unset, use the
4312
default arithmetic proof flag.
4313
4314
.. note::
4315
4316
Integer primes are by definition *positive*! This is
4317
different than Magma, but the same as in PARI. See also the
4318
:meth:`is_irreducible()` method.
4319
4320
EXAMPLES::
4321
4322
sage: z = 2^31 - 1
4323
sage: z.is_prime()
4324
True
4325
sage: z = 2^31
4326
sage: z.is_prime()
4327
False
4328
sage: z = 7
4329
sage: z.is_prime()
4330
True
4331
sage: z = -7
4332
sage: z.is_prime()
4333
False
4334
sage: z.is_irreducible()
4335
True
4336
4337
::
4338
4339
sage: z = 10^80 + 129
4340
sage: z.is_prime(proof=False)
4341
True
4342
sage: z.is_prime(proof=True)
4343
True
4344
4345
IMPLEMENTATION: Calls the PARI ``isprime`` function.
4346
"""
4347
from sage.structure.proof.proof import get_flag
4348
proof = get_flag(proof, "arithmetic")
4349
if proof:
4350
return bool(self._pari_().isprime())
4351
else:
4352
return bool(self._pari_().ispseudoprime())
4353
4354
def is_irreducible(self):
4355
r"""
4356
Returns ``True`` if self is irreducible, i.e. +/-
4357
prime
4358
4359
EXAMPLES::
4360
4361
sage: z = 2^31 - 1
4362
sage: z.is_irreducible()
4363
True
4364
sage: z = 2^31
4365
sage: z.is_irreducible()
4366
False
4367
sage: z = 7
4368
sage: z.is_irreducible()
4369
True
4370
sage: z = -7
4371
sage: z.is_irreducible()
4372
True
4373
"""
4374
n = self if self >= 0 else -self
4375
return bool(n._pari_().isprime())
4376
4377
def is_pseudoprime(self):
4378
r"""
4379
Returns ``True`` if self is a pseudoprime
4380
4381
EXAMPLES::
4382
4383
sage: z = 2^31 - 1
4384
sage: z.is_pseudoprime()
4385
True
4386
sage: z = 2^31
4387
sage: z.is_pseudoprime()
4388
False
4389
"""
4390
return bool(self._pari_().ispseudoprime())
4391
4392
def is_perfect_power(self):
4393
r"""
4394
Returns ``True`` if self is a perfect power.
4395
4396
EXAMPLES::
4397
4398
sage: z = 8
4399
sage: z.is_perfect_power()
4400
True
4401
sage: 144.is_perfect_power()
4402
True
4403
sage: 10.is_perfect_power()
4404
False
4405
sage: (-8).is_perfect_power()
4406
True
4407
sage: (-4).is_perfect_power()
4408
False
4409
4410
This is a test to make sure we work around a bug in GMP,
4411
see trac \#4612.
4412
4413
::
4414
4415
sage: [ -a for a in srange(100) if not (-a^3).is_perfect_power() ]
4416
[]
4417
"""
4418
cdef mpz_t tmp
4419
cdef int res
4420
if mpz_sgn(self.value) < 0:
4421
if mpz_cmp_si(self.value, -1) == 0:
4422
return True
4423
mpz_init(tmp)
4424
mpz_neg(tmp, self.value)
4425
while mpz_perfect_square_p(tmp):
4426
mpz_sqrt(tmp, tmp)
4427
res = mpz_perfect_power_p(tmp)
4428
mpz_clear(tmp)
4429
if res:
4430
return True
4431
else:
4432
return False
4433
return mpz_perfect_power_p(self.value)
4434
4435
def is_norm(self, K, element=False, proof=True):
4436
r"""
4437
See ``QQ(self).is_norm()``.
4438
4439
EXAMPLES::
4440
4441
sage: K = NumberField(x^2 - 2, 'beta')
4442
sage: n = 4
4443
sage: n.is_norm(K)
4444
True
4445
sage: 5.is_norm(K)
4446
False
4447
sage: 7.is_norm(QQ)
4448
True
4449
sage: n.is_norm(K, element=True)
4450
(True, -4*beta + 6)
4451
sage: n.is_norm(K, element=True)[1].norm()
4452
4
4453
sage: n = 5
4454
sage: n.is_norm(K, element=True)
4455
(False, None)
4456
sage: n = 7
4457
sage: n.is_norm(QQ, element=True)
4458
(True, 7)
4459
4460
"""
4461
from sage.rings.rational_field import QQ
4462
return QQ(self).is_norm(K, element=element, proof=proof)
4463
4464
def _bnfisnorm(self, K, proof=True, extra_primes=0):
4465
r"""
4466
See ``QQ(self)._bnfisnorm()``.
4467
4468
EXAMPLES::
4469
4470
sage: 3._bnfisnorm(QuadraticField(-1, 'i'))
4471
(1, 3)
4472
sage: 7._bnfisnorm(CyclotomicField(7))
4473
(-zeta7 + 1, 1) # 64-bit
4474
(-zeta7^5 + zeta7^4, 1) # 32-bit
4475
"""
4476
from sage.rings.rational_field import QQ
4477
return QQ(self)._bnfisnorm(K, proof=proof, extra_primes=extra_primes)
4478
4479
4480
def jacobi(self, b):
4481
r"""
4482
Calculate the Jacobi symbol `\left(\frac{self}{b}\right)`.
4483
4484
EXAMPLES::
4485
4486
sage: z = -1
4487
sage: z.jacobi(17)
4488
1
4489
sage: z.jacobi(19)
4490
-1
4491
sage: z.jacobi(17*19)
4492
-1
4493
sage: (2).jacobi(17)
4494
1
4495
sage: (3).jacobi(19)
4496
-1
4497
sage: (6).jacobi(17*19)
4498
-1
4499
sage: (6).jacobi(33)
4500
0
4501
sage: a = 3; b = 7
4502
sage: a.jacobi(b) == -b.jacobi(a)
4503
True
4504
"""
4505
cdef long tmp
4506
if PY_TYPE_CHECK(b, int):
4507
tmp = b
4508
if (tmp & 1) == 0:
4509
raise ValueError, "Jacobi symbol not defined for even b."
4510
return mpz_kronecker_si(self.value, tmp)
4511
if not PY_TYPE_CHECK(b, Integer):
4512
b = Integer(b)
4513
if mpz_even_p((<Integer>b).value):
4514
raise ValueError, "Jacobi symbol not defined for even b."
4515
return mpz_jacobi(self.value, (<Integer>b).value)
4516
4517
def kronecker(self, b):
4518
r"""
4519
Calculate the Kronecker symbol `\left(\frac{self}{b}\right)`
4520
with the Kronecker extension `(self/2)=(2/self)` when `self` is odd,
4521
or `(self/2)=0` when `self` is even.
4522
4523
EXAMPLES::
4524
4525
sage: z = 5
4526
sage: z.kronecker(41)
4527
1
4528
sage: z.kronecker(43)
4529
-1
4530
sage: z.kronecker(8)
4531
-1
4532
sage: z.kronecker(15)
4533
0
4534
sage: a = 2; b = 5
4535
sage: a.kronecker(b) == b.kronecker(a)
4536
True
4537
"""
4538
if PY_TYPE_CHECK(b, int):
4539
return mpz_kronecker_si(self.value, b)
4540
if not PY_TYPE_CHECK(b, Integer):
4541
b = Integer(b)
4542
return mpz_kronecker(self.value, (<Integer>b).value)
4543
4544
def radical(self, *args, **kwds):
4545
r"""
4546
Return the product of the prime divisors of self. Computing
4547
the radical of zero gives an error.
4548
4549
EXAMPLES::
4550
4551
sage: Integer(10).radical()
4552
10
4553
sage: Integer(20).radical()
4554
10
4555
sage: Integer(-100).radical()
4556
10
4557
sage: Integer(0).radical()
4558
Traceback (most recent call last):
4559
...
4560
ArithmeticError: Radical of 0 not defined.
4561
"""
4562
if self.is_zero():
4563
raise ArithmeticError, "Radical of 0 not defined."
4564
return self.factor(*args, **kwds).radical_value()
4565
4566
def squarefree_part(self, long bound=-1):
4567
r"""
4568
Return the square free part of `x` (=self), i.e., the unique integer
4569
`z` that `x = z y^2`, with `y^2` a perfect square and `z` square-free.
4570
4571
Use ``self.radical()`` for the product of the primes that divide self.
4572
4573
If self is 0, just returns 0.
4574
4575
EXAMPLES::
4576
4577
sage: squarefree_part(100)
4578
1
4579
sage: squarefree_part(12)
4580
3
4581
sage: squarefree_part(17*37*37)
4582
17
4583
sage: squarefree_part(-17*32)
4584
-34
4585
sage: squarefree_part(1)
4586
1
4587
sage: squarefree_part(-1)
4588
-1
4589
sage: squarefree_part(-2)
4590
-2
4591
sage: squarefree_part(-4)
4592
-1
4593
4594
::
4595
4596
sage: a = 8 * 5^6 * 101^2
4597
sage: a.squarefree_part(bound=2).factor()
4598
2 * 5^6 * 101^2
4599
sage: a.squarefree_part(bound=5).factor()
4600
2 * 101^2
4601
sage: a.squarefree_part(bound=1000)
4602
2
4603
sage: a = 7^3 * next_prime(2^100)^2 * next_prime(2^200)
4604
sage: a / a.squarefree_part(bound=1000)
4605
49
4606
"""
4607
cdef Integer z
4608
cdef long even_part, p, p2
4609
cdef char switch_p
4610
if mpz_sgn(self.value) == 0:
4611
return self
4612
if 0 <= bound < 2:
4613
return self
4614
elif 2 <= bound <= 10000:
4615
z = PY_NEW(Integer)
4616
even_part = mpz_scan1(self.value, 0)
4617
mpz_fdiv_q_2exp(z.value, self.value, even_part ^ (even_part&1))
4618
sig_on()
4619
if bound >= 3:
4620
while mpz_divisible_ui_p(z.value, 9):
4621
mpz_divexact_ui(z.value, z.value, 9)
4622
if bound >= 5:
4623
while mpz_divisible_ui_p(z.value, 25):
4624
mpz_divexact_ui(z.value, z.value, 25)
4625
for p from 7 <= p <= bound by 2:
4626
switch_p = p % 30
4627
if switch_p in [1, 7, 11, 13, 17, 19, 23, 29]:
4628
p2 = p*p
4629
while mpz_divisible_ui_p(z.value, p2):
4630
mpz_divexact_ui(z.value, z.value, p2)
4631
sig_off()
4632
return z
4633
else:
4634
if bound == -1:
4635
F = self.factor()
4636
else:
4637
F = self._factor_trial_division(bound)
4638
n = one
4639
for pp, e in F:
4640
if e % 2 != 0:
4641
n = n * pp
4642
return n * F.unit()
4643
4644
def next_probable_prime(self):
4645
"""
4646
Returns the next probable prime after self, as determined by PARI.
4647
4648
EXAMPLES::
4649
4650
sage: (-37).next_probable_prime()
4651
2
4652
sage: (100).next_probable_prime()
4653
101
4654
sage: (2^512).next_probable_prime()
4655
13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084171
4656
sage: 0.next_probable_prime()
4657
2
4658
sage: 126.next_probable_prime()
4659
127
4660
sage: 144168.next_probable_prime()
4661
144169
4662
"""
4663
return Integer( self._pari_().nextprime(True) )
4664
4665
def next_prime(self, proof=None):
4666
r"""
4667
Returns the next prime after self.
4668
4669
INPUT:
4670
4671
- ``proof`` - bool or None (default: None, see
4672
proof.arithmetic or sage.structure.proof) Note that the global Sage
4673
default is proof=True
4674
4675
EXAMPLES::
4676
4677
sage: Integer(100).next_prime()
4678
101
4679
4680
Use Proof = False, which is way faster::
4681
4682
sage: b = (2^1024).next_prime(proof=False)
4683
4684
::
4685
4686
sage: Integer(0).next_prime()
4687
2
4688
sage: Integer(1001).next_prime()
4689
1009
4690
"""
4691
if mpz_cmp(self.value, PARI_PSEUDOPRIME_LIMIT) < 0:
4692
return Integer( self._pari_().nextprime(True) )
4693
if proof is None:
4694
from sage.structure.proof.proof import get_flag
4695
proof = get_flag(proof, "arithmetic")
4696
if not proof:
4697
return Integer( self._pari_().nextprime(True) )
4698
n = self
4699
if n % 2 == 0:
4700
n += 1
4701
else:
4702
n += 2
4703
while not n.is_prime(): # PARI's isprime() is provably correct
4704
n += 2
4705
return integer_ring.ZZ(n)
4706
4707
def additive_order(self):
4708
"""
4709
Return the additive order of self.
4710
4711
EXAMPLES::
4712
4713
sage: ZZ(0).additive_order()
4714
1
4715
sage: ZZ(1).additive_order()
4716
+Infinity
4717
"""
4718
import sage.rings.infinity
4719
if self.is_zero():
4720
return one
4721
else:
4722
return sage.rings.infinity.infinity
4723
4724
def multiplicative_order(self):
4725
r"""
4726
Return the multiplicative order of self.
4727
4728
EXAMPLES::
4729
4730
sage: ZZ(1).multiplicative_order()
4731
1
4732
sage: ZZ(-1).multiplicative_order()
4733
2
4734
sage: ZZ(0).multiplicative_order()
4735
+Infinity
4736
sage: ZZ(2).multiplicative_order()
4737
+Infinity
4738
"""
4739
import sage.rings.infinity
4740
if mpz_cmp_si(self.value, 1) == 0:
4741
return one
4742
elif mpz_cmp_si(self.value, -1) == 0:
4743
return Integer(2)
4744
else:
4745
return sage.rings.infinity.infinity
4746
4747
def is_squarefree(self):
4748
"""
4749
Returns True if this integer is not divisible by the square of any
4750
prime and False otherwise.
4751
4752
EXAMPLES::
4753
4754
sage: Integer(100).is_squarefree()
4755
False
4756
sage: Integer(102).is_squarefree()
4757
True
4758
sage: Integer(0).is_squarefree()
4759
False
4760
"""
4761
return self._pari_().issquarefree()
4762
4763
def _pari_(self):
4764
"""
4765
Returns the PARI version of this integer.
4766
4767
EXAMPLES::
4768
4769
sage: n = 9390823
4770
sage: m = n._pari_(); m
4771
9390823
4772
sage: type(m)
4773
<type 'sage.libs.pari.gen.gen'>
4774
4775
TESTS::
4776
4777
sage: n = 10^10000000
4778
sage: m = n._pari_() ## crash from trac 875
4779
sage: m % 1234567
4780
1041334
4781
"""
4782
return self._pari_c()
4783
4784
cdef _pari_c(self):
4785
cdef PariInstance P
4786
P = sage.libs.pari.gen.pari
4787
return P.new_gen_from_mpz_t(self.value)
4788
4789
def _interface_init_(self, I=None):
4790
"""
4791
Return canonical string to coerce this integer to any other math
4792
software, i.e., just the string representation of this integer in
4793
base 10.
4794
4795
EXAMPLES::
4796
4797
sage: n = 9390823
4798
sage: n._interface_init_()
4799
'9390823'
4800
"""
4801
return str(self)
4802
4803
property __array_interface__:
4804
def __get__(self):
4805
"""
4806
Used for NumPy conversion.
4807
4808
EXAMPLES::
4809
4810
sage: import numpy
4811
sage: numpy.array([1, 2, 3])
4812
array([1, 2, 3])
4813
sage: numpy.array([1, 2, 3]).dtype
4814
dtype('int32') # 32-bit
4815
dtype('int64') # 64-bit
4816
4817
sage: numpy.array(2**40).dtype
4818
dtype('int64')
4819
sage: numpy.array(2**400).dtype
4820
dtype('object')
4821
4822
sage: numpy.array([1,2,3,0.1]).dtype
4823
dtype('float64')
4824
"""
4825
if mpz_fits_slong_p(self.value):
4826
return numpy_long_interface
4827
elif sizeof(long) == 4 and mpz_sizeinbase(self.value, 2) <= 63:
4828
return numpy_int64_interface
4829
else:
4830
return numpy_object_interface
4831
4832
def _magma_init_(self, magma):
4833
"""
4834
Return string that evaluates in Magma to this element.
4835
4836
For small integers we just use base 10. For large integers we use
4837
base 16, but use Magma's StringToInteger command, which (for no
4838
good reason) is much faster than 0x[string literal]. We only use
4839
base 16 for integers with at least 10000 binary digits, since e.g.,
4840
for a large list of small integers the overhead of calling
4841
StringToInteger can be a killer.
4842
4843
EXAMPLES:
4844
sage: (117)._magma_init_(magma) # optional - magma
4845
'117'
4846
4847
Large integers use hex:
4848
sage: m = 3^(2^20) # optional - magma
4849
sage: s = m._magma_init_(magma) # optional - magma
4850
sage: 'StringToInteger' in s # optional - magma
4851
True
4852
sage: magma(m).sage() == m # optional - magma
4853
True
4854
"""
4855
if self.ndigits(2) > 10000:
4856
return 'StringToInteger("%s",16)'%self.str(16)
4857
else:
4858
return str(self)
4859
4860
def _sage_input_(self, sib, coerced):
4861
r"""
4862
Produce an expression which will reproduce this value when
4863
evaluated.
4864
4865
EXAMPLES::
4866
4867
sage: sage_input(1, verify=True)
4868
# Verified
4869
1
4870
sage: sage_input(1, preparse=False)
4871
ZZ(1)
4872
sage: sage_input(-12435, verify=True)
4873
# Verified
4874
-12435
4875
sage: sage_input(0, verify=True)
4876
# Verified
4877
0
4878
sage: sage_input(-3^70, verify=True)
4879
# Verified
4880
-2503155504993241601315571986085849
4881
sage: sage_input(-37, preparse=False)
4882
-ZZ(37)
4883
sage: sage_input(-37 * polygen(ZZ), preparse=False)
4884
R = ZZ['x']
4885
x = R.gen()
4886
-37*x
4887
sage: from sage.misc.sage_input import SageInputBuilder
4888
sage: (-314159)._sage_input_(SageInputBuilder(preparse=False), False)
4889
{unop:- {call: {atomic:ZZ}({atomic:314159})}}
4890
sage: (314159)._sage_input_(SageInputBuilder(preparse=False), True)
4891
{atomic:314159}
4892
"""
4893
if coerced or sib.preparse():
4894
return sib.int(self)
4895
else:
4896
if self < 0:
4897
return -sib.name('ZZ')(sib.int(-self))
4898
else:
4899
return sib.name('ZZ')(sib.int(self))
4900
4901
def sqrtrem(self):
4902
r"""
4903
Return (s, r) where s is the integer square root of self and
4904
r is the remainder such that `\text{self} = s^2 + r`.
4905
Raises ``ValueError`` if self is negative.
4906
4907
EXAMPLES::
4908
4909
sage: 25.sqrtrem()
4910
(5, 0)
4911
sage: 27.sqrtrem()
4912
(5, 2)
4913
sage: 0.sqrtrem()
4914
(0, 0)
4915
4916
::
4917
4918
sage: Integer(-102).sqrtrem()
4919
Traceback (most recent call last):
4920
...
4921
ValueError: square root of negative integer not defined.
4922
4923
"""
4924
if mpz_sgn(self.value) < 0:
4925
raise ValueError, "square root of negative integer not defined."
4926
cdef Integer s = PY_NEW(Integer)
4927
cdef Integer r = PY_NEW(Integer)
4928
mpz_sqrtrem(s.value, r.value, self.value)
4929
return s, r
4930
4931
def isqrt(self):
4932
r"""
4933
Returns the integer floor of the square root of self, or raises an
4934
``ValueError`` if self is negative.
4935
4936
EXAMPLE::
4937
4938
sage: a = Integer(5)
4939
sage: a.isqrt()
4940
2
4941
4942
::
4943
4944
sage: Integer(-102).isqrt()
4945
Traceback (most recent call last):
4946
...
4947
ValueError: square root of negative integer not defined.
4948
"""
4949
if mpz_sgn(self.value) < 0:
4950
raise ValueError, "square root of negative integer not defined."
4951
4952
cdef Integer x = PY_NEW(Integer)
4953
4954
sig_on()
4955
mpz_sqrt(x.value, self.value)
4956
sig_off()
4957
4958
return x
4959
4960
def sqrt_approx(self, prec=None, all=False):
4961
"""
4962
EXAMPLES::
4963
4964
sage: 5.sqrt_approx(prec=200)
4965
doctest:...: DeprecationWarning: This function is deprecated. Use sqrt with a given number of bits of precision instead.
4966
2.2360679774997896964091736687312762354406183596115257242709
4967
sage: 5.sqrt_approx()
4968
2.23606797749979
4969
sage: 4.sqrt_approx()
4970
2
4971
"""
4972
from sage.misc.misc import deprecation
4973
deprecation("This function is deprecated. Use sqrt with a given number of bits of precision instead.")
4974
try:
4975
return self.sqrt(extend=False,all=all)
4976
except ValueError:
4977
pass
4978
if prec is None:
4979
prec = max(53, 2*(mpz_sizeinbase(self.value, 2)+2))
4980
return self.sqrt(prec=prec, all=all)
4981
4982
def sqrt(self, prec=None, extend=True, all=False):
4983
"""
4984
The square root function.
4985
4986
INPUT:
4987
4988
- ``prec`` - integer (default: None): if None, returns
4989
an exact square root; otherwise returns a numerical square root if
4990
necessary, to the given bits of precision.
4991
4992
- ``extend`` - bool (default: True); if True, return a
4993
square root in an extension ring, if necessary. Otherwise, raise a
4994
ValueError if the square is not in the base ring. Ignored if prec
4995
is not None.
4996
4997
- ``all`` - bool (default: False); if True, return all
4998
square roots of self, instead of just one.
4999
5000
EXAMPLES::
5001
5002
sage: Integer(144).sqrt()
5003
12
5004
sage: sqrt(Integer(144))
5005
12
5006
sage: Integer(102).sqrt()
5007
sqrt(102)
5008
5009
::
5010
5011
sage: n = 2
5012
sage: n.sqrt(all=True)
5013
[sqrt(2), -sqrt(2)]
5014
sage: n.sqrt(prec=10)
5015
1.4
5016
sage: n.sqrt(prec=100)
5017
1.4142135623730950488016887242
5018
sage: n.sqrt(prec=100,all=True)
5019
[1.4142135623730950488016887242, -1.4142135623730950488016887242]
5020
sage: n.sqrt(extend=False)
5021
Traceback (most recent call last):
5022
...
5023
ValueError: square root of 2 not an integer
5024
sage: Integer(144).sqrt(all=True)
5025
[12, -12]
5026
sage: Integer(0).sqrt(all=True)
5027
[0]
5028
sage: type(Integer(5).sqrt())
5029
<type 'sage.symbolic.expression.Expression'>
5030
sage: type(Integer(5).sqrt(prec=53))
5031
<type 'sage.rings.real_mpfr.RealNumber'>
5032
sage: type(Integer(-5).sqrt(prec=53))
5033
<type 'sage.rings.complex_number.ComplexNumber'>
5034
"""
5035
if mpz_sgn(self.value) == 0:
5036
return [self] if all else self
5037
5038
if mpz_sgn(self.value) < 0:
5039
if not extend:
5040
raise ValueError, "square root of negative number not an integer"
5041
from sage.functions.other import _do_sqrt
5042
return _do_sqrt(self, prec=prec, all=all)
5043
5044
cdef int non_square
5045
cdef Integer z = PY_NEW(Integer)
5046
cdef mpz_t tmp
5047
sig_on()
5048
mpz_init(tmp)
5049
mpz_sqrtrem(z.value, tmp, self.value)
5050
non_square = mpz_sgn(tmp) != 0
5051
mpz_clear(tmp)
5052
sig_off()
5053
5054
if non_square:
5055
if not extend:
5056
raise ValueError, "square root of %s not an integer"%self
5057
from sage.functions.other import _do_sqrt
5058
return _do_sqrt(self, prec=prec, all=all)
5059
5060
if prec:
5061
from sage.functions.other import _do_sqrt
5062
return _do_sqrt(self, prec=prec, all=all)
5063
5064
if all:
5065
return [z, -z]
5066
return z
5067
5068
def _xgcd(self, Integer n, bint minimal=0):
5069
r"""
5070
Computes extended gcd of self and `n`.
5071
5072
INPUT:
5073
5074
- ``self`` - integer
5075
5076
- ``n`` - integer
5077
5078
- ``minimal`` - boolean (default false), whether to
5079
compute minimal cofactors (see below)
5080
5081
OUTPUT:
5082
5083
A triple (g, s, t), where `g` is the non-negative gcd of self
5084
and `n`, and `s` and `t` are cofactors satisfying the Bezout identity
5085
5086
.. math::
5087
5088
g = s \cdot \mbox{\rm self} + t \cdot n.
5089
5090
.. note::
5091
5092
If minimal is False, then there is no guarantee that the returned
5093
cofactors will be minimal in any sense; the only guarantee is that
5094
the Bezout identity will be satisfied (see examples below).
5095
5096
If minimal is True, the cofactors will satisfy the following
5097
conditions. If either self or `n` are zero, the trivial solution
5098
is returned. If both self and `n` are nonzero, the function returns
5099
the unique solution such that `0 \leq s < |n|/g` (which then must
5100
also satisfy `0 \leq |t| \leq |\mbox{\rm self}|/g`).
5101
5102
EXAMPLES::
5103
5104
sage: Integer(5)._xgcd(7)
5105
(1, 3, -2)
5106
sage: 5*3 + 7*-2
5107
1
5108
sage: g,s,t = Integer(58526524056)._xgcd(101294172798)
5109
sage: g
5110
22544886
5111
sage: 58526524056 * s + 101294172798 * t
5112
22544886
5113
5114
Without minimality guarantees, weird things can happen::
5115
5116
sage: Integer(3)._xgcd(21)
5117
(3, 1, 0)
5118
sage: Integer(3)._xgcd(24)
5119
(3, 1, 0)
5120
sage: Integer(3)._xgcd(48)
5121
(3, 1, 0)
5122
5123
Weirdness disappears with minimal option::
5124
5125
sage: Integer(3)._xgcd(21, minimal=True)
5126
(3, 1, 0)
5127
sage: Integer(3)._xgcd(24, minimal=True)
5128
(3, 1, 0)
5129
sage: Integer(3)._xgcd(48, minimal=True)
5130
(3, 1, 0)
5131
sage: Integer(21)._xgcd(3, minimal=True)
5132
(3, 0, 1)
5133
5134
Try minimal option with various edge cases::
5135
5136
sage: Integer(5)._xgcd(0, minimal=True)
5137
(5, 1, 0)
5138
sage: Integer(-5)._xgcd(0, minimal=True)
5139
(5, -1, 0)
5140
sage: Integer(0)._xgcd(5, minimal=True)
5141
(5, 0, 1)
5142
sage: Integer(0)._xgcd(-5, minimal=True)
5143
(5, 0, -1)
5144
sage: Integer(0)._xgcd(0, minimal=True)
5145
(0, 1, 0)
5146
5147
Exhaustive tests, checking minimality conditions::
5148
5149
sage: for a in srange(-20, 20):
5150
... for b in srange(-20, 20):
5151
... if a == 0 or b == 0: continue
5152
... g, s, t = a._xgcd(b)
5153
... assert g > 0
5154
... assert a % g == 0 and b % g == 0
5155
... assert a*s + b*t == g
5156
... g, s, t = a._xgcd(b, minimal=True)
5157
... assert g > 0
5158
... assert a % g == 0 and b % g == 0
5159
... assert a*s + b*t == g
5160
... assert s >= 0 and s < abs(b)/g
5161
... assert abs(t) <= abs(a)/g
5162
5163
AUTHORS:
5164
5165
- David Harvey (2007-12-26): added minimality option
5166
"""
5167
cdef Integer g = PY_NEW(Integer)
5168
cdef Integer s = PY_NEW(Integer)
5169
cdef Integer t = PY_NEW(Integer)
5170
5171
sig_on()
5172
mpz_gcdext(g.value, s.value, t.value, self.value, n.value)
5173
sig_off()
5174
5175
# Note: the GMP documentation for mpz_gcdext (or mpn_gcdext for that
5176
# matter) makes absolutely no claims about any minimality conditions
5177
# satisfied by the returned cofactors. They guarantee a non-negative
5178
# gcd, but that's it. So we have to do some work ourselves.
5179
5180
if not minimal:
5181
return g, s, t
5182
5183
# handle degenerate cases n == 0 and self == 0
5184
5185
if not mpz_sgn(n.value):
5186
mpz_set_ui(t.value, 0)
5187
mpz_abs(g.value, self.value)
5188
mpz_set_si(s.value, 1 if mpz_sgn(self.value) >= 0 else -1)
5189
return g, s, t
5190
5191
if not mpz_sgn(self.value):
5192
mpz_set_ui(s.value, 0)
5193
mpz_abs(g.value, n.value)
5194
mpz_set_si(t.value, 1 if mpz_sgn(n.value) >= 0 else -1)
5195
return g, s, t
5196
5197
# both n and self are nonzero, so we need to do a division and
5198
# make the appropriate adjustment
5199
5200
cdef mpz_t u1, u2
5201
mpz_init(u1)
5202
mpz_init(u2)
5203
mpz_divexact(u1, n.value, g.value)
5204
mpz_divexact(u2, self.value, g.value)
5205
if mpz_sgn(u1) > 0:
5206
mpz_fdiv_qr(u1, s.value, s.value, u1)
5207
else:
5208
mpz_cdiv_qr(u1, s.value, s.value, u1)
5209
mpz_addmul(t.value, u1, u2)
5210
mpz_clear(u2)
5211
mpz_clear(u1)
5212
5213
return g, s, t
5214
5215
cpdef _shift_helper(Integer self, y, int sign):
5216
"""
5217
Function used to compute left and right shifts of integers.
5218
Shifts self y bits to the left if sign is 1, and to the right
5219
if sign is -1.
5220
5221
WARNING: This function does no error checking. In particular,
5222
it assumes that sign is either 1 or -1,
5223
5224
EXAMPLES::
5225
5226
sage: n = 1234
5227
sage: factor(n)
5228
2 * 617
5229
sage: n._shift_helper(1, 1)
5230
2468
5231
sage: n._shift_helper(1, -1)
5232
617
5233
sage: n._shift_helper(100, 1)
5234
1564280840681635081446931755433984
5235
sage: n._shift_helper(100, -1)
5236
0
5237
sage: n._shift_helper(-100, 1)
5238
0
5239
sage: n._shift_helper(-100, -1)
5240
1564280840681635081446931755433984
5241
"""
5242
cdef long n
5243
5244
if PyInt_CheckExact(y):
5245
# For a Python int, we can just use the Python/C API.
5246
n = PyInt_AS_LONG(y)
5247
else:
5248
# If it's not already an Integer, try to convert it.
5249
if not PY_TYPE_CHECK(y, Integer):
5250
try:
5251
y = Integer(y)
5252
except TypeError:
5253
raise TypeError, "unsupported operands for %s: %s, %s"%(("<<" if sign == 1 else ">>"), self, y)
5254
except ValueError:
5255
return bin_op(self, y, operator.lshift if sign == 1 else operator.rshift)
5256
5257
# If y wasn't a Python int, it's now an Integer, so set n
5258
# accordingly.
5259
if mpz_fits_slong_p((<Integer>y).value):
5260
n = mpz_get_si((<Integer>y).value)
5261
elif sign * mpz_sgn((<Integer>y).value) < 0:
5262
# Doesn't fit in a long so shifting to the right by
5263
# this much will be 0.
5264
return PY_NEW(Integer)
5265
else:
5266
# Doesn't fit in a long so shifting to the left by
5267
# this much will raise appropriate overflow error
5268
n = y
5269
5270
# Decide which way we're shifting
5271
n *= sign
5272
5273
# Now finally call into MPIR to do the shifting.
5274
cdef Integer z = PY_NEW(Integer)
5275
if n < 0:
5276
mpz_fdiv_q_2exp(z.value, self.value, -n)
5277
else:
5278
mpz_mul_2exp(z.value, self.value, n)
5279
return z
5280
5281
def __lshift__(x, y):
5282
"""
5283
Shift x to the left by y bits.
5284
5285
EXAMPLES::
5286
5287
sage: 32 << 2
5288
128
5289
sage: 32 << int(2)
5290
128
5291
sage: int(32) << 2
5292
128
5293
sage: 1 << 2.5
5294
Traceback (most recent call last):
5295
...
5296
TypeError: unsupported operands for <<: 1, 2.5000...
5297
5298
sage: 32 << (4/2)
5299
128
5300
5301
A negative shift to the left is treated as a right shift::
5302
5303
sage: 128 << -2
5304
32
5305
sage: 128 << (-2^100)
5306
0
5307
"""
5308
# note that x need not be self -- int(3) << ZZ(2) will
5309
# dispatch this function
5310
if not PY_TYPE_CHECK(x, Integer):
5311
return x << int(y)
5312
return (<Integer>x)._shift_helper(y, 1)
5313
5314
def __rshift__(x, y):
5315
"""
5316
Shift x to the right by y bits.
5317
5318
EXAMPLES::
5319
5320
sage: 32 >> 2
5321
8
5322
sage: 32 >> int(2)
5323
8
5324
sage: int(32) >> 2
5325
8
5326
sage: 1 >> 2.5
5327
Traceback (most recent call last):
5328
...
5329
TypeError: unsupported operands for >>: 1, 2.5000...
5330
sage: 10^5 >> 10^100
5331
0
5332
5333
A negative shift to the right is treated as a left shift::
5334
5335
sage: 8 >> -2
5336
32
5337
"""
5338
# note that x need not be self -- int(3) >> ZZ(2) will
5339
# dispatch this function
5340
if not PY_TYPE_CHECK(x, Integer):
5341
return x >> int(y)
5342
return (<Integer>x)._shift_helper(y, -1)
5343
5344
cdef _and(Integer self, Integer other):
5345
cdef Integer x = PY_NEW(Integer)
5346
mpz_and(x.value, self.value, other.value)
5347
return x
5348
5349
def __and__(x, y):
5350
"""
5351
Return the bitwise and two integers.
5352
5353
EXAMPLES::
5354
5355
sage: n = Integer(6); m = Integer(2)
5356
sage: n & m
5357
2
5358
sage: n.__and__(m)
5359
2
5360
"""
5361
if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
5362
return (<Integer>x)._and(y)
5363
return bin_op(x, y, operator.and_)
5364
5365
cdef _or(Integer self, Integer other):
5366
cdef Integer x = PY_NEW(Integer)
5367
mpz_ior(x.value, self.value, other.value)
5368
return x
5369
5370
def __or__(x, y):
5371
"""
5372
Return the bitwise or of the integers x and y.
5373
5374
EXAMPLES::
5375
5376
sage: n = 8; m = 4
5377
sage: n.__or__(m)
5378
12
5379
"""
5380
if PY_TYPE_CHECK(x, Integer) and PY_TYPE_CHECK(y, Integer):
5381
return (<Integer>x)._or(y)
5382
return bin_op(x, y, operator.or_)
5383
5384
def __invert__(self):
5385
"""
5386
Return the multiplicative inverse of self, as a rational number.
5387
5388
EXAMPLE::
5389
5390
sage: n = 10
5391
sage: 1/n
5392
1/10
5393
sage: n.__invert__()
5394
1/10
5395
"""
5396
return one / self
5397
5398
def inverse_of_unit(self):
5399
"""
5400
Return inverse of self if self is a unit in the integers, i.e.,
5401
self is -1 or 1. Otherwise, raise a ZeroDivisionError.
5402
5403
EXAMPLES::
5404
5405
sage: (1).inverse_of_unit()
5406
1
5407
sage: (-1).inverse_of_unit()
5408
-1
5409
sage: 5.inverse_of_unit()
5410
Traceback (most recent call last):
5411
...
5412
ZeroDivisionError: Inverse does not exist.
5413
sage: 0.inverse_of_unit()
5414
Traceback (most recent call last):
5415
...
5416
ZeroDivisionError: Inverse does not exist.
5417
"""
5418
if mpz_cmpabs_ui(self.value, 1) == 0:
5419
return self
5420
else:
5421
raise ZeroDivisionError, "Inverse does not exist."
5422
5423
def inverse_mod(self, n):
5424
"""
5425
Returns the inverse of self modulo `n`, if this inverse exists.
5426
Otherwise, raises a ``ZeroDivisionError`` exception.
5427
5428
INPUT:
5429
5430
- ``self`` - Integer
5431
5432
- ``n`` - Integer, or ideal of integer ring
5433
5434
OUTPUT:
5435
5436
- ``x`` - Integer such that x\*self = 1 (mod m), or
5437
raises ZeroDivisionError.
5438
5439
IMPLEMENTATION:
5440
5441
Call the mpz_invert GMP library function.
5442
5443
EXAMPLES::
5444
5445
sage: a = Integer(189)
5446
sage: a.inverse_mod(10000)
5447
4709
5448
sage: a.inverse_mod(-10000)
5449
4709
5450
sage: a.inverse_mod(1890)
5451
Traceback (most recent call last):
5452
...
5453
ZeroDivisionError: Inverse does not exist.
5454
sage: a = Integer(19)**100000
5455
sage: b = a*a
5456
sage: c = a.inverse_mod(b)
5457
Traceback (most recent call last):
5458
...
5459
ZeroDivisionError: Inverse does not exist.
5460
5461
We check that #10625 is fixed::
5462
5463
sage: ZZ(2).inverse_mod(ZZ.ideal(3))
5464
2
5465
5466
We check that #9955 is fixed::
5467
5468
sage: Rational(3)%Rational(-1)
5469
0
5470
"""
5471
cdef int r
5472
if isinstance(n, sage.rings.ideal.Ideal_pid) and n.ring() == the_integer_ring:
5473
n = n.gen()
5474
cdef Integer m = as_Integer(n)
5475
cdef Integer ans = <Integer>PY_NEW(Integer)
5476
if mpz_cmpabs_ui(m.value, 1) == 0:
5477
return zero
5478
sig_on()
5479
r = mpz_invert(ans.value, self.value, m.value)
5480
sig_off()
5481
if r == 0:
5482
raise ZeroDivisionError, "Inverse does not exist."
5483
return ans
5484
5485
def gcd(self, n):
5486
"""
5487
Return the greatest common divisor of self and `n`.
5488
5489
EXAMPLE::
5490
5491
sage: gcd(-1,1)
5492
1
5493
sage: gcd(0,1)
5494
1
5495
sage: gcd(0,0)
5496
0
5497
sage: gcd(2,2^6)
5498
2
5499
sage: gcd(21,2^6)
5500
1
5501
"""
5502
if not isinstance(n, Integer) and not isinstance(n, int):
5503
left, right = canonical_coercion(self, n)
5504
return left.gcd(right)
5505
cdef Integer m = as_Integer(n)
5506
cdef Integer g = PY_NEW(Integer)
5507
sig_on()
5508
mpz_gcd(g.value, self.value, m.value)
5509
sig_off()
5510
return g
5511
5512
def crt(self, y, m, n):
5513
"""
5514
Return the unique integer between `0` and `mn` that is congruent to
5515
the integer modulo `m` and to `y` modulo `n`. We assume that `m` and
5516
`n` are coprime.
5517
5518
EXAMPLES::
5519
5520
sage: n = 17
5521
sage: m = n.crt(5, 23, 11); m
5522
247
5523
sage: m%23
5524
17
5525
sage: m%11
5526
5
5527
"""
5528
cdef object g, s, t
5529
cdef Integer _y, _m, _n
5530
_y = Integer(y); _m = Integer(m); _n = Integer(n)
5531
g, s, t = _m.xgcd(_n)
5532
if not g.is_one():
5533
raise ArithmeticError, "CRT requires that gcd of moduli is 1."
5534
# Now s*m + t*n = 1, so the answer is x + (y-x)*s*m, where x=self.
5535
return (self + (_y - self) * s * _m) % (_m * _n)
5536
5537
def test_bit(self, long index):
5538
r"""
5539
Return the bit at ``index``.
5540
5541
If the index is negative, returns 0.
5542
5543
Although internally a sign-magnitude representation is used
5544
for integers, this method pretends to use a two's complement
5545
representation. This is illustrated with a negative integer
5546
below.
5547
5548
EXAMPLES::
5549
5550
sage: w = 6
5551
sage: w.str(2)
5552
'110'
5553
sage: w.test_bit(2)
5554
1
5555
sage: w.test_bit(-1)
5556
0
5557
sage: x = -20
5558
sage: x.str(2)
5559
'-10100'
5560
sage: x.test_bit(4)
5561
0
5562
sage: x.test_bit(5)
5563
1
5564
sage: x.test_bit(6)
5565
1
5566
"""
5567
if index < 0:
5568
return 0
5569
else:
5570
return mpz_tstbit(self.value, index)
5571
5572
def popcount(self):
5573
"""
5574
Return the number of 1 bits in the binary representation.
5575
If self<0, we return Infinity.
5576
5577
EXAMPLES::
5578
5579
sage: n = 123
5580
sage: n.str(2)
5581
'1111011'
5582
sage: n.popcount()
5583
6
5584
5585
sage: n = -17
5586
sage: n.popcount()
5587
+Infinity
5588
"""
5589
if self<0:
5590
return sage.rings.infinity.Infinity
5591
return int(mpz_popcount(self.value))
5592
5593
5594
def conjugate(self):
5595
"""
5596
Return the complex conjugate of this integer, which is the
5597
integer itself.
5598
5599
EXAMPLES:
5600
sage: n = 205
5601
sage: n.conjugate()
5602
205
5603
"""
5604
return self
5605
5606
def binomial(self, m, algorithm='mpir'):
5607
"""
5608
Return the binomial coefficient "self choose m".
5609
5610
INPUT:
5611
5612
- ``m`` -- an integer
5613
5614
- ``algorithm`` -- ``'mpir'`` (default) or ``'pari'``; ``'mpir'`` is
5615
faster for small ``m``, and ``'pari'`` tends to be faster for
5616
large ``m``
5617
5618
OUTPUT:
5619
5620
- integer
5621
5622
EXAMPLES::
5623
5624
sage: 10.binomial(2)
5625
45
5626
sage: 10.binomial(2, algorithm='pari')
5627
45
5628
sage: 10.binomial(-2)
5629
0
5630
"""
5631
cdef Integer x
5632
if m < 0:
5633
return the_integer_ring.zero()
5634
if algorithm == 'mpir':
5635
x = PY_NEW(Integer)
5636
sig_on()
5637
mpz_bin_ui(x.value, self.value, m)
5638
sig_off()
5639
return x
5640
elif algorithm == 'pari':
5641
return the_integer_ring(self._pari_().binomial(m))
5642
else:
5643
raise ValueError("algorithm must be one of: 'pari', 'mpir'")
5644
5645
5646
ONE = Integer(1)
5647
Py_INCREF(ONE)
5648
5649
cpdef LCM_list(v):
5650
"""
5651
Return the LCM of a list v of integers. Elements of v are converted
5652
to Sage integers if they aren't already.
5653
5654
This function is used, e.g., by rings/arith.py
5655
5656
INPUT:
5657
5658
- ``v`` - list or tuple
5659
5660
OUTPUT: integer
5661
5662
EXAMPLES::
5663
5664
sage: from sage.rings.integer import LCM_list
5665
sage: w = LCM_list([3,9,30]); w
5666
90
5667
sage: type(w)
5668
<type 'sage.rings.integer.Integer'>
5669
5670
The inputs are converted to Sage integers.
5671
5672
::
5673
5674
sage: w = LCM_list([int(3), int(9), '30']); w
5675
90
5676
sage: type(w)
5677
<type 'sage.rings.integer.Integer'>
5678
"""
5679
cdef int i, n = len(v)
5680
cdef Integer z = <Integer>PY_NEW(Integer)
5681
5682
for i from 0 <= i < n:
5683
if not isinstance(v[i], Integer):
5684
if not isinstance(v, list):
5685
v = list(v)
5686
v[i] = Integer(v[i])
5687
5688
if n == 0:
5689
return one
5690
elif n == 1:
5691
return v[0].abs()
5692
5693
sig_on()
5694
mpz_lcm(z.value, (<Integer>v[0]).value, (<Integer>v[1]).value)
5695
for i from 2 <= i < n:
5696
mpz_lcm(z.value, z.value, (<Integer>v[i]).value)
5697
sig_off()
5698
5699
return z
5700
5701
def GCD_list(v):
5702
r"""
5703
Return the GCD of a list v of integers. Elements of v are converted
5704
to Sage integers if they aren't already.
5705
5706
This function is used, e.g., by rings/arith.py
5707
5708
INPUT:
5709
5710
- ``v`` - list or tuple
5711
5712
OUTPUT: integer
5713
5714
EXAMPLES::
5715
5716
sage: from sage.rings.integer import GCD_list
5717
sage: w = GCD_list([3,9,30]); w
5718
3
5719
sage: type(w)
5720
<type 'sage.rings.integer.Integer'>
5721
5722
Check that the bug reported in trac #3118 has been fixed::
5723
5724
sage: sage.rings.integer.GCD_list([2,2,3])
5725
1
5726
5727
The inputs are converted to Sage integers.
5728
5729
::
5730
5731
sage: w = GCD_list([int(3), int(9), '30']); w
5732
3
5733
sage: type(w)
5734
<type 'sage.rings.integer.Integer'>
5735
"""
5736
cdef int i, n = len(v)
5737
cdef Integer z = <Integer>PY_NEW(Integer)
5738
5739
for i from 0 <= i < n:
5740
if not isinstance(v[i], Integer):
5741
if not isinstance(v, list):
5742
v = list(v)
5743
v[i] = Integer(v[i])
5744
5745
if n == 0:
5746
return one
5747
elif n == 1:
5748
return v[0].abs()
5749
5750
sig_on()
5751
mpz_gcd(z.value, (<Integer>v[0]).value, (<Integer>v[1]).value)
5752
for i from 2 <= i < n:
5753
if mpz_cmp_ui(z.value, 1) == 0:
5754
break
5755
mpz_gcd(z.value, z.value, (<Integer>v[i]).value)
5756
sig_off()
5757
5758
return z
5759
5760
def make_integer(s):
5761
"""
5762
Create a Sage integer from the base-32 Python *string* s. This is
5763
used in unpickling integers.
5764
5765
EXAMPLES::
5766
5767
sage: from sage.rings.integer import make_integer
5768
sage: make_integer('-29')
5769
-73
5770
sage: make_integer(29)
5771
Traceback (most recent call last):
5772
...
5773
TypeError: expected string or Unicode object, sage.rings.integer.Integer found
5774
"""
5775
cdef Integer r = PY_NEW(Integer)
5776
r._reduce_set(s)
5777
return r
5778
5779
cdef class int_to_Z(Morphism):
5780
def __init__(self):
5781
import integer_ring
5782
import sage.categories.homset
5783
from sage.structure.parent import Set_PythonType
5784
Morphism.__init__(self, sage.categories.homset.Hom(Set_PythonType(int), integer_ring.ZZ))
5785
cpdef Element _call_(self, a):
5786
cdef Integer r
5787
r = <Integer>PY_NEW(Integer)
5788
mpz_set_si(r.value, PyInt_AS_LONG(a))
5789
return r
5790
def _repr_type(self):
5791
return "Native"
5792
5793
cdef class long_to_Z(Morphism):
5794
"""
5795
EXAMPLES::
5796
5797
sage: f = ZZ.coerce_map_from(long); f
5798
Native morphism:
5799
From: Set of Python objects of type 'long'
5800
To: Integer Ring
5801
sage: f(1rL)
5802
1
5803
sage: f(-10000000000000000000001r)
5804
-10000000000000000000001
5805
"""
5806
def __init__(self):
5807
import integer_ring
5808
import sage.categories.homset
5809
from sage.structure.parent import Set_PythonType
5810
Morphism.__init__(self, sage.categories.homset.Hom(Set_PythonType(long), integer_ring.ZZ))
5811
cpdef Element _call_(self, a):
5812
cdef Integer r
5813
r = <Integer>PY_NEW(Integer)
5814
mpz_set_pylong(r.value, a)
5815
return r
5816
def _repr_type(self):
5817
return "Native"
5818
5819
############### INTEGER CREATION CODE #####################
5820
5821
# We need a couple of internal GMP datatypes.
5822
5823
# This may be potentially very dangerous as it reaches
5824
# deeply into the internal structure of GMP which may not
5825
# be consistent across future versions of GMP.
5826
# See extensive note in the fast_tp_new() function below.
5827
5828
include "../ext/python_rich_object.pxi"
5829
cdef extern from "gmp.h":
5830
ctypedef long mp_limb_t
5831
ctypedef mp_limb_t* mp_ptr #"mp_ptr"
5832
5833
# We allocate _mp_d directly (mpz_t is typedef of this in GMP)
5834
ctypedef struct __mpz_struct "__mpz_struct":
5835
int _mp_alloc
5836
int _mp_size
5837
mp_ptr _mp_d
5838
5839
# sets the three free, alloc, and realloc function pointers to the
5840
# memory management functions set in GMP. Accepts NULL pointer.
5841
# Potentially dangerous if changed by calling
5842
# mp_set_memory_functions again after we initialized this module.
5843
void mp_get_memory_functions (void *(**alloc) (size_t), void *(**realloc)(void *, size_t, size_t), void (**free) (void *, size_t))
5844
5845
# GMP's configuration of how many Bits are stuffed into a limb
5846
cdef int GMP_LIMB_BITS
5847
5848
# This variable holds the size of any Integer object in bytes.
5849
cdef int sizeof_Integer
5850
5851
# We use a global Integer element to steal all the references
5852
# from. DO NOT INITIALIZE IT AGAIN and DO NOT REFERENCE IT!
5853
cdef Integer global_dummy_Integer
5854
global_dummy_Integer = Integer()
5855
5856
# Accessing the .value attribute of an Integer object causes Cython to
5857
# refcount it. This is problematic, because that causes overhead and
5858
# more importantly an infinite loop in the destructor. If you refcount
5859
# in the destructor and the refcount reaches zero (which is true
5860
# every time) the destructor is called.
5861
#
5862
# To avoid this we calculate the byte offset of the value member and
5863
# remember it in this variable.
5864
#
5865
# Eventually this may be rendered obsolete by a change in Cython allowing
5866
# non-reference counted extension types.
5867
cdef long mpz_t_offset
5868
mpz_t_offset_python = None
5869
5870
5871
# stores the GMP alloc function
5872
cdef void * (* mpz_alloc)(size_t)
5873
5874
# stores the GMP free function
5875
cdef void (* mpz_free)(void *, size_t)
5876
5877
# A global pool for performance when integers are rapidly created and destroyed.
5878
# It operates on the following principles:
5879
#
5880
# - The pool starts out empty.
5881
# - When an new integer is needed, one from the pool is returned
5882
# if available, otherwise a new Integer object is created
5883
# - When an integer is collected, it will add it to the pool
5884
# if there is room, otherwise it will be deallocated.
5885
5886
5887
cdef int integer_pool_size = 100
5888
5889
cdef PyObject** integer_pool
5890
cdef int integer_pool_count = 0
5891
5892
# used for profiling the pool
5893
cdef int total_alloc = 0
5894
cdef int use_pool = 0
5895
5896
# The signature of tp_new is
5897
# PyObject* tp_new(RichPyTypeObject *t, PyObject *a, PyObject *k).
5898
# However we only use t in this implementation.
5899
#
5900
# t in this case is the Integer TypeObject.
5901
5902
cdef PyObject* fast_tp_new(RichPyTypeObject *t, PyObject *a, PyObject *k):
5903
5904
global integer_pool, integer_pool_count, total_alloc, use_pool
5905
5906
cdef RichPyObject* new
5907
5908
# for profiling pool usage
5909
# total_alloc += 1
5910
5911
# If there is a ready integer in the pool, we will
5912
# decrement the counter and return that.
5913
5914
if integer_pool_count > 0:
5915
5916
# for profiling pool usage
5917
# use_pool += 1
5918
5919
integer_pool_count -= 1
5920
new = <RichPyObject *> integer_pool[integer_pool_count]
5921
5922
# Otherwise, we have to create one.
5923
5924
else:
5925
5926
# allocate enough room for the Integer, sizeof_Integer is
5927
# sizeof(Integer). The use of PyObject_MALLOC directly
5928
# assumes that Integers are not garbage collected, i.e.
5929
# they do not possess references to other Python
5930
# objects (as indicated by the Py_TPFLAGS_HAVE_GC flag).
5931
# See below for a more detailed description.
5932
new = PyObject_MALLOC( sizeof_Integer )
5933
5934
# Now set every member as set in z, the global dummy Integer
5935
# created before this tp_new started to operate.
5936
5937
memcpy(new, (<void*>global_dummy_Integer), sizeof_Integer )
5938
5939
# This line is only needed if Python is compiled in debugging
5940
# mode './configure --with-pydebug'. If that is the case a Python
5941
# object has a bunch of debugging fields which are initialized
5942
# with this macro. For speed reasons, we don't call it if Python
5943
# is not compiled in debug mode. So uncomment the following line
5944
# if you are debugging Python.
5945
5946
#PyObject_INIT(new, (<RichPyObject*>global_dummy_Integer).ob_type)
5947
5948
# We take the address 'new' and move mpz_t_offset bytes (chars)
5949
# to the address of 'value'. We treat that address as a pointer
5950
# to a mpz_t struct and allocate memory for the _mp_d element of
5951
# that struct. We allocate one limb.
5952
#
5953
# What is done here is potentially very dangerous as it reaches
5954
# deeply into the internal structure of GMP. Consequently things
5955
# may break if a new release of GMP changes some internals. To
5956
# emphasize this, this is what the GMP manual has to say about
5957
# the documentation for the struct we are using:
5958
#
5959
# "This chapter is provided only for informational purposes and the
5960
# various internals described here may change in future GMP releases.
5961
# Applications expecting to be compatible with future releases should use
5962
# only the documented interfaces described in previous chapters."
5963
#
5964
# If this line is used Sage is not such an application.
5965
#
5966
# The clean version of the following line is:
5967
#
5968
# mpz_init( <mpz_t>(<char *>new + mpz_t_offset) )
5969
#
5970
# We save time both by avoiding an extra function call and
5971
# because the rest of the mpz struct was already initialized
5972
# fully using the memcpy above.
5973
5974
(<__mpz_struct *>( <char *>new + mpz_t_offset) )._mp_d = <mp_ptr>mpz_alloc(GMP_LIMB_BITS >> 3)
5975
5976
# The global_dummy_Integer may have a reference count larger than
5977
# one, but it is expected that newly created objects have a
5978
# reference count of one. This is potentially unneeded if
5979
# everybody plays nice, because the gobal_dummy_Integer has only
5980
# one reference in that case.
5981
5982
# Objects from the pool have reference count zero, so this
5983
# needs to be set in this case.
5984
5985
new.ob_refcnt = 1
5986
5987
return new
5988
5989
cdef void fast_tp_dealloc(PyObject* o):
5990
5991
# If there is room in the pool for a used integer object,
5992
# then put it in rather than deallocating it.
5993
5994
global integer_pool, integer_pool_count
5995
5996
if integer_pool_count < integer_pool_size:
5997
5998
# Here we free any extra memory used by the mpz_t by
5999
# setting it to a single limb.
6000
if (<__mpz_struct *>( <char *>o + mpz_t_offset))._mp_alloc > 10:
6001
_mpz_realloc(<mpz_t *>( <char *>o + mpz_t_offset), 10)
6002
6003
# It's cheap to zero out an integer, so do it here.
6004
(<__mpz_struct *>( <char *>o + mpz_t_offset))._mp_size = 0
6005
6006
# And add it to the pool.
6007
integer_pool[integer_pool_count] = o
6008
integer_pool_count += 1
6009
return
6010
6011
# Again, we move to the mpz_t and clear it. See above, why this is evil.
6012
# The clean version of this line would be:
6013
# mpz_clear(<mpz_t>(<char *>o + mpz_t_offset))
6014
6015
mpz_free((<__mpz_struct *>( <char *>o + mpz_t_offset) )._mp_d, 0)
6016
6017
# Free the object. This assumes that Py_TPFLAGS_HAVE_GC is not
6018
# set. If it was set another free function would need to be
6019
# called.
6020
6021
PyObject_FREE(o)
6022
6023
6024
hook_fast_tp_functions()
6025
from sage.misc.allocator cimport hook_tp_functions
6026
cdef hook_fast_tp_functions():
6027
"""
6028
Initialize the fast integer creation functions.
6029
"""
6030
global global_dummy_Integer, mpz_t_offset, sizeof_Integer, integer_pool
6031
6032
integer_pool = <PyObject**>sage_malloc(integer_pool_size * sizeof(PyObject*))
6033
6034
cdef RichPyObject* o
6035
o = <RichPyObject*>global_dummy_Integer
6036
6037
# calculate the offset of the GMP mpz_t to avoid casting to/from
6038
# an Integer which includes reference counting. Reference counting
6039
# is bad in constructors and destructors as it potentially calls
6040
# the destructor.
6041
# Eventually this may be rendered obsolete by a change in Cython allowing
6042
# non-reference counted extension types.
6043
mpz_t_offset = <char *>(&global_dummy_Integer.value) - <char *>o
6044
global mpz_t_offset_python
6045
mpz_t_offset_python = mpz_t_offset
6046
6047
# store how much memory needs to be allocated for an Integer.
6048
sizeof_Integer = o.ob_type.tp_basicsize
6049
6050
# get the functions to do memory management for the GMP elements
6051
# WARNING: if the memory management functions are changed after
6052
# this initialisation, we are/you are doomed.
6053
6054
mp_get_memory_functions(&mpz_alloc, NULL, &mpz_free)
6055
6056
# Finally replace the functions called when an Integer needs
6057
# to be constructed/destructed.
6058
hook_tp_functions(global_dummy_Integer, NULL, <newfunc>&fast_tp_new, NULL, &fast_tp_dealloc, False)
6059
6060
cdef integer(x):
6061
if PY_TYPE_CHECK(x, Integer):
6062
return x
6063
return Integer(x)
6064
6065
6066
def free_integer_pool():
6067
cdef int i
6068
cdef PyObject *o
6069
6070
global integer_pool_count, integer_pool_size
6071
6072
for i from 0 <= i < integer_pool_count:
6073
o = integer_pool[i]
6074
mpz_clear(<__mpz_struct *>(<char*>o + mpz_t_offset))
6075
# Free the object. This assumes that Py_TPFLAGS_HAVE_GC is not
6076
# set. If it was set another free function would need to be
6077
# called.
6078
PyObject_FREE(o)
6079
6080
integer_pool_size = 0
6081
integer_pool_count = 0
6082
sage_free(integer_pool)
6083
6084
6085