1
2 #include <NTL/ZZ_pX.h>
3
4 #include <cstdio>
5
6 NTL_CLIENT
7
8
clean_data(double * t)9 double clean_data(double *t)
10 {
11 double x, y, z;
12 long i, ix, iy, n;
13
14 x = t[0]; ix = 0;
15 y = t[0]; iy = 0;
16
17 for (i = 1; i < 5; i++) {
18 if (t[i] < x) {
19 x = t[i];
20 ix = i;
21 }
22 if (t[i] > y) {
23 y = t[i];
24 iy = i;
25 }
26 }
27
28 z = 0; n = 0;
29 for (i = 0; i < 5; i++) {
30 if (i != ix && i != iy) z+= t[i], n++;
31 }
32
33 z = z/n;
34
35 return z;
36 }
37
print_flag()38 void print_flag()
39 {
40
41
42 #if (defined(NTL_CRT_ALTCODE))
43 printf("CRT_ALTCODE ");
44 #else
45 printf("DEFAULT ");
46 #endif
47
48
49 printf("\n");
50
51 }
52
53
main()54 int main()
55 {
56
57 SetSeed(ZZ(0));
58
59 long n, k;
60
61 n = 1024;
62 k = 30*NTL_SP_NBITS;
63
64 ZZ p;
65
66 RandomLen(p, k);
67 if (!IsOdd(p)) p++;
68
69
70 ZZ_p::init(p); // initialization
71
72 ZZ_pX f, g, h, r1, r2, r3;
73
74 random(g, n); // g = random polynomial of degree < n
75 random(h, n); // h = " "
76 random(f, n); // f = " "
77
78 SetCoeff(f, n); // Sets coefficient of X^n to 1
79
80 // For doing arithmetic mod f quickly, one must pre-compute
81 // some information.
82
83 ZZ_pXModulus F;
84 build(F, f);
85
86 PlainMul(r1, g, h); // this uses classical arithmetic
87 PlainRem(r1, r1, f);
88
89 MulMod(r2, g, h, F); // this uses the FFT
90
91 MulMod(r3, g, h, f); // uses FFT, but slower
92
93 // compare the results...
94
95 if (r1 != r2) {
96 printf("999999999999999 ");
97 print_flag();
98 return 0;
99 }
100 else if (r1 != r3) {
101 printf("999999999999999 ");
102 print_flag();
103 return 0;
104 }
105
106 double t;
107 long i;
108 long iter;
109
110 ZZ_pX a, b, c;
111 random(a, n);
112 random(b, n);
113 long da = deg(a);
114 long db = deg(b);
115 long dc = da + db;
116 long l = NextPowerOfTwo(dc+1);
117
118 FFTRep arep, brep, crep;
119 ToFFTRep(arep, a, l, 0, da);
120 ToFFTRep(brep, b, l, 0, db);
121
122 mul(crep, arep, brep);
123
124 ZZ_pXModRep modrep;
125 FromFFTRep(modrep, crep);
126
127 FromZZ_pXModRep(c, modrep, 0, dc);
128
129 iter = 1;
130
131 do {
132 t = GetTime();
133 for (i = 0; i < iter; i++) {
134 FromZZ_pXModRep(c, modrep, 0, dc);
135 }
136 t = GetTime() - t;
137 iter = 2*iter;
138 } while(t < 1);
139
140 iter = iter/2;
141
142 iter = long((3/t)*iter) + 1;
143
144 double tvec[5];
145 long w;
146
147 for (w = 0; w < 5; w++) {
148 t = GetTime();
149 for (i = 0; i < iter; i++) {
150 FromZZ_pXModRep(c, modrep, 0, dc);
151 }
152 t = GetTime() - t;
153 tvec[w] = t;
154 }
155
156
157 t = clean_data(tvec);
158
159 t = floor((t/iter)*1e12);
160
161 if (t < 0 || t >= 1e15)
162 printf("999999999999999 ");
163 else
164 printf("%015.0f ", t);
165
166 printf(" [%ld] ", iter);
167
168 print_flag();
169
170 return 0;
171 }
172