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