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
_mpd_getkernel(mpd_uint_t n,int sign,int modnum)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 *
_mpd_init_fnt_params(mpd_size_t n,int sign,int modnum)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
_mpd_init_w3table(mpd_uint_t w3table[3],int sign,int modnum)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 }
133