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