"""
Utility classes for multi-modular algorithms.
"""
include "../ext/interrupt.pxi"
include "../ext/stdsage.pxi"
include "../ext/gmp.pxi"
from sage.rings.integer_ring import ZZ
from sage.rings.arith import random_prime
from types import GeneratorType
from sage.rings.fast_arith cimport arith_llong
cdef arith_llong ai
ai = arith_llong()
MAX_MODULUS = 2**23
cdef class MultiModularBasis_base:
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 is generated automatically in the second case::
sage: from sage.ext.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]
"""
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.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([1099511627791])
Traceback (most recent call last):
...
ValueError: given modulus cannot be manipulated in a single machine word
"""
mpz_init(self.product)
mpz_init(self.half_product)
cdef _realloc_to_new_count(self, new_count):
self.moduli = <mod_int*>sage_realloc(self.moduli, sizeof(mod_int)*new_count)
self.partial_products = <mpz_t*>sage_realloc(self.partial_products, sizeof(mpz_t)*new_count)
self.C = <mod_int*>sage_realloc(self.C, sizeof(mod_int)*new_count)
if self.moduli == NULL or self.partial_products == NULL or self.C == NULL:
raise MemoryError, "out of memory allocating multi-modular prime list"
def __dealloc__(self):
"""
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([1099511627791])
Traceback (most recent call last):
...
ValueError: given modulus cannot be manipulated in a single machine word
sage: mm = MultiModularBasis_base(1099511627791); mm
MultiModularBasis with moduli [4561, 17351, 28499]
sage: del mm
"""
sage_free(self.moduli)
for i from 0 <= i < self.n:
mpz_clear(self.partial_products[i])
sage_free(self.partial_products)
sage_free(self.C)
mpz_clear(self.product)
mpz_clear(self.half_product)
def __init__(self, val, l_bound=0, u_bound=0):
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.ext.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):
...
ValueError: given modulus cannot be manipulated in a single machine word
sage: mm = MultiModularBasis_base(0); mm
MultiModularBasis with moduli [28499]
sage: mm = MultiModularBasis_base([6, 10])
Traceback (most recent call last):
...
ValueError: the values provided are not relatively prime
"""
if l_bound == 0:
self._l_bound = 2**10
elif l_bound < 2:
raise ValueError("minimum value for lower bound is 2, given: %s"%(l_bound))
else:
self._l_bound = l_bound
if u_bound == 0:
self._u_bound = 2**15
elif u_bound > MAX_MODULUS:
raise ValueError("upper bound cannot be greater than MAX_MODULUS: %s, given: %s"%(MAX_MODULUS, u_bound))
else:
self._u_bound = u_bound
from sage.functions.prime_pi import prime_pi
self._num_primes = prime_pi(self._u_bound) - prime_pi(self._l_bound-1)
cdef int i
if isinstance(val, (list, tuple, GeneratorType)):
try:
self.extend_with_primes(val, check=True)
except ArithmeticError:
raise ValueError("the values provided are not relatively prime")
else:
self._extend_moduli_to_height(val)
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
while True:
if len(known_primes) >= self._num_primes:
raise RuntimeError, "there are not enough primes in the interval [%s, %s] to complete this multimodular computation"%(self._l_bound, self._u_bound)
p = random_prime(self._u_bound, lbound =self._l_bound)
if p not in known_primes:
return 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.ext.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]
"""
if isinstance(plist, GeneratorType):
plist = list(plist)
elif not isinstance(plist, (tuple, list)):
raise ValueError, "plist should be a list, tuple or a generator, got: %s"%(str(type(plist)))
cdef Py_ssize_t len_plist = len(plist)
if len_plist == 0:
return self.n
if check:
for p in plist:
if p > MAX_MODULUS:
raise ValueError, "given modulus cannot be manipulated in a single machine word"
self._realloc_to_new_count(self.n + len_plist)
cdef Py_ssize_t i
for i in range(len_plist):
self.moduli[self.n+i] = plist[i]
mpz_init(self.partial_products[self.n + i])
if partial_products:
mpz_set(self.partial_products[self.n + i],
(<Integer>partial_products[i]).value)
cdef int old_count = self.n
self.n += len_plist
if not partial_products:
self._refresh_products(old_count)
else:
self._refresh_prod()
self._refresh_precomputations(old_count)
return self.n
def __cmp__(self, other):
"""
EXAMPLES::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([10007])
sage: nn = MultiModularBasis_base([10007])
sage: mm == nn
True
sage: mm == 1
False
"""
if not PY_TYPE_CHECK(other, MultiModularBasis_base):
return cmp(type(other), MultiModularBasis_base)
return cmp((self.list(), self._u_bound, self._l_bound),
(other.list(), (<MultiModularBasis_base>other)._u_bound,
(<MultiModularBasis_base>other)._l_bound))
def __setstate__(self, state):
"""
Initialize a new MultiModularBasis_base object from a state stored
in a pickle.
TESTS::
sage: from sage.ext.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]
"""
if len(state) != 3:
raise ValueError, "cannot read state tuple"
nmoduli, lbound, ubound = state
self._realloc_to_new_count(len(nmoduli))
self._l_bound = lbound
self._u_bound = ubound
self.extend_with_primes(nmoduli, check=False)
def __getstate__(self):
"""
Return a tuple describing the state of this object for pickling.
TESTS::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([10007, 10009])
sage: mm.__getstate__()
([10007, 10009], 1024L, 32768L)
"""
return (self.list(), self._l_bound, self._u_bound)
def _extend_moduli_to_height(self, height):
"""
EXAMPLES::
sage: from sage.ext.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.ext.multi_modular.MultiModularBasis_base(0)
sage: m._extend_moduli_to_height(prod(prime_range(50)))
sage: m = sage.ext.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 #11358 is fixed)::
sage: set_random_seed(0); m = sage.ext.multi_modular.MultiModularBasis_base(0); m._extend_moduli_to_height(10**10000)
sage: len(set(m)) == len(m)
True
sage: len(m)
2438
"""
cdef Integer h
h = ZZ(height)
if h < self._l_bound:
h = ZZ(self._l_bound)
self._extend_moduli_to_height_c(h.value)
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)
"""
cdef mpz_t height
mpz_init(height)
mpz_mul_2exp(height, height0, 1)
if self.n > 0 and mpz_cmp(height, self.partial_products[self.n-1]) <= 0:
mpz_clear(height)
return self.n
cdef int i
new_moduli = []
new_partial_products = []
cdef Integer M
cdef mod_int p
if self.n == 0:
M = ZZ(1)
else:
M = PY_NEW(Integer)
mpz_set(M.value, self.partial_products[self.n-1])
known_primes = set(self)
while mpz_cmp(height, M.value) > 0:
p = self._new_random_prime(known_primes)
new_moduli.append(p)
known_primes.add(p)
M *= p
new_partial_products.append(M)
mpz_clear(height)
return self.extend_with_primes(new_moduli, new_partial_products,
check=False)
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.ext.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]
"""
if count <= self.n:
return self.n
new_moduli = []
cdef int i
cdef mod_int p
known_primes = set(self)
for i from self.n <= i < count:
p = self._new_random_prime(known_primes)
known_primes.add(p)
new_moduli.append(p)
return self.extend_with_primes(new_moduli, check=False)
def _extend_moduli(self, count):
"""
Expand the list of prime moduli with `count` new random primes.
EXAMPLES::
sage: from sage.ext.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]
"""
self._extend_moduli_to_count(self.n + count)
cdef void _refresh_products(self, int start):
r"""
Compute and store $\prod_j=1^{i-1} m_j$ of i > start.
"""
cdef mpz_t z
mpz_init(z)
if start == 0:
mpz_set_ui(self.partial_products[0], self.moduli[0])
start += 1
for i from start <= i < self.n:
mpz_set_ui(z, self.moduli[i])
mpz_mul(self.partial_products[i], self.partial_products[i-1], z)
mpz_clear(z)
self._refresh_prod()
cdef void _refresh_prod(self):
mpz_set(self.product, self.partial_products[self.n-1])
mpz_fdiv_q_ui(self.half_product, self.product, 2)
cdef void _refresh_precomputations(self, int start) except *:
r"""
Compute and store $\prod_j=1^{i-1} m_j^{-1} (mod m_i)$ of i >= start.
"""
if start == 0:
start = 1
self.C[0] = 1
for i from start <= i < self.n:
self.C[i] = ai.c_inverse_mod_longlong(mpz_fdiv_ui(self.partial_products[i-1], self.moduli[i]), self.moduli[i])
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)
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"""
Performs 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
m = self.moduli + offset
for i from 0 <= i < len:
b[i] = mpz_fdiv_ui(z, m[i])
return 0
cdef int mpz_reduce_vec_tail(self, mpz_t* z, mod_int** b, int vn, int offset, int len) except -1:
r"""
Performs 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 from 0 <= i < len:
mi = m[i]
for j from 0 <= j < 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
m = self.moduli + offset
mpz_init(u)
if offset == 0:
s = 1
mpz_init_set_ui(z, b[0])
if b[0] == 0:
while s < len and b[s] == 0:
s += 1
else:
s = 0
for i from s <= i < len:
mpz_set_ui(u, ((b[i] + m[i] - mpz_fdiv_ui(z, m[i])) * self.C[i]) % m[i])
mpz_mul(u, u, self.partial_products[i-1])
mpz_add(z, z, u)
if mpz_cmp(z, self.half_product) > 0:
mpz_sub(z, z, self.product)
mpz_clear(u)
return 0
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
m = self.moduli + offset
mpz_init(u)
if offset == 0:
s = 1
else:
s = 0
for j from 0 <= j < vc:
i = s
if offset == 0:
mpz_set_ui(z[j], b[0][j])
if b[0][j] == 0:
while i < len and b[i][j] == 0:
i += 1
while i < len:
mpz_set_ui(u, ((b[i][j] + m[i] - mpz_fdiv_ui(z[j], m[i])) * self.C[i]) % m[i])
mpz_mul(u, u, self.partial_products[i-1])
mpz_add(z[j], z[j], u)
i += 1
if mpz_cmp(z[j], self.half_product) > 0:
mpz_sub(z[j], z[j], self.product)
cdef Integer zz
zz = PY_NEW(Integer)
mpz_set(zz.value, self.half_product)
mpz_clear(u)
return 0
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.ext.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
n = len(b)
if n > self.n:
raise IndexError, "beyond bound for multi-modular prime list"
cdef mod_int* bs
bs = <mod_int*>sage_malloc(sizeof(mod_int)*n)
if bs == NULL:
raise MemoryError, "out of memory allocating multi-modular prime list"
for i from 0 <= i < n:
bs[i] = b[i]
cdef Integer z
z = PY_NEW(Integer)
self.mpz_crt_tail(z.value, bs, 0, n)
sage_free(bs)
return z
def precomputation_list(self):
"""
Returns 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.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([46307, 10007]); mm
MultiModularBasis with moduli [46307, 10007]
sage: mm.precomputation_list()
[1, 4013]
"""
cdef int i
return [ZZ(self.C[i]) for i from 0 <= i < self.n]
def partial_product(self, n):
"""
Returns a list containing precomputed partial products.
EXAMPLES::
sage: from sage.ext.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
"""
if n >= self.n:
raise IndexError, "beyond bound for multi-modular prime list"
if n < 0:
raise IndexError, "negative index not valid"
cdef Integer z
z = PY_NEW(Integer)
mpz_set(z.value, self.partial_products[n])
return z
def prod(self):
"""
Returns the product of the prime moduli.
EXAMPLES::
sage: from sage.ext.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
"""
if self.n == 0:
return 1
cdef Integer z
z = PY_NEW(Integer)
mpz_set(z.value, self.partial_products[self.n-1])
return z
def list(self):
"""
Return a list with the prime moduli.
EXAMPLES::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([46307, 10007])
sage: mm.list()
[46307, 10007]
"""
cdef Py_ssize_t i
return [ZZ(self.moduli[i]) for i in range(self.n)]
def __len__(self):
"""
Returns the number of moduli stored.
EXAMPLES::
sage: from sage.ext.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
"""
return self.n
def __iter__(self):
"""
Returns an iterator over the prime moduli.
EXAMPLES::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([10007, 10009])
sage: t = iter(mm); t
<listiterator object at ...>
sage: list(mm.__iter__())
[10007, 10009]
"""
return self.list().__iter__()
def __getitem__(self, ix):
"""
Return the moduli stored at index `ix` as a Python long.
EXAMPLES::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: mm = MultiModularBasis_base([10007, 10009])
sage: mm[1]
10009L
sage: mm[-1]
Traceback (most recent call last):
...
IndexError: index out of range
#sage: mm[:1]
MultiModularBasis with moduli [10007]
"""
if isinstance(ix, slice):
return self.__class__(self.list()[ix], l_bound = self._l_bound,
u_bound = self._u_bound)
cdef int i
i = ix
if i != ix:
raise ValueError, "index must be an integer"
if i < 0 or i >= self.n:
raise IndexError, "index out of range"
return self.moduli[i]
def __repr__(self):
"""
Return a string representation of this object.
EXAMPLES::
sage: from sage.ext.multi_modular import MultiModularBasis_base
sage: MultiModularBasis_base([10007])
MultiModularBasis with moduli [10007]
"""
return "MultiModularBasis with moduli "+str(self.list())
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"""
Performs 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
"""
self.mpz_reduce_tail(z, b, 0, self.n)
cdef int mpz_reduce_vec(self, mpz_t* z, mod_int** b, int vn) except -1:
r"""
Performs 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]
"""
self.mpz_crt_vec_tail(z, b, vn, 0, self.n)
cdef class MutableMultiModularBasis(MultiModularBasis):
"""
Class used for performing multi-modular methods,
with the possibility of removing bad primes.
"""
def next_prime(self):
"""
Picks a new random prime between the bounds given during the
initialization of this object, updates the precomputed data,
and returns the new prime modulus.
EXAMPLES::
sage: from sage.ext.multi_modular import MutableMultiModularBasis
sage: mm = MutableMultiModularBasis([10007])
sage: mm.next_prime()
4561L
sage: mm
MultiModularBasis with moduli [10007, 4561]
"""
return self.next_prime_c()
cdef mod_int next_prime_c(self) except -1:
self._extend_moduli(1)
return self.moduli[self.n-1]
def replace_prime(self, ix):
"""
Replace the prime moduli at the given index with a different one,
update the precomputed data accordingly, and return the new prime.
EXAMPLES::
sage: from sage.ext.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)
4561L
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
"""
return self.replace_prime_c(ix)
cdef mod_int replace_prime_c(self, int ix) except -1:
r"""
Replace $m_{ix}$ in the list of moduli with a new
prime number greater than all others in the list,
and recompute all precomputations.
INPUT:
- ``ix`` - index into list of moduli
OUTPUT:
p -- the new prime modulus
"""
cdef mod_int new_p
if ix < 0 or ix >= self.n:
raise IndexError, "index out of range"
new_p = self._new_random_prime(set(self))
self.moduli[ix] = new_p
self._refresh_products(ix)
self._refresh_precomputations(ix)
return new_p
