source: sage/ext/multi_modular.pyx @ 4878:5ab27888269a

Revision 4878:5ab27888269a, 19.8 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

Fix bug that Thea Gegenberg reported in vectors mod n.
Line 
1"""
2Utility classes for multi-modular algorithms.
3"""
4
5######################################################################
6#       Copyright (C) 2006 William Stein
7#
8#  Distributed under the terms of the GNU General Public License (GPL)
9#
10#  The full text of the GPL is available at:
11#                  http://www.gnu.org/licenses/
12######################################################################
13
14
15include "../ext/interrupt.pxi"
16include "../ext/stdsage.pxi"
17include "../ext/gmp.pxi"
18
19
20from sage.rings.integer_ring import ZZ
21from sage.rings.arith import next_prime # should I just use probable primes?
22
23# should I have mod_int versions of these functions?
24# c_inverse_mod_longlong modular inverse used exactly once in _refresh_precomputations
25from sage.ext.arith cimport arith_llong
26cdef arith_llong ai
27ai = arith_llong()
28
29# We use both integer and double operations, hence the min.
30# MAX_MODULUS = min(int(sqrt(int(MOD_INT_OVERFLOW))-1), int(2)**int(20))
31
32# Hard coded because currently matrix_modn_dense is implemented using C ints
33# which are always 32-bit.   Once this gets firxed, i.e., there is a better
34# matrix_modn class, then this can change.
35MAX_MODULUS = 46341
36
37# TODO: have one global instance for sharing, copy for MutableMultiModularBasis
38cdef class MultiModularBasis_base:
39    r"""
40    This class stores a list of machine-sized prime numbers,
41    and can do reduction and Chinese Remainder Theorem lifting
42    modulo these primes.
43   
44    Lifting implemented via Garner's algorithm, which has the advantage
45    that all reductions are word-sized. For each $i$ precompute
46       
47       $\prod_j=1^{i-1} m_j$ and $\prod_j=1^{i-1} m_j^{-1} (mod m_i)$
48    """
49
50    def __new__(self, height):
51        r"""
52        Allocate the space for the moduli and precomputation lists
53        and initalize the first element of that list.
54        """
55        cdef mod_int p
56        p = next_prime(MAX_MODULUS/7)
57
58        self.n = 1
59
60        self.moduli = <mod_int*>sage_malloc(sizeof(mod_int))
61        self.partial_products = <mpz_t*>sage_malloc(sizeof(mpz_t))
62        self.C = <mod_int*>sage_malloc(sizeof(mod_int))
63
64        if self.moduli == NULL or self.partial_products == NULL or self.C == NULL:
65            raise MemoryError, "out of memory allocating multi-modular prime list"
66
67        self.moduli[0] = p
68        mpz_init_set_ui(self.partial_products[0], p)
69        self.C[0] = 1
70
71        mpz_init(self.product)
72        mpz_init(self.half_product)
73       
74    def __dealloc__(self):
75        sage_free(self.moduli)
76        sage_free(self.partial_products)
77        sage_free(self.C)
78
79    def __init__(self, height):
80        r"""
81        Initialize a multi-modular basis and perform precomputations.
82       
83        INPUT:
84            height -- as integer
85                          determines how many primes are computed
86                          (their product will be at least 2*height)
87                      as list
88                          a list of prime moduli to start with
89        """
90        cdef int i
91        if isinstance(height, list):
92            m = height
93            self.n = len(m)
94            self.moduli = <mod_int*>sage_realloc(self.moduli, sizeof(mod_int) * self.n)
95            self.partial_products = <mpz_t*>sage_realloc(self.partial_products, sizeof(mpz_t) * self.n)
96            self.C = <mod_int*>sage_realloc(self.C, sizeof(mod_int) * self.n)
97            if self.moduli == NULL or self.partial_products == NULL or self.C == NULL:
98                raise MemoryError, "out of memory allocating multi-modular prime list"
99            for i from 0 <= i < len(m):
100                if m[i] > MAX_MODULUS:
101                    raise ValueError, "given modulus cannot be manipulated in a single machine word"
102                self.moduli[i] = m[i]
103            for i from 1 <= i < len(m):
104                mpz_init(self.partial_products[i])
105            self._refresh_products(0)
106            self._refresh_precomputations(0)
107        else:
108            self._extend_moduli_to_height(height)
109
110        self._refresh_prod()       
111       
112       
113    def _extend_moduli_to_height(self, height):
114        cdef Integer h
115        h = ZZ(height)
116        self._extend_moduli_to_height_c(h.value)
117   
118    cdef int _extend_moduli_to_height_c(self, mpz_t height0) except -1:
119        r"""
120        Expand the list of primes and perform precomputations.
121       
122        INPUT:
123            height -- determines how many primes are computed
124                      (their product must be at least height)
125        """
126        cdef mpz_t height
127        mpz_init(height)
128        mpz_mul_2exp(height, height0, 1)
129        if mpz_cmp(height, self.partial_products[self.n-1]) <= 0:
130            return self.n
131        cdef int i
132        new_moduli = []
133        new_partial_products = []
134        cdef Integer p, M
135        M = PY_NEW(Integer)
136        mpz_set(M.value, self.partial_products[self.n-1])
137        p = ZZ(self.last_prime())
138        while mpz_cmp(height, M.value) > 0:
139            p = next_prime(p)
140            new_moduli.append(p)
141            M *= p
142            new_partial_products.append(M)
143        mpz_clear(height)
144           
145        cdef int new_count, old_count
146        old_count = self.n
147        new_count = self.n + len(new_moduli)
148
149        self.moduli = <mod_int*>sage_realloc(self.moduli, sizeof(mod_int)*new_count)
150        self.partial_products = <mpz_t*>sage_realloc(self.partial_products, sizeof(mpz_t)*new_count)
151        self.C = <mod_int*>sage_realloc(self.C, sizeof(mod_int)*new_count)
152        if self.moduli == NULL or self.partial_products == NULL or self.C == NULL:
153            raise MemoryError, "out of memory allocating multi-modular prime list"
154       
155        for i from self.n <= i < new_count:
156            self.moduli[i] = new_moduli[i-self.n]
157            mpz_init_set(self.partial_products[i], (<Integer>new_partial_products[i-self.n]).value)
158           
159        self.n = new_count
160        self._refresh_precomputations(old_count)
161        return new_count
162       
163    def _extend_moduli_to_count(self, int count):
164        r"""
165        Expand the list of primes and perform precomputations.
166       
167        INPUT:
168            count -- the minimum number of moduli in the resulting list
169        """
170        if count <= self.n:
171            return self.n
172        self.moduli = <mod_int*>sage_realloc(self.moduli, sizeof(mod_int) * count)
173        self.partial_products = <mpz_t*>sage_realloc(self.partial_products, sizeof(mpz_t) * count)
174        self.C = <mod_int*>sage_realloc(self.C, sizeof(mod_int) * count)
175        if self.moduli == NULL or self.partial_products == NULL or self.C == NULL:
176            raise MemoryError, "out of memory allocating multi-modular prime list"
177           
178        cdef int i
179        cdef Integer p
180        p = ZZ(self.last_prime())
181        for i from self.n <= i < count:
182            p = next_prime(p)
183            self.moduli[i] = p
184            mpz_init(self.partial_products[i])
185            mpz_mul(self.partial_products[i], self.partial_products[i-1], p.value)
186        old_count = self.n
187        self.n = count
188        self._refresh_precomputations(old_count)
189        return count
190
191    def _extend_moduli(self, count):
192        self._extend_moduli_to_count(self.n + count)
193
194    cdef void _refresh_products(self, int start):
195        r"""
196        Compute and store $\prod_j=1^{i-1} m_j$ of i > start.
197        """
198        cdef mpz_t z
199        mpz_init(z)
200        if start == 0:
201            mpz_set_ui(self.partial_products[0], self.moduli[0])
202            start += 1
203        for i from start <= i < self.n:
204            mpz_set_ui(z, self.moduli[i])
205            mpz_mul(self.partial_products[i], self.partial_products[i-1], z)
206        mpz_clear(z)
207        self._refresh_prod()
208
209    cdef void _refresh_prod(self):
210        # record the product and half product for balancing the lifts.
211        mpz_set(self.product, self.partial_products[self.n-1])
212        mpz_fdiv_q_ui(self.half_product, self.product, 2)
213       
214    cdef void _refresh_precomputations(self, int start):
215        r"""
216        Compute and store $\prod_j=1^{i-1} m_j^{-1} (mod m_i)$ of i >= start.
217        """
218        if start == 0:
219            start = 1 # first one is trivial, never used
220        for i from start <= i < self.n:
221            self.C[i] = ai.c_inverse_mod_longlong(mpz_fdiv_ui(self.partial_products[i-1], self.moduli[i]), self.moduli[i])
222       
223    cdef int min_moduli_count(self, mpz_t height) except -1:
224        r"""
225        Compute the minimum number of primes needed to uniquely determin
226        an integer mod height.
227        """
228        self._extend_moduli_to_height_c(height)
229       
230        cdef int count
231        count = self.n * mpz_sizeinbase(height, 2) / mpz_sizeinbase(self.partial_products[self.n-1], 2) # an estimate
232        count = max(min(count, self.n), 1)
233        while count > 1 and mpz_cmp(height, self.partial_products[count-1]) < 0:
234           count -= 1
235        while mpz_cmp(height, self.partial_products[count-1]) > 0:
236           count += 1
237           
238        return count
239
240    cdef int moduli_list_c(self, mod_int** moduli):
241        memcpy(moduli[0], self.moduli, sizeof(mod_int)*self.n)
242        return self.n
243       
244    cdef mod_int last_prime(self):
245        return self.moduli[self.n-1]
246       
247    cdef int mpz_reduce_tail(self, mpz_t z, mod_int* b, int offset, int len) except -1:
248        r"""
249        Performs reduction mod $m_i$ for offset <= i < len
250       
251        b[i] = z mod $m_{i+offset}$ for 0 <= i < len
252       
253        INPUT:
254            z -- the integer being reduced
255            b -- array to hold the reductions mod each m_i.
256                 It MUST be allocated and have length at least len
257            offset -- first prime in list to reduce against
258            len    -- number of primes in list to reduce against
259        """
260        cdef int i
261        cdef mod_int* m
262        m = self.moduli + offset
263        for i from 0 <= i < len:
264            b[i] = mpz_fdiv_ui(z, m[i])
265        return 0
266       
267    cdef int mpz_reduce_vec_tail(self, mpz_t* z, mod_int** b, int vn, int offset, int len) except -1:
268        r"""
269        Performs reduction mod $m_i$ for offset <= i < len
270       
271        b[i][j] = z[j] mod $m_{i+offset}$ for 0 <= i < len
272       
273        INPUT:
274            z  -- an array of integers being reduced
275            b  -- array to hold the reductions mod each m_i.
276                 It MUST be fully allocated and each
277                 have length at least len
278            vn -- length of z and each b[i]
279            offset -- first prime in list to reduce against
280            len    -- number of primes in list to reduce against
281        """
282        cdef int i, j
283        cdef mod_int* m
284        m = self.moduli + offset
285        for i from 0 <= i < len:
286            mi = m[i]
287            for j from 0 <= j < vn:
288                b[i][j] = mpz_fdiv_ui(z[j], mi)
289        return 0
290       
291    cdef int mpz_crt_tail(self, mpz_t z, mod_int* b, int offset, int len) except -1:
292        r"""
293        Calculate lift mod $\prod_{i=0}^{offset+len-1} m_i$.
294       
295        z = b[i] mod $m_{i+offset}$ for 0 <= i < len
296       
297        In the case that offset > 0,
298        z remains unchanged mod $\prod_{i=0}^{offset-1} m_i$
299
300        INPUT:
301            z  -- a placeholder for the constructed integer
302                  z MUST be initalized IF and ONLY IF offset > 0
303            b  -- array holding the reductions mod each m_i.
304                  It MUST have length at least len
305            offset -- first prime in list to reduce against
306            len    -- number of primes in list to reduce against
307        """
308        cdef int i, s
309        cdef mpz_t u
310        cdef mod_int* m
311        m = self.moduli + offset
312        mpz_init(u)
313        if offset == 0:
314            s = 1
315            mpz_init_set_ui(z, b[0])
316            if b[0] == 0:
317                while s < len and b[s] == 0: # fast forward to first non-zero
318                    s += 1
319        else:
320            s = 0
321        for i from s <= i < len:
322            mpz_set_ui(u, ((b[i] + m[i] - mpz_fdiv_ui(z, m[i])) * self.C[i]) % m[i])
323            mpz_mul(u, u, self.partial_products[i-1])
324            mpz_add(z, z, u)
325           
326        # normalize to be between -prod/2 and prod/2.
327        if mpz_cmp(z, self.half_product) > 0:
328            mpz_sub(z, z, self.product)
329        mpz_clear(u)
330        return 0
331
332    cdef int mpz_crt_vec_tail(self, mpz_t* z, mod_int** b, int vc, int offset, int len) except -1:
333        r"""
334        Calculate lift mod $\prod_{i=0}^{offset+len-1} m_i$.
335       
336        z[j] = b[i][j] mod $m_{i+offset}$ for 0 <= i < len
337       
338        In the case that offset > 0,
339        z[j] remains unchanged mod $\prod_{i=0}^{offset-1} m_i$
340
341        INPUT:
342            z  -- a placeholder for the constructed integers
343                  z MUST be allocated and have length at least vc
344                  z[j] MUST be initalized IF and ONLY IF offset > 0
345            b  -- array holding the reductions mod each m_i.
346                  MUST have length at least len
347            vn -- length of z and each b[i]
348            offset -- first prime in list to reduce against
349            len    -- number of primes in list to reduce against
350        """
351        cdef int i, j
352        cdef mpz_t u
353        cdef mod_int* m
354
355        m = self.moduli + offset
356        mpz_init(u)
357        if offset == 0:
358            s = 1
359        else:
360            s = 0
361
362        for j from 0 <= j < vc:
363            i = s
364            if offset == 0:
365                mpz_init_set_ui(z[j], b[0][j])
366                if b[0][j] == 0:
367                    while i < len and b[i][j] == 0: # fast forward to first non-zero
368                        i += 1
369            while i < len:
370                mpz_set_ui(u, ((b[i][j] + m[i] - mpz_fdiv_ui(z[j], m[i])) * self.C[i]) % m[i]) # u = ((b_i - z) * C_i) % m_i
371                mpz_mul(u, u, self.partial_products[i-1])
372                mpz_add(z[j], z[j], u)
373                i += 1
374
375            # normalize to be between -prod/2 and prod/2.
376            if mpz_cmp(z[j], self.half_product) > 0:
377                mpz_sub(z[j], z[j], self.product)
378               
379
380        cdef Integer zz
381        zz = PY_NEW(Integer)
382        mpz_init_set(zz.value, self.product)
383        mpz_set(zz.value, self.half_product)
384   
385        mpz_clear(u)
386        return 0
387       
388    def crt(self, b):
389        r"""
390        Calculate lift mod $\prod_{i=0}^{len(b)-1} m_i$.
391       
392        In the case that offset > 0,
393        z[j] remains unchanged mod $\prod_{i=0}^{offset-1} m_i$
394
395        INPUT:
396            b  -- a list of length at most self.n
397        OUTPUT:
398            Integer z where z = b[i] mod $m_i$ for 0 <= i < len(b)
399        """
400        cdef int i, n
401        n = len(b)
402        if n > self.n:
403            raise IndexError, "beyond bound for multi-modular prime list"
404        cdef mod_int* bs
405        bs = <mod_int*>sage_malloc(sizeof(mod_int)*n)
406        if bs == NULL:
407            raise MemoryError, "out of memory allocating multi-modular prime list"
408        for i from 0 <= i < n:
409            bs[i] = b[i]
410        cdef Integer z
411        z = PY_NEW(Integer)
412        self.mpz_crt_tail(z.value, bs, 0, n)
413        sage_free(bs)
414        return z
415
416    def precomputation_list(self):
417        cdef int i
418        return [ZZ(self.C[i]) for i from 0 <= i < self.n]
419       
420    def partial_product(self, n):
421        if n > self.n:
422            raise IndexError, "beyond bound for multi-modular prime list"
423        cdef Integer z
424        z = PY_NEW(Integer)
425        mpz_init_set(z.value, self.partial_products[n-1])
426        return z
427       
428    def prod(self):
429        cdef Integer z
430        z = PY_NEW(Integer)
431        mpz_init_set(z.value, self.partial_products[self.n-1])
432        return z
433           
434    def list(self):
435        cdef int i
436        cdef mod_int* moduli
437        moduli = <mod_int*>sage_malloc(sizeof(mod_int) * self.n)
438        n = self.moduli_list_c(&moduli)
439        m = [ZZ(moduli[i]) for i from 0 <= i < n]
440        sage_free(moduli)
441        return m
442       
443    def __len__(self):
444        return self.n
445       
446    def __iter__(self):
447        return self.list().__iter__()
448       
449    def __getitem__(self, ix):
450   
451        if isinstance(ix, slice):
452            return type(self)(self.list()[ix])
453   
454        cdef int i
455        i = ix
456        if i < 0 or i >= self.n:
457            raise IndexError, "index out of range"
458        return self.moduli[i]
459   
460    def __repr__(self):
461        return str(type(self))+str(self.list())
462       
463       
464cdef class MultiModularBasis(MultiModularBasis_base):
465    """
466    Class used for storing a MultiModular bases of a fixed length.
467    """
468       
469       
470    cdef int mpz_reduce(self, mpz_t z, mod_int* b) except -1:
471        r"""
472        Performs reduction mod $m_i$ for each modulus $m_i$
473       
474        b[i] = z mod $m_i$ for 0 <= i < len(self)
475       
476        INPUT:
477            z -- the integer being reduced
478            b -- array to hold the reductions mod each m_i.
479                 It MUST be allocated and have length at least len
480        """
481        self.mpz_reduce_tail(z, b, 0, self.n)
482       
483    cdef int mpz_reduce_vec(self, mpz_t* z, mod_int** b, int vn) except -1:
484        r"""
485        Performs reduction mod $m_i$ for each modulus $m_i$
486       
487        b[i][j] = z[j] mod $m_i$ for 0 <= i < len(self)
488       
489        INPUT:
490            z  -- an array of integers being reduced
491            b  -- array to hold the reductions mod each m_i.
492                 It MUST be fully allocated and each
493                 have length at least len
494            vn -- length of z and each b[i]
495        """
496        self.mpz_reduce_vec_tail(z, b, vn, 0, self.n)
497       
498    cdef int mpz_crt(self, mpz_t z, mod_int* b) except -1:
499        r"""
500        Calculate lift mod $\prod m_i$.
501       
502        z = b[i] mod $m_{i+offset}$ for 0 <= i < len(self)
503       
504        INPUT:
505            z  -- a placeholder for the constructed integer
506                  z MUST NOT be initalized
507            b  -- array holding the reductions mod each $m_i$.
508                  It MUST have length at least len(self)
509        """
510        self.mpz_crt_tail(z, b, 0, self.n)
511       
512    cdef int mpz_crt_vec(self, mpz_t* z, mod_int** b, int vn) except -1:
513        r"""
514        Calculate lift mod $\prod m_i$.
515       
516        z[j] = b[i][j] mod $m_i$ for 0 <= i < len(self)
517       
518        INPUT:
519            z  -- a placeholder for the constructed integers
520                  z MUST be allocated and have length at least vn,
521                  but each z[j] MUST NOT be initalized
522            b  -- array holding the reductions mod each $m_i$.
523                  It MUST have length at least len(self)
524            vn -- length of z and each b[i]
525        """
526        self.mpz_crt_vec_tail(z, b, vn, 0, self.n)
527       
528       
529cdef class MutableMultiModularBasis(MultiModularBasis):
530    """
531    Class used for performing multi-modular methods,
532    with the possiblity of removing bad primes.
533    """
534       
535    cdef mod_int last_prime(self):
536        if self.moduli[self.n-1] > self.__last_prime:
537            self.__last_prime = self.moduli[self.n-1]
538        return self.__last_prime
539       
540    def next_prime(self):
541        return self.next_prime_c()
542
543    cdef mod_int next_prime_c(self) except -1:
544        self._extend_moduli(1)
545        return self.moduli[self.n-1]
546       
547    def replace_prime(self, ix):
548        return self.replace_prime_c(ix)
549       
550    cdef mod_int replace_prime_c(self, int ix) except -1:
551        r"""
552        Replace $m_{ix} in the list of moduli with a new
553        prime number greater than all others in the list,
554        and recompute all precomputations.
555       
556        INPUT:
557            ix -- index into list of moduli
558           
559        OUTPUT:
560            p -- the new prime modulus
561        """
562        cdef int i
563        cdef mod_int new_p
564       
565        if ix < 0 or ix >= self.n:
566            raise IndexError, "index out of range"
567
568        new_p = next_prime(self.last_prime())
569        self.__last_prime = new_p
570        self.moduli[ix] = new_p
571
572        self._refresh_products(ix)
573        self._refresh_precomputations(ix)
574        return new_p
Note: See TracBrowser for help on using the repository browser.