Home | History | Annotate | Download | only in libmpdec 
      1 /*
      2  * Copyright (c) 2008-2016 Stefan Krah. All rights reserved.
      3  *
      4  * Redistribution and use in source and binary forms, with or without
      5  * modification, are permitted provided that the following conditions
      6  * are met:
      7  *
      8  * 1. Redistributions of source code must retain the above copyright
      9  *    notice, this list of conditions and the following disclaimer.
     10  *
     11  * 2. Redistributions in binary form must reproduce the above copyright
     12  *    notice, this list of conditions and the following disclaimer in the
     13  *    documentation and/or other materials provided with the distribution.
     14  *
     15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
     16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     17  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     18  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
     19  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     20  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     21  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     22  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     23  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     24  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     25  * SUCH DAMAGE.
     26  */
     27 
     28 
     29 #include "mpdecimal.h"
     30 #include <stdio.h>
     31 #include <assert.h>
     32 #include "numbertheory.h"
     33 #include "umodarith.h"
     34 #include "crt.h"
     35 
     36 
     37 /* Bignum: Chinese Remainder Theorem, extends the maximum transform length. */
     38 
     39 
     40 /* Multiply P1P2 by v, store result in w. */
     41 static inline void
     42 _crt_mulP1P2_3(mpd_uint_t w[3], mpd_uint_t v)
     43 {
     44     mpd_uint_t hi1, hi2, lo;
     45 
     46     _mpd_mul_words(&hi1, &lo, LH_P1P2, v);
     47     w[0] = lo;
     48 
     49     _mpd_mul_words(&hi2, &lo, UH_P1P2, v);
     50     lo = hi1 + lo;
     51     if (lo < hi1) hi2++;
     52 
     53     w[1] = lo;
     54     w[2] = hi2;
     55 }
     56 
     57 /* Add 3 words from v to w. The result is known to fit in w. */
     58 static inline void
     59 _crt_add3(mpd_uint_t w[3], mpd_uint_t v[3])
     60 {
     61     mpd_uint_t carry;
     62     mpd_uint_t s;
     63 
     64     s = w[0] + v[0];
     65     carry = (s < w[0]);
     66     w[0] = s;
     67 
     68     s = w[1] + (v[1] + carry);
     69     carry = (s < w[1]);
     70     w[1] = s;
     71 
     72     w[2] = w[2] + (v[2] + carry);
     73 }
     74 
     75 /* Divide 3 words in u by v, store result in w, return remainder. */
     76 static inline mpd_uint_t
     77 _crt_div3(mpd_uint_t *w, const mpd_uint_t *u, mpd_uint_t v)
     78 {
     79     mpd_uint_t r1 = u[2];
     80     mpd_uint_t r2;
     81 
     82     if (r1 < v) {
     83         w[2] = 0;
     84     }
     85     else {
     86         _mpd_div_word(&w[2], &r1, u[2], v); /* GCOV_NOT_REACHED */
     87     }
     88 
     89     _mpd_div_words(&w[1], &r2, r1, u[1], v);
     90     _mpd_div_words(&w[0], &r1, r2, u[0], v);
     91 
     92     return r1;
     93 }
     94 
     95 
     96 /*
     97  * Chinese Remainder Theorem:
     98  * Algorithm from Joerg Arndt, "Matters Computational",
     99  * Chapter 37.4.1 [http://www.jjj.de/fxt/]
    100  *
    101  * See also Knuth, TAOCP, Volume 2, 4.3.2, exercise 7.
    102  */
    103 
    104 /*
    105  * CRT with carry: x1, x2, x3 contain numbers modulo p1, p2, p3. For each
    106  * triple of members of the arrays, find the unique z modulo p1*p2*p3, with
    107  * zmax = p1*p2*p3 - 1.
    108  *
    109  * In each iteration of the loop, split z into result[i] = z % MPD_RADIX
    110  * and carry = z / MPD_RADIX. Let N be the size of carry[] and cmax the
    111  * maximum carry.
    112  *
    113  * Limits for the 32-bit build:
    114  *
    115  *   N    = 2**96
    116  *   cmax = 7711435591312380274
    117  *
    118  * Limits for the 64 bit build:
    119  *
    120  *   N    = 2**192
    121  *   cmax = 627710135393475385904124401220046371710
    122  *
    123  * The following statements hold for both versions:
    124  *
    125  *   1) cmax + zmax < N, so the addition does not overflow.
    126  *
    127  *   2) (cmax + zmax) / MPD_RADIX == cmax.
    128  *
    129  *   3) If c <= cmax, then c_next = (c + zmax) / MPD_RADIX <= cmax.
    130  */
    131 void
    132 crt3(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3, mpd_size_t rsize)
    133 {
    134     mpd_uint_t p1 = mpd_moduli[P1];
    135     mpd_uint_t umod;
    136 #ifdef PPRO
    137     double dmod;
    138     uint32_t dinvmod[3];
    139 #endif
    140     mpd_uint_t a1, a2, a3;
    141     mpd_uint_t s;
    142     mpd_uint_t z[3], t[3];
    143     mpd_uint_t carry[3] = {0,0,0};
    144     mpd_uint_t hi, lo;
    145     mpd_size_t i;
    146 
    147     for (i = 0; i < rsize; i++) {
    148 
    149         a1 = x1[i];
    150         a2 = x2[i];
    151         a3 = x3[i];
    152 
    153         SETMODULUS(P2);
    154         s = ext_submod(a2, a1, umod);
    155         s = MULMOD(s, INV_P1_MOD_P2);
    156 
    157         _mpd_mul_words(&hi, &lo, s, p1);
    158         lo = lo + a1;
    159         if (lo < a1) hi++;
    160 
    161         SETMODULUS(P3);
    162         s = dw_submod(a3, hi, lo, umod);
    163         s = MULMOD(s, INV_P1P2_MOD_P3);
    164 
    165         z[0] = lo;
    166         z[1] = hi;
    167         z[2] = 0;
    168 
    169         _crt_mulP1P2_3(t, s);
    170         _crt_add3(z, t);
    171         _crt_add3(carry, z);
    172 
    173         x1[i] = _crt_div3(carry, carry, MPD_RADIX);
    174     }
    175 
    176     assert(carry[0] == 0 && carry[1] == 0 && carry[2] == 0);
    177 }
    178 
    179 
    180