Coverage for local/lib/python2.7/site-packages/sage/arith/multi_modular.pyx : 73%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" Utility classes for multi-modular algorithms """
#***************************************************************************** # Copyright (C) 2006 William Stein # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #***************************************************************************** from __future__ import absolute_import
from cysignals.memory cimport check_allocarray, check_reallocarray, sig_free
from sage.libs.gmp.mpz cimport * from sage.rings.integer cimport Integer, smallInteger from sage.arith.all import random_prime from types import GeneratorType from sage.ext.stdsage cimport PY_NEW from cpython.object cimport PyObject_RichCompare
# should I have mod_int versions of these functions? # c_inverse_mod_longlong modular inverse used exactly once in _refresh_precomputations from sage.rings.fast_arith cimport arith_llong cdef arith_llong ai ai = arith_llong()
# This is the maximum modulus for the code in this module, i.e. the # largest prime that can be used as modulus is # previous_prime(MAX_MODULUS) MAX_MODULUS = MOD_INT_MAX
cdef class MultiModularBasis_base(object): r""" This class stores a list of machine-sized prime numbers, and can do reduction and Chinese Remainder Theorem lifting modulo these primes.
Lifting implemented via Garner's algorithm, which has the advantage that all reductions are word-sized. For each `i`, precompute `\prod_j=1^{i-1} m_j` and `\prod_j=1^{i-1} m_j^{-1} (mod m_i)`.
This class can be initialized in two ways, either with a list of prime moduli or an upper bound for the product of the prime moduli. The prime moduli are generated automatically in the second case.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([3, 5, 7]); mm MultiModularBasis with moduli [3, 5, 7]
sage: height = 52348798724 sage: mm = MultiModularBasis_base(height); mm MultiModularBasis with moduli [4561, 17351, 28499] sage: mm = MultiModularBasis_base(height); mm MultiModularBasis with moduli [32573, 4339, 30859] sage: mm = MultiModularBasis_base(height); mm MultiModularBasis with moduli [16451, 14323, 28631]
sage: mm.prod()//height 128
TESTS::
sage: mm = MultiModularBasis_base((3,5,7)); mm MultiModularBasis with moduli [3, 5, 7] sage: mm = MultiModularBasis_base(primes(10,20)); mm MultiModularBasis with moduli [11, 13, 17, 19]
There is no overflow if the modulus is below ``MAX_MODULUS``::
sage: from sage.arith.multi_modular import MAX_MODULUS sage: p0 = previous_prime(MAX_MODULUS) sage: p1 = previous_prime(p0) sage: MultiModularBasis_base([p0, p1]).crt([p0-1, p1-1]) -1
If we add another bit to the prime length then there is an overflow, as expected::
sage: p0 = previous_prime(2*MAX_MODULUS) sage: p1 = previous_prime(p0) sage: MultiModularBasis_base([p0, p1]).crt([p0-1, p1-1]) Traceback (most recent call last): ... OverflowError: given modulus 6074000981 is larger than 3037000498 """ def __cinit__(self): r""" Allocate the space for the moduli and precomputation lists and initialize the first element of that list.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([1099511627791]) Traceback (most recent call last): ... OverflowError: given modulus 1099511627791 is larger than 3037000498 """
cdef _realloc_to_new_count(self, new_count):
def __dealloc__(self): """ TESTS::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base(1099511627791); mm MultiModularBasis with moduli [4561, 17351, 28499] sage: del mm """
def __init__(self, val, unsigned long l_bound=2**10, unsigned long u_bound=2**15): r""" Initialize a multi-modular basis and perform precomputations.
INPUT:
- ``val`` -- as integer determines how many primes are computed (their product will be at least 2*val) as list, tuple or generator a list of prime moduli to start with - ``l_bound`` -- an integer: lower bound for the random primes (default: 2^10) - ``u_bound`` -- an integer: upper bound for the random primes (default: 2^15)
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([1009, 10007]); mm MultiModularBasis with moduli [1009, 10007] sage: mm.prod() 10097063
sage: height = 10097063 sage: mm = MultiModularBasis_base(height); mm MultiModularBasis with moduli [...]
sage: mm.prod()//height > 2 True
sage: mm = MultiModularBasis_base([1000000000000000000000000000057]) Traceback (most recent call last): ... OverflowError: given modulus 1000000000000000000000000000057 is larger than 3037000498
sage: mm = MultiModularBasis_base(0); mm MultiModularBasis with moduli [28499]
sage: mm = MultiModularBasis_base([6, 10]) Traceback (most recent call last): ... ArithmeticError: The inverse of 6 modulo 10 is not defined. """ raise ValueError(f"minimum value for lower bound is 2, given: {l_bound}") raise ValueError(f"upper bound cannot be greater than {MAX_MODULUS}, given: {u_bound}")
else:
cdef mod_int _new_random_prime(self, set known_primes) except 1: """ Choose a new random prime for inclusion in the list of moduli, or raise a ``RuntimeError`` if there are no more primes.
INPUT:
- ``known_primes`` -- Python set of already known primes in the allowed interval; we will not return a prime in known_primes. """ cdef Py_ssize_t i cdef mod_int p
def extend_with_primes(self, plist, partial_products = None, check=True): """ Extend the stored list of moduli with the given primes in ``plist``.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([1009, 10007]); mm MultiModularBasis with moduli [1009, 10007] sage: mm.extend_with_primes([10037, 10039]) 4 sage: mm MultiModularBasis with moduli [1009, 10007, 10037, 10039] """ raise TypeError("plist should be a list, tuple or a generator, got: %s"%(str(type(plist))))
cdef Py_ssize_t i
else:
def __richcmp__(self, other, int op): """ EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007]) sage: nn = MultiModularBasis_base([10007]) sage: mm == nn True """ return NotImplemented
def __setstate__(self, state): """ Initialize a new :class:`MultiModularBasis_base` object from a state stored in a pickle.
TESTS::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007, 10009]) sage: mm == loads(dumps(mm)) True
sage: mm = MultiModularBasis_base([]) sage: mm.__setstate__(([10007, 10009], 2^10, 2^15))
sage: mm MultiModularBasis with moduli [10007, 10009] """
def __getstate__(self): """ Return a tuple describing the state of this object for pickling.
TESTS::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007, 10009]) sage: mm.__getstate__() ([10007, 10009], 1024L, 32768L) """
def _extend_moduli_to_height(self, height): """ EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base(0); mm MultiModularBasis with moduli [4561]
sage: mm._extend_moduli_to_height(10000) sage: mm MultiModularBasis with moduli [4561, 17351]
sage: mm = MultiModularBasis_base([46307]); mm MultiModularBasis with moduli [46307]
sage: mm._extend_moduli_to_height(10^30); mm MultiModularBasis with moduli [46307, 28499, 32573, 4339, 30859, 16451, 14323, 28631]
TESTS:
Verify that :trac:`11358` is fixed::
sage: set_random_seed(0); m = sage.arith.multi_modular.MultiModularBasis_base(0) sage: m._extend_moduli_to_height(prod(prime_range(50))) sage: m = sage.arith.multi_modular.MultiModularBasis_base([],2,100) sage: m._extend_moduli_to_height(prod(prime_range(90))) sage: m._extend_moduli_to_height(prod(prime_range(150))) Traceback (most recent call last): ... RuntimeError: there are not enough primes in the interval [2, 100] to complete this multimodular computation
Another check (which fails horribly before :trac:`11358` is fixed)::
sage: set_random_seed(0); m = sage.arith.multi_modular.MultiModularBasis_base(0); m._extend_moduli_to_height(10**10000) sage: len(set(m)) == len(m) True sage: len(m) 2438 """
cdef int _extend_moduli_to_height_c(self, mpz_t height0) except -1: r""" Expand the list of primes and perform precomputations.
INPUT:
- ``height`` -- determines how many primes are computed (their product must be at least 2*height) """ # real height we use is twice the given, set height to 2*height0 cdef mpz_t height # check if we already have enough prime moduli mpz_clear(height) return self.n
# find new prime moduli cdef int i cdef Integer M # keeps current height cdef mod_int p # keeps current prime moduli
else:
def _extend_moduli_to_count(self, int count): r""" Expand the list of primes and perform precomputations.
INPUT:
- ``count`` -- the minimum number of moduli in the resulting list
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307]); mm MultiModularBasis with moduli [46307] sage: mm._extend_moduli_to_count(3) 3 sage: mm MultiModularBasis with moduli [46307, 4561, 17351] """ return self.n
cdef int i cdef mod_int p
def _extend_moduli(self, count): """ Expand the list of prime moduli with `count` new random primes.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307]); mm MultiModularBasis with moduli [46307] sage: mm._extend_moduli(2); mm MultiModularBasis with moduli [46307, 4561, 17351] """
cdef void _refresh_products(self, int start): r""" Compute and store `\prod_j=1^{i-1} m_j` for i > start. """ cdef mpz_t z
cdef void _refresh_prod(self): # record the product and half product for balancing the lifts.
cdef void _refresh_precomputations(self, int start) except *: r""" Compute and store `\prod_j=1^{i-1} m_j^{-1} (mod m_i)` for i >= start. """
cdef int min_moduli_count(self, mpz_t height) except -1: r""" Compute the minimum number of primes needed to uniquely determine an integer mod height. """ self._extend_moduli_to_height_c(height)
cdef int count count = self.n * mpz_sizeinbase(height, 2) / mpz_sizeinbase(self.partial_products[self.n-1], 2) # an estimate count = max(min(count, self.n), 1) while count > 1 and mpz_cmp(height, self.partial_products[count-1]) < 0: count -= 1 while mpz_cmp(height, self.partial_products[count-1]) > 0: count += 1
return count
cdef mod_int last_prime(self): return self.moduli[self.n-1]
cdef int mpz_reduce_tail(self, mpz_t z, mod_int* b, int offset, int len) except -1: r""" Perform reduction mod `m_i` for offset <= i < len.
`b[i] = z mod m_{i+offset}` for 0 <= i < len
INPUT:
- ``z`` - the integer being reduced - ``b`` - array to hold the reductions mod each m_i. It MUST be allocated and have length at least len - ``offset`` - first prime in list to reduce against - ``len`` - number of primes in list to reduce against """ cdef int i cdef mod_int* m
cdef int mpz_reduce_vec_tail(self, mpz_t* z, mod_int** b, int vn, int offset, int len) except -1: r""" Perform reduction mod `m_i` for offset <= i < len.
`b[i][j] = z[j] mod m_{i+offset}` for 0 <= i < len
INPUT:
- ``z`` - an array of integers being reduced - ``b`` - array to hold the reductions mod each m_i. It MUST be fully allocated and each have length at least len - ``vn`` - length of z and each b[i] - ``offset`` - first prime in list to reduce against - ``len`` - number of primes in list to reduce against """ cdef int i, j cdef mod_int* m m = self.moduli + offset for i in range(len): mi = m[i] for j in range(vn): b[i][j] = mpz_fdiv_ui(z[j], mi) return 0
cdef int mpz_crt_tail(self, mpz_t z, mod_int* b, int offset, int len) except -1: r""" Calculate lift mod `\prod_{i=0}^{offset+len-1} m_i`.
z = b[i] mod `m_{i+offset}` for 0 <= i < len
In the case that offset > 0, z remains unchanged mod `\prod_{i=0}^{offset-1} m_i`
INPUT:
- ``z`` - a placeholder for the constructed integer z MUST be initialized IF and ONLY IF offset > 0 - ``b`` - array holding the reductions mod each m_i. It MUST have length at least len - ``offset`` - first prime in list to reduce against - ``len`` - number of primes in list to reduce against """ cdef int i, s cdef mpz_t u cdef mod_int* m while s < len and b[s] == 0: # fast forward to first non-zero s += 1 else: s = 0
# normalize to be between -prod/2 and prod/2.
cdef int mpz_crt_vec_tail(self, mpz_t* z, mod_int** b, int vc, int offset, int len) except -1: r""" Calculate lift mod `\prod_{i=0}^{offset+len-1} m_i`.
`z[j] = b[i][j] mod m_{i+offset}` for 0 <= i < len
In the case that offset > 0, z[j] remains unchanged mod `\prod_{i=0}^{offset-1} m_i`
INPUT:
- ``z`` - a placeholder for the constructed integers z MUST be allocated and have length at least vc z[j] MUST be initialized IF and ONLY IF offset > 0 - ``b`` - array holding the reductions mod each m_i. MUST have length at least len - ``vn`` - length of z and each b[i] - ``offset`` - first prime in list to reduce against - ``len`` - number of primes in list to reduce against """ cdef int i, j cdef mpz_t u cdef mod_int* m
else: s = 0
# normalize to be between -prod/2 and prod/2.
cdef Integer zz
def crt(self, b): r""" Calculate lift mod `\prod_{i=0}^{len(b)-1} m_i`.
In the case that offset > 0, z[j] remains unchanged mod `\prod_{i=0}^{offset-1} m_i`
INPUT:
- ``b`` - a list of length at most self.n
OUTPUT:
Integer z where `z = b[i] mod m_i` for 0 <= i < len(b)
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007, 10009, 10037, 10039, 17351]) sage: res = mm.crt([3,5,7,9]); res 8474803647063985 sage: res % 10007 3 sage: res % 10009 5 sage: res % 10037 7 sage: res % 10039 9
""" cdef int i, n raise IndexError("beyond bound for multi-modular prime list") cdef mod_int* bs cdef Integer z
def precomputation_list(self): """ Return a list of the precomputed coefficients `\prod_j=1^{i-1} m_j^{-1} (mod m_i)` where `m_i` are the prime moduli.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307, 10007]); mm MultiModularBasis with moduli [46307, 10007] sage: mm.precomputation_list() [1, 4013] """
def partial_product(self, n): """ Return a list containing precomputed partial products.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307, 10007]); mm MultiModularBasis with moduli [46307, 10007] sage: mm.partial_product(0) 46307 sage: mm.partial_product(1) 463394149
TESTS::
sage: mm.partial_product(2) Traceback (most recent call last): ... IndexError: beyond bound for multi-modular prime list sage: mm.partial_product(-2) Traceback (most recent call last): ... IndexError: negative index not valid
""" cdef Integer z
def prod(self): """ Return the product of the prime moduli.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307]); mm MultiModularBasis with moduli [46307] sage: mm.prod() 46307 sage: mm = MultiModularBasis_base([46307, 10007]); mm MultiModularBasis with moduli [46307, 10007] sage: mm.prod() 463394149
TESTS::
sage: mm = MultiModularBasis_base([]); mm MultiModularBasis with moduli [] sage: len(mm) 0 sage: mm.prod() 1 """ cdef Integer z
def list(self): """ Return a list with the prime moduli.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([46307, 10007]) sage: mm.list() [46307, 10007] """
def __len__(self): """ Returns the number of moduli stored.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007]) sage: len(mm) 1 sage: mm._extend_moduli_to_count(2) 2 sage: len(mm) 2 """
def __iter__(self): """ Return an iterator over the prime moduli.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007, 10009]) sage: t = iter(mm); t <list...iterator object at ...> sage: list(mm.__iter__()) [10007, 10009] """
def __getitem__(self, ix): """ Return the moduli stored at index `ix` as a Python long.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: mm = MultiModularBasis_base([10007, 10009]) sage: mm[1] 10009 # 64-bit 10009L # 32-bit sage: mm[-1] Traceback (most recent call last): ... IndexError: index out of range
sage: mm[:1] MultiModularBasis with moduli [10007] """
raise TypeError("index must be an integer")
def __repr__(self): """ Return a string representation of this object.
EXAMPLES::
sage: from sage.arith.multi_modular import MultiModularBasis_base sage: MultiModularBasis_base([10007]) MultiModularBasis with moduli [10007] """
cdef class MultiModularBasis(MultiModularBasis_base): """ Class used for storing a MultiModular bases of a fixed length. """ cdef int mpz_reduce(self, mpz_t z, mod_int* b) except -1: r""" Perform reduction mod `m_i` for each modulus `m_i`.
`b[i] = z mod m_i` for 0 <= i < len(self)
INPUT:
- ``z`` -- the integer being reduced - ``b`` -- array to hold the reductions mod each m_i. It MUST be allocated and have length at least len """
cdef int mpz_reduce_vec(self, mpz_t* z, mod_int** b, int vn) except -1: r""" Perform reduction mod `m_i` for each modulus `m_i`.
`b[i][j] = z[j] mod m_i` for 0 <= i < len(self)
INPUT:
- ``z`` -- an array of integers being reduced - ``b`` -- array to hold the reductions mod each m_i. It MUST be fully allocated and each have length at least len - ``vn`` -- length of z and each b[i] """ self.mpz_reduce_vec_tail(z, b, vn, 0, self.n)
cdef int mpz_crt(self, mpz_t z, mod_int* b) except -1: r""" Calculate lift mod `\prod m_i`.
`z = b[i] mod m_{i+offset}` for 0 <= i < len(self)
INPUT:
- ``z`` -- a placeholder for the constructed integer z MUST NOT be initialized - ``b`` -- array holding the reductions mod each `m_i`. It MUST have length at least len(self) """ self.mpz_crt_tail(z, b, 0, self.n)
cdef int mpz_crt_vec(self, mpz_t* z, mod_int** b, int vn) except -1: r""" Calculate lift mod `\prod m_i`.
`z[j] = b[i][j] mod m_i` for 0 <= i < len(self)
INPUT:
- ``z`` -- a placeholder for the constructed integers z MUST be allocated and have length at least vn, but each z[j] MUST NOT be initialized - ``b`` -- array holding the reductions mod each `m_i`. It MUST have length at least len(self) - ``vn`` -- length of z and each b[i] """
cdef class MutableMultiModularBasis(MultiModularBasis): """ Class used for performing multi-modular methods, with the possibility of removing bad primes. """ cpdef mod_int next_prime(self) except -1: """ Pick a new random prime between the bounds given during the initialization of this object, update the precomputed data, and return the new prime modulus.
EXAMPLES::
sage: from sage.arith.multi_modular import MutableMultiModularBasis sage: mm = MutableMultiModularBasis([10007]) sage: mm.next_prime() 4561 # 64-bit 4561L # 32-bit sage: mm MultiModularBasis with moduli [10007, 4561] """
cpdef mod_int replace_prime(self, int ix) except -1: """ Replace the prime moduli at the given index with a different one, update the precomputed data accordingly, and return the new prime.
INPUT:
- ``ix`` -- index into list of moduli
OUTPUT: the new prime modulus
EXAMPLES::
sage: from sage.arith.multi_modular import MutableMultiModularBasis sage: mm = MutableMultiModularBasis([10007, 10009, 10037, 10039]) sage: mm MultiModularBasis with moduli [10007, 10009, 10037, 10039] sage: mm.prod() 10092272478850909 sage: mm.precomputation_list() [1, 5004, 6536, 6060] sage: mm.partial_product(2) 1005306552331 sage: mm.replace_prime(1) 4561 # 64-bit 4561L # 32-bit sage: mm MultiModularBasis with moduli [10007, 4561, 10037, 10039] sage: mm.prod() 4598946425820661 sage: mm.precomputation_list() [1, 2314, 3274, 3013] sage: mm.partial_product(2) 458108021299 """ cdef mod_int new_p
raise IndexError("index out of range")
|