1 #include <NTL/lzz_pX.h>
2
3 NTL_CLIENT
4
5 #define ITER (500)
6
multest()7 void multest()
8 {
9 cerr << "mul";
10 for (long iter = 0; iter < ITER; iter++) {
11 if (iter % 100 == 0) cerr << ".";
12
13 long da = RandomBnd(5000) + 100;
14 long db = RandomBnd(5000) + 100;
15
16 zz_pX a, b, c1, c2;
17
18 random(a, da);
19 random(b, db);
20
21 if (deg(a) < 80 || deg(b) < 80) {
22 cerr << "*";
23 continue;
24 }
25
26 FFTMul(c1, a, b);
27 PlainMul(c2, a, b);
28
29 if (c1 != c2) {
30 cerr << "******* oops\n";
31 break;
32 }
33 }
34
35 cerr << "\n";
36 }
37
38
sqrtest()39 void sqrtest()
40 {
41 cerr << "sqr";
42 for (long iter = 0; iter < ITER; iter++) {
43 if (iter % 100 == 0) cerr << ".";
44
45 long da = RandomBnd(5000) + 100;
46 long db = RandomBnd(5000) + 100;
47
48 zz_pX a, b, c1, c2;
49
50 random(a, da);
51
52 if (deg(a) < 80) {
53 cerr << "*";
54 continue;
55 }
56
57 FFTSqr(c1, a);
58 PlainMul(c2, a, a);
59
60 if (c1 != c2) {
61 cerr << "******* oops\n";
62 break;
63 }
64 }
65
66 cerr << "\n";
67 }
68
69
70
71
mulmodtest()72 void mulmodtest()
73 {
74 cerr << "mulmod";
75 for (long iter = 0; iter < ITER; iter++) {
76 if (iter % 100 == 0) cerr << ".";
77
78 long n = RandomBnd(5000) + 300;
79 long da = RandomBnd(n)+1;
80 long db = RandomBnd(n)+1;
81
82 if (RandomBnd(2)) { da = n; db = n; }
83
84 zz_pX f;
85 random(f, n);
86 SetCoeff(f, n);
87 zz_pXModulus F(f);
88
89 zz_pX a, b, c1, c2;
90 random(a, da);
91 random(b, db);
92
93 MulMod(c1, a, b, F);
94 PlainMul(c2, a, b);
95 rem(c2, c2, f);
96
97 if (c1 != c2) {
98 cerr << "******** oops\n";
99 break;
100 }
101 }
102
103 cerr << "\n";
104 }
105
106
sqrmodtest()107 void sqrmodtest()
108 {
109 cerr << "sqrmod";
110 for (long iter = 0; iter < ITER; iter++) {
111 if (iter % 100 == 0) cerr << ".";
112
113 long n = RandomBnd(5000) + 300;
114 long da = RandomBnd(n)+1;
115 long db = RandomBnd(n)+1;
116
117 if (RandomBnd(2)) { da = n; db = n; }
118
119 zz_pX f;
120 random(f, n);
121 SetCoeff(f, n);
122 zz_pXModulus F(f);
123
124 zz_pX a, b, c1, c2;
125 random(a, da);
126 random(b, db);
127
128 SqrMod(c1, a, F);
129
130 PlainMul(c2, a, a);
131 rem(c2, c2, f);
132
133 if (c1 != c2) {
134 cerr << "******** oops\n";
135 break;
136 }
137 }
138
139 cerr << "\n";
140 }
141
142
143
mulmod1test()144 void mulmod1test()
145 {
146 cerr << "mulmod1";
147 for (long iter = 0; iter < ITER; iter++) {
148 if (iter % 100 == 0) cerr << ".";
149
150 long n = RandomBnd(5000) + 300;
151 long da = RandomBnd(n)+1;
152 long db = RandomBnd(n)+1;
153
154 if (RandomBnd(2)) { da = n; db = n; }
155
156 zz_pX f;
157 random(f, n);
158 SetCoeff(f, n);
159 zz_pXModulus F(f);
160
161 zz_pX a, b, c1, c2;
162 random(a, da);
163 random(b, db);
164
165 zz_pXMultiplier bb;
166 build(bb, b, F);
167
168 MulMod(c1, a, bb, F);
169
170 PlainMul(c2, a, b);
171 rem(c2, c2, f);
172
173 if (c1 != c2) {
174 cerr << "******** oops\n";
175 break;
176 }
177 }
178
179 cerr << "\n";
180 }
181
182
183 namespace NTL {
184
185 void CopyReverse(zz_pX& x, const zz_pX& a, long lo, long hi);
186
187 }
188
189
190
191 struct zz_pXTransMultiplier {
192 zz_pX f0, fbi, b;
193 long shamt, shamt_fbi, shamt_b;
194 };
195
196
197
198
build(zz_pXTransMultiplier & B,const zz_pX & b,const zz_pXModulus & F)199 void build(zz_pXTransMultiplier& B, const zz_pX& b, const zz_pXModulus& F)
200 {
201 long db = deg(b);
202
203 if (db >= F.n) LogicError("build TransMultiplier: bad args");
204
205 zz_pX t;
206
207 LeftShift(t, b, F.n-1);
208 div(t, t, F);
209
210 // we optimize for low degree b
211
212 long d;
213
214 d = deg(t);
215 if (d < 0)
216 B.shamt_fbi = 0;
217 else
218 B.shamt_fbi = F.n-2 - d;
219
220 CopyReverse(B.fbi, t, 0, d);
221
222 // The following code optimizes the case when
223 // f = X^n + low degree poly
224
225 trunc(t, F.f, F.n);
226 d = deg(t);
227 if (d < 0)
228 B.shamt = 0;
229 else
230 B.shamt = d;
231
232 CopyReverse(B.f0, t, 0, d);
233
234 if (db < 0)
235 B.shamt_b = 0;
236 else
237 B.shamt_b = db;
238
239 CopyReverse(B.b, b, 0, db);
240 }
241
242
243
TransMulMod(zz_pX & x,const zz_pX & a,const zz_pXTransMultiplier & B,const zz_pXModulus & F)244 void TransMulMod(zz_pX& x, const zz_pX& a, const zz_pXTransMultiplier& B,
245 const zz_pXModulus& F)
246 {
247 if (deg(a) >= F.n) LogicError("TransMulMod: bad args");
248
249 zz_pX t1, t2;
250
251 mul(t1, a, B.b);
252 RightShift(t1, t1, B.shamt_b);
253
254 mul(t2, a, B.f0);
255 RightShift(t2, t2, B.shamt);
256 trunc(t2, t2, F.n-1);
257
258 mul(t2, t2, B.fbi);
259 if (B.shamt_fbi > 0) LeftShift(t2, t2, B.shamt_fbi);
260 trunc(t2, t2, F.n-1);
261 LeftShift(t2, t2, 1);
262
263 sub(x, t1, t2);
264 }
265
266
267
UpdateMap(vec_zz_p & x,const vec_zz_p & a,const zz_pXTransMultiplier & B,const zz_pXModulus & F)268 void UpdateMap(vec_zz_p& x, const vec_zz_p& a,
269 const zz_pXTransMultiplier& B, const zz_pXModulus& F)
270 {
271 zz_pX xx;
272 TransMulMod(xx, to_zz_pX(a), B, F);
273 x = xx.rep;
274 }
275
276
277
updatetest()278 void updatetest()
279 {
280 cerr << "update";
281 for (long iter = 0; iter < ITER; iter++) {
282 if (iter % 100 == 0) cerr << ".";
283
284 long n = RandomBnd(5000) + 300;
285 long da = RandomBnd(n)+1;
286 long db = RandomBnd(n)+1;
287
288 if (RandomBnd(2)) { da = n; db = n; }
289
290 zz_pX f;
291 random(f, n);
292 SetCoeff(f, n);
293 zz_pXModulus F(f);
294
295 zz_pX a, b;
296 random(a, da);
297 random(b, db);
298
299 zz_pXMultiplier bb1;
300 build(bb1, b, F);
301
302 zz_pXTransMultiplier bb2;
303 build(bb2, b, F);
304
305 Vec<zz_p> x1, x2;
306
307 UpdateMap(x1, a.rep, bb1, F);
308 UpdateMap(x2, a.rep, bb2, F);
309
310
311 if (x1 != x2) {
312 cerr << "******** oops\n";
313 break;
314 }
315 }
316
317 cerr << "\n";
318 }
319
divremtest()320 void divremtest()
321 {
322 cerr << "divrem";
323 for (long iter = 0; iter < ITER; iter++) {
324 if (iter % 100 == 0) cerr << ".";
325
326 long n = RandomBnd(5000) + 300;
327 long dq = RandomBnd(n);
328
329
330 zz_pX f;
331 random(f, n);
332 SetCoeff(f, n);
333 zz_pXModulus F(f);
334
335 zz_pX a, q, r, q1, r1;
336
337 random(a, 2*n-1);
338
339 DivRem(q, r, a, F);
340 rem(r1, a, F);
341 div(q1, a, F);
342
343 if (deg(r) >= n || a != q*f + r || q != q1 || r != r1) {
344 cerr << "******** oops\n";
345 break;
346 }
347 }
348
349 cerr << "\n";
350 }
351
main()352 int main()
353 {
354 long p;
355 p = GenPrime_long(NTL_SP_NBITS);
356
357 zz_p::init(p);
358
359 multest();
360 sqrtest();
361 mulmodtest();
362 sqrmodtest();
363 mulmod1test();
364 divremtest();
365 updatetest();
366
367 zz_p::FFTInit(0);
368
369 cerr << "FFT Prime\n";
370
371 multest();
372 sqrtest();
373 mulmodtest();
374 sqrmodtest();
375 mulmod1test();
376 divremtest();
377 updatetest();
378
379 }
380
381