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 | } |
133 | |