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 <stdlib.h>
     31 #include <assert.h>
     32 #include "bits.h"
     33 #include "umodarith.h"
     34 #include "numbertheory.h"
     35 
     36 
     37 /* Bignum: Initialize the Number Theoretic Transform. */
     38 
     39 
     40 /*
     41  * Return the nth root of unity in F(p). This corresponds to e**((2*pi*i)/n)
     42  * in the Fourier transform. We have w**n == 1 (mod p).
     43  *    n := transform length.
     44  *    sign := -1 for forward transform, 1 for backward transform.
     45  *    modnum := one of {P1, P2, P3}.
     46  */
     47 mpd_uint_t
     48 _mpd_getkernel(mpd_uint_t n, int sign, int modnum)
     49 {
     50     mpd_uint_t umod, p, r, xi;
     51 #ifdef PPRO
     52     double dmod;
     53     uint32_t dinvmod[3];
     54 #endif
     55 
     56     SETMODULUS(modnum);
     57     r = mpd_roots[modnum]; /* primitive root of F(p) */
     58     p = umod;
     59     xi = (p-1) / n;
     60 
     61     if (sign == -1)
     62         return POWMOD(r, (p-1-xi));
     63     else
     64         return POWMOD(r, xi);
     65 }
     66 
     67 /*
     68  * Initialize and return transform parameters.
     69  *    n := transform length.
     70  *    sign := -1 for forward transform, 1 for backward transform.
     71  *    modnum := one of {P1, P2, P3}.
     72  */
     73 struct fnt_params *
     74 _mpd_init_fnt_params(mpd_size_t n, int sign, int modnum)
     75 {
     76     struct fnt_params *tparams;
     77     mpd_uint_t umod;
     78 #ifdef PPRO
     79     double dmod;
     80     uint32_t dinvmod[3];
     81 #endif
     82     mpd_uint_t kernel, w;
     83     mpd_uint_t i;
     84     mpd_size_t nhalf;
     85 
     86     assert(ispower2(n));
     87     assert(sign == -1 || sign == 1);
     88     assert(P1 <= modnum && modnum <= P3);
     89 
     90     nhalf = n/2;
     91     tparams = mpd_sh_alloc(sizeof *tparams, nhalf, sizeof (mpd_uint_t));
     92     if (tparams == NULL) {
     93         return NULL;
     94     }
     95 
     96     SETMODULUS(modnum);
     97     kernel = _mpd_getkernel(n, sign, modnum);
     98 
     99     tparams->modnum = modnum;
    100     tparams->modulus = umod;
    101     tparams->kernel = kernel;
    102 
    103     /* wtable[] := w**0, w**1, ..., w**(nhalf-1) */
    104     w = 1;
    105     for (i = 0; i < nhalf; i++) {
    106         tparams->wtable[i] = w;
    107         w = MULMOD(w, kernel);
    108     }
    109 
    110     return tparams;
    111 }
    112 
    113 /* Initialize wtable of size three. */
    114 void
    115 _mpd_init_w3table(mpd_uint_t w3table[3], int sign, int modnum)
    116 {
    117     mpd_uint_t umod;
    118 #ifdef PPRO
    119     double dmod;
    120     uint32_t dinvmod[3];
    121 #endif
    122     mpd_uint_t kernel;
    123 
    124     SETMODULUS(modnum);
    125     kernel = _mpd_getkernel(3, sign, modnum);
    126 
    127     w3table[0] = 1;
    128     w3table[1] = kernel;
    129     w3table[2] = POWMOD(kernel, 2);
    130 }
    131 
    132 
    133