Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(276352)

Side by Side Diff: Modules/_decimal/libmpdec/fourstep.c

Issue 7652: Merge C version of decimal into py3k.
Patch Set: Created 7 years, 5 months ago
Left:
Right:
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View unified diff | Download patch
« no previous file with comments | « Modules/_decimal/libmpdec/fnt.h ('k') | Modules/_decimal/libmpdec/fourstep.h » ('j') | no next file with comments »
Toggle Intra-line Diffs ('i') | Expand Comments ('e') | Collapse Comments ('c') | Show Comments Hide Comments ('s')
OLDNEW
(Empty)
1 /*
2 * Copyright (c) 2008-2010 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 <assert.h>
31 #include "numbertheory.h"
32 #include "sixstep.h"
33 #include "transpose.h"
34 #include "umodarith.h"
35 #include "fourstep.h"
36
37
38 /* Bignum: Cache efficient Matrix Fourier Transform for arrays of the
39 form 3 * 2**n (See literature/matrix-transform.txt). */
40
41
42 #ifndef PPRO
43 static inline void
44 std_size3_ntt(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3,
45 mpd_uint_t w3table[3], mpd_uint_t umod)
46 {
47 mpd_uint_t r1, r2;
48 mpd_uint_t w;
49 mpd_uint_t s, tmp;
50
51
52 /* k = 0 -> w = 1 */
53 s = *x1;
54 s = addmod(s, *x2, umod);
55 s = addmod(s, *x3, umod);
56
57 r1 = s;
58
59 /* k = 1 */
60 s = *x1;
61
62 w = w3table[1];
63 tmp = MULMOD(*x2, w);
64 s = addmod(s, tmp, umod);
65
66 w = w3table[2];
67 tmp = MULMOD(*x3, w);
68 s = addmod(s, tmp, umod);
69
70 r2 = s;
71
72 /* k = 2 */
73 s = *x1;
74
75 w = w3table[2];
76 tmp = MULMOD(*x2, w);
77 s = addmod(s, tmp, umod);
78
79 w = w3table[1];
80 tmp = MULMOD(*x3, w);
81 s = addmod(s, tmp, umod);
82
83 *x3 = s;
84 *x2 = r2;
85 *x1 = r1;
86 }
87 #else /* PPRO */
88 static inline void
89 ppro_size3_ntt(mpd_uint_t *x1, mpd_uint_t *x2, mpd_uint_t *x3, mpd_uint_t w3tabl e[3],
90 mpd_uint_t umod, double *dmod, uint32_t dinvmod[3])
91 {
92 mpd_uint_t r1, r2;
93 mpd_uint_t w;
94 mpd_uint_t s, tmp;
95
96
97 /* k = 0 -> w = 1 */
98 s = *x1;
99 s = addmod(s, *x2, umod);
100 s = addmod(s, *x3, umod);
101
102 r1 = s;
103
104 /* k = 1 */
105 s = *x1;
106
107 w = w3table[1];
108 tmp = ppro_mulmod(*x2, w, dmod, dinvmod);
109 s = addmod(s, tmp, umod);
110
111 w = w3table[2];
112 tmp = ppro_mulmod(*x3, w, dmod, dinvmod);
113 s = addmod(s, tmp, umod);
114
115 r2 = s;
116
117 /* k = 2 */
118 s = *x1;
119
120 w = w3table[2];
121 tmp = ppro_mulmod(*x2, w, dmod, dinvmod);
122 s = addmod(s, tmp, umod);
123
124 w = w3table[1];
125 tmp = ppro_mulmod(*x3, w, dmod, dinvmod);
126 s = addmod(s, tmp, umod);
127
128 *x3 = s;
129 *x2 = r2;
130 *x1 = r1;
131 }
132 #endif
133
134
135 /* forward transform, sign = -1; transform length = 3 * 2**n */
136 int
137 four_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
138 {
139 mpd_size_t R = 3; /* number of rows */
140 mpd_size_t C = n / 3; /* number of columns */
141 mpd_uint_t w3table[3];
142 mpd_uint_t kernel, w0, w1, wstep;
143 mpd_uint_t *s, *p0, *p1, *p2;
144 mpd_uint_t umod;
145 #ifdef PPRO
146 double dmod;
147 uint32_t dinvmod[3];
148 #endif
149 mpd_size_t i, k;
150
151
152 assert(n >= 48);
153 assert(n <= 3*MPD_MAXTRANSFORM_2N);
154
155
156 /* Length R transform on the columns. */
157 SETMODULUS(modnum);
158 _mpd_init_w3table(w3table, -1, modnum);
159 for (p0=a, p1=p0+C, p2=p0+2*C; p0<a+C; p0++,p1++,p2++) {
160
161 SIZE3_NTT(p0, p1, p2, w3table);
162 }
163
164 /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
165 kernel = _mpd_getkernel(n, -1, modnum);
166 for (i = 1; i < R; i++) {
167 w0 = 1; /* r**(i*0): initial value for k=0 */
168 w1 = POWMOD(kernel, i); /* r**(i*1): initial value for k=1 */
169 wstep = MULMOD(w1, w1); /* r**(2*i) */
170 for (k = 0; k < C-1; k += 2) {
171 mpd_uint_t x0 = a[i*C+k];
172 mpd_uint_t x1 = a[i*C+k+1];
173 MULMOD2(&x0, w0, &x1, w1);
174 MULMOD2C(&w0, &w1, wstep); /* r**(i*(k+2)) = r**(i*k) * r**(2*i) */
175 a[i*C+k] = x0;
176 a[i*C+k+1] = x1;
177 }
178 }
179
180 /* Length C transform on the rows. */
181 for (s = a; s < a+n; s += C) {
182 if (!six_step_fnt(s, C, modnum)) {
183 return 0;
184 }
185 }
186
187 #if 0 /* An unordered transform is sufficient for convolution. */
188 /* Transpose the matrix. */
189 transpose_3xpow2(a, R, C);
190 #endif
191
192 return 1;
193 }
194
195 /* backward transform, sign = 1; transform length = 3 * 2**n */
196 int
197 inv_four_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
198 {
199 mpd_size_t R = 3; /* number of rows */
200 mpd_size_t C = n / 3; /* number of columns */
201 mpd_uint_t w3table[3];
202 mpd_uint_t kernel, w0, w1, wstep;
203 mpd_uint_t *s, *p0, *p1, *p2;
204 mpd_uint_t umod;
205 #ifdef PPRO
206 double dmod;
207 uint32_t dinvmod[3];
208 #endif
209 mpd_size_t i, k;
210
211
212 assert(n >= 48);
213 assert(n <= 3*MPD_MAXTRANSFORM_2N);
214
215
216 #if 0 /* An unordered transform is sufficient for convolution. */
217 /* Transpose the matrix, producing an R*C matrix. */
218 transpose_3xpow2(a, C, R);
219 #endif
220
221 /* Length C transform on the rows. */
222 for (s = a; s < a+n; s += C) {
223 if (!inv_six_step_fnt(s, C, modnum)) {
224 return 0;
225 }
226 }
227
228 /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
229 SETMODULUS(modnum);
230 kernel = _mpd_getkernel(n, 1, modnum);
231 for (i = 1; i < R; i++) {
232 w0 = 1;
233 w1 = POWMOD(kernel, i);
234 wstep = MULMOD(w1, w1);
235 for (k = 0; k < C; k += 2) {
236 mpd_uint_t x0 = a[i*C+k];
237 mpd_uint_t x1 = a[i*C+k+1];
238 MULMOD2(&x0, w0, &x1, w1);
239 MULMOD2C(&w0, &w1, wstep);
240 a[i*C+k] = x0;
241 a[i*C+k+1] = x1;
242 }
243 }
244
245 /* Length R transform on the columns. */
246 _mpd_init_w3table(w3table, 1, modnum);
247 for (p0=a, p1=p0+C, p2=p0+2*C; p0<a+C; p0++,p1++,p2++) {
248
249 SIZE3_NTT(p0, p1, p2, w3table);
250 }
251
252 return 1;
253 }
254
255
OLDNEW
« no previous file with comments | « Modules/_decimal/libmpdec/fnt.h ('k') | Modules/_decimal/libmpdec/fourstep.h » ('j') | no next file with comments »

RSS Feeds Recent Issues | This issue
This is Rietveld 894c83f36cb7+