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

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

Issue 7652: Merge C version of decimal into py3k.
Patch Set: Created 7 years, 8 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/README.txt ('k') | Modules/_decimal/libmpdec/sixstep.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 <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
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 /* An unordered transform is sufficient for convolution. */
119 /* Transpose the matrix. */
120 if (!transpose_pow2(a, R, C)) {
121 return 0;
122 }
123 #endif
124
125 return 1;
126 }
127
128
129 /* reverse transform, sign = 1 */
130 int
131 inv_six_step_fnt(mpd_uint_t *a, mpd_size_t n, int modnum)
132 {
133 struct fnt_params *tparams;
134 mpd_size_t log2n, C, R;
135 mpd_uint_t kernel;
136 mpd_uint_t umod;
137 #ifdef PPRO
138 double dmod;
139 uint32_t dinvmod[3];
140 #endif
141 mpd_uint_t *x, w0, w1, wstep;
142 mpd_size_t i, k;
143
144
145 assert(ispower2(n));
146 assert(n >= 16);
147 assert(n <= MPD_MAXTRANSFORM_2N);
148
149 log2n = mpd_bsr(n);
150 C = ((mpd_size_t)1) << (log2n / 2); /* number of columns */
151 R = ((mpd_size_t)1) << (log2n - (log2n / 2)); /* number of rows */
152
153
154 #if 0 /* An unordered transform is sufficient for convolution. */
155 /* Transpose the matrix, producing an R*C matrix. */
156 if (!transpose_pow2(a, C, R)) {
157 return 0;
158 }
159 #endif
160
161 /* Length C transform on the rows. */
162 if ((tparams = _mpd_init_fnt_params(C, 1, modnum)) == NULL) {
163 return 0;
164 }
165 for (x = a; x < a+n; x += C) {
166 fnt_dif2(x, C, tparams);
167 }
168
169 /* Multiply each matrix element (addressed by i*C+k) by r**(i*k). */
170 SETMODULUS(modnum);
171 kernel = _mpd_getkernel(n, 1, modnum);
172 for (i = 1; i < R; i++) {
173 w0 = 1;
174 w1 = POWMOD(kernel, i);
175 wstep = MULMOD(w1, w1);
176 for (k = 0; k < C; k += 2) {
177 mpd_uint_t x0 = a[i*C+k];
178 mpd_uint_t x1 = a[i*C+k+1];
179 MULMOD2(&x0, w0, &x1, w1);
180 MULMOD2C(&w0, &w1, wstep);
181 a[i*C+k] = x0;
182 a[i*C+k+1] = x1;
183 }
184 }
185
186 /* Transpose the matrix. */
187 if (!transpose_pow2(a, R, C)) {
188 mpd_free(tparams);
189 return 0;
190 }
191
192 /* Length R transform on the rows. */
193 if (R != C) {
194 mpd_free(tparams);
195 if ((tparams = _mpd_init_fnt_params(R, 1, modnum)) == NULL) {
196 return 0;
197 }
198 }
199 for (x = a; x < a+n; x += R) {
200 fnt_dif2(x, R, tparams);
201 }
202 mpd_free(tparams);
203
204 /* Transpose the matrix. */
205 if (!transpose_pow2(a, C, R)) {
206 return 0;
207 }
208
209 return 1;
210 }
211
212
OLDNEW
« no previous file with comments | « Modules/_decimal/libmpdec/README.txt ('k') | Modules/_decimal/libmpdec/sixstep.h » ('j') | no next file with comments »

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