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 <stdlib.h>
32 #include <assert.h>
33 #include "bits.h"
34 #include "difradix2.h"
35 #include "numbertheory.h"
36 #include "transpose.h"
37 #include "umodarith.h"
38 #include "sixstep.h"
39
40
41 /* Bignum: Cache efficient Matrix Fourier Transform for arrays of the
42 form 2**n (See literature/six-step.txt). */
43
44
45 /* forward transform with sign = -1 */
46 int
six_step_fnt(mpd_uint_t * a,mpd_size_t n,int modnum)47 six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
48 {
49 struct fnt_params *tparams;
50 mpd_size_t log2n, C, R;
51 mpd_uint_t kernel;
52 mpd_uint_t umod;
53 #ifdef PPRO
54 double dmod;
55 uint32_t dinvmod[3];
56 #endif
57 mpd_uint_t *x, w0, w1, wstep;
58 mpd_size_t i, k;
59
60
61 assert(ispower2(n));
62 assert(n >= 16);
63 assert(n <= MPD_MAXTRANSFORM_2N);
64
65 log2n = mpd_bsr(n);
66 C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
67 R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
68
69
70 /* Transpose the matrix. */
71 if (!transpose_pow2(a, R, C)) {
72 return 0;
73 }
74
75 /* Length R transform on the rows. */
76 if ((tparams = _mpd_init_fnt_params(R, -1, modnum)) == NULL) {
77 return 0;
78 }
79 for (x = a; x < a+n; x += R) {
80 fnt_dif2(x, R, tparams);
81 }
82
83 /* Transpose the matrix. */
84 if (!transpose_pow2(a, C, R)) {
85 mpd_free(tparams);
86 return 0;
87 }
88
89 /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
90 SETMODULUS(modnum);
91 kernel = _mpd_getkernel(n, -1, modnum);
92 for (i = 1; i < R; i++) {
93 w0 = 1; /* r**(i*0): initial value for k=0 */
94 w1 = POWMOD(kernel, i); /* r**(i*1): initial value for k=1 */
95 wstep = MULMOD(w1, w1); /* r**(2*i) */
96 for (k = 0; k < C; k += 2) {
97 mpd_uint_t x0 = a[i*C+k];
98 mpd_uint_t x1 = a[i*C+k+1];
99 MULMOD2(&x0, w0, &x1, w1);
100 MULMOD2C(&w0, &w1, wstep); /* r**(i*(k+2)) = r**(i*k) * r**(2*i) */
101 a[i*C+k] = x0;
102 a[i*C+k+1] = x1;
103 }
104 }
105
106 /* Length C transform on the rows. */
107 if (C != R) {
108 mpd_free(tparams);
109 if ((tparams = _mpd_init_fnt_params(C, -1, modnum)) == NULL) {
110 return 0;
111 }
112 }
113 for (x = a; x < a+n; x += C) {
114 fnt_dif2(x, C, tparams);
115 }
116 mpd_free(tparams);
117
118 #if 0
119 /* An unordered transform is sufficient for convolution. */
120 /* Transpose the matrix. */
121 if (!transpose_pow2(a, R, C)) {
122 return 0;
123 }
124 #endif
125
126 return 1;
127 }
128
129
130 /* reverse transform, sign = 1 */
131 int
inv_six_step_fnt(mpd_uint_t * a,mpd_size_t n,int modnum)132 inv_six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
133 {
134 struct fnt_params *tparams;
135 mpd_size_t log2n, C, R;
136 mpd_uint_t kernel;
137 mpd_uint_t umod;
138 #ifdef PPRO
139 double dmod;
140 uint32_t dinvmod[3];
141 #endif
142 mpd_uint_t *x, w0, w1, wstep;
143 mpd_size_t i, k;
144
145
146 assert(ispower2(n));
147 assert(n >= 16);
148 assert(n <= MPD_MAXTRANSFORM_2N);
149
150 log2n = mpd_bsr(n);
151 C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
152 R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
153
154
155 #if 0
156 /* An unordered transform is sufficient for convolution. */
157 /* Transpose the matrix, producing an R*C matrix. */
158 if (!transpose_pow2(a, C, R)) {
159 return 0;
160 }
161 #endif
162
163 /* Length C transform on the rows. */
164 if ((tparams = _mpd_init_fnt_params(C, 1, modnum)) == NULL) {
165 return 0;
166 }
167 for (x = a; x < a+n; x += C) {
168 fnt_dif2(x, C, tparams);
169 }
170
171 /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
172 SETMODULUS(modnum);
173 kernel = _mpd_getkernel(n, 1, modnum);
174 for (i = 1; i < R; i++) {
175 w0 = 1;
176 w1 = POWMOD(kernel, i);
177 wstep = MULMOD(w1, w1);
178 for (k = 0; k < C; k += 2) {
179 mpd_uint_t x0 = a[i*C+k];
180 mpd_uint_t x1 = a[i*C+k+1];
181 MULMOD2(&x0, w0, &x1, w1);
182 MULMOD2C(&w0, &w1, wstep);
183 a[i*C+k] = x0;
184 a[i*C+k+1] = x1;
185 }
186 }
187
188 /* Transpose the matrix. */
189 if (!transpose_pow2(a, R, C)) {
190 mpd_free(tparams);
191 return 0;
192 }
193
194 /* Length R transform on the rows. */
195 if (R != C) {
196 mpd_free(tparams);
197 if ((tparams = _mpd_init_fnt_params(R, 1, modnum)) == NULL) {
198 return 0;
199 }
200 }
201 for (x = a; x < a+n; x += R) {
202 fnt_dif2(x, R, tparams);
203 }
204 mpd_free(tparams);
205
206 /* Transpose the matrix. */
207 if (!transpose_pow2(a, C, R)) {
208 return 0;
209 }
210
211 return 1;
212 }
213
214
215