Annotation of OpenXM_contrib2/asir2000/fft/polmul.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/fft/polmul.c,v 1.1.1.1 1999/11/10 08:12:27 noro Exp $ */
2: #include "dft.h"
3: extern struct PrimesS Primes[];
4:
5: /*
6: #define TIMING
7: */
8:
9: #ifdef TIMING
10:
11: #define MAXTIMING 20
12:
13: #include <sys/time.h>
14: #include <sys/resource.h>
15: #if 0
16: #if define(hitm)
17: #include <machine/xclock.h>
18: #endif
19: #endif
20:
21: static struct rusage ru_time[MAXTIMING];
22:
23: struct timing {
24: struct timeval user, sys;
25: };
26: static struct timing time_duration[MAXTIMING];
27:
28: #define RECORD_TIME(i) getrusage( RUSAGE_SELF, &ru_time[i] )
29: #endif /* TIMING */
30:
31: void FFT_primes(Num, p_prime, p_primroot, p_d)
32: int Num, *p_prime, *p_primroot, *p_d;
33: {
34: *p_prime = Primes[Num].prime;
35: *p_primroot = Primes[Num].primroot;
36: *p_d = Primes[Num].d;
37: }
38:
39: int FFT_pol_square( d1, C1, Prod, MinMod, wk)
40: unsigned int d1;
41: int MinMod;
42: unsigned int C1[], Prod[], wk[];
43: {
44:
45: unsigned int Low0bits;
46: unsigned int Proot, Prime;
47: double Pinv;
48: void MNpol_square_DFT();
49:
50: if ( MinMod < 0 || MinMod > NPrimes - 1 ) return 2;
51: Prime = Primes[MinMod].prime;
52: Proot = Primes[MinMod].primroot;
53: Low0bits = Primes[MinMod].d;
54: Pinv = (((double)1.0)/((double)Prime));
55:
56: MNpol_square_DFT( d1, C1, Prod, Proot, Low0bits, Prime, Pinv, wk );
57:
58: return 0;
59: }
60:
61: int FFT_pol_product( d1, C1, d2, C2, Prod, MinMod, wk)
62: unsigned int d1, d2;
63: int MinMod;
64: unsigned int C1[], C2[], Prod[], wk[];
65: {
66:
67: unsigned int Low0bits;
68: unsigned int Proot, Prime;
69: double Pinv;
70: void MNpol_product_DFT();
71:
72: if ( MinMod < 0 || MinMod > NPrimes - 1 ) return 2;
73: Prime = Primes[MinMod].prime;
74: Proot = Primes[MinMod].primroot;
75: Low0bits = Primes[MinMod].d;
76: Pinv = (((double)1.0)/((double)Prime));
77: MNpol_product_DFT( d1, C1, d2, C2, Prod, Proot, Low0bits, Prime, Pinv, wk );
78: return 0;
79: }
80:
81: struct oEGT {
82: double exectime,gctime;
83: };
84:
85: extern struct oEGT eg_fore,eg_back;
86:
87: void get_eg(struct oEGT *);
88: void add_eg(struct oEGT *,struct oEGT *, struct oEGT *);
89:
90: void MNpol_product_DFT( d1, C1, d2, C2, Prod, a, low0s, P, Pinv, wk )
91: unsigned int d1, d2, low0s;
92: unsigned int C1[], C2[], Prod[], a, P, wk[];
93: double Pinv;
94: /*
95: * The amount of space of wk[] must be >= (11/2)*2^{\lceil \log_2(d1+d2+1) \rceil}.
96: */
97: {
98: int i, d, n, log2n, halfn;
99: unsigned int *dft1, *dft2, *dftprod, *powa, *true_wk, ninv;
100: struct oEGT eg0,eg1;
101:
102: d = d1 + d2;
103: /**/
104: for ( n = 2, log2n = 1; n <= d; n <<= 1) log2n++;
105: halfn = n >> 1;
106: dft1 = wk, dft2 = &wk[n];
107: dftprod = &dft2[n];
108: powa = &dftprod[n];
109: true_wk = &powa[halfn]; /* 2*n areas */
110: /**/
111: C_PREP_ALPHA( a, low0s, log2n, halfn -1, powa, P, Pinv );
112: #ifdef TIMING
113: RECORD_TIME(1);
114: #endif
115: /**/
116: get_eg(&eg0);
117: C_DFT_FORE( C1, d1 + 1, 1, log2n, powa, dft1, 1, P, Pinv, true_wk );
118: #ifdef TIMING
119: RECORD_TIME(2);
120: #endif
121: C_DFT_FORE( C2, d2 + 1, 1, log2n, powa, dft2, 1, P, Pinv, true_wk );
122: get_eg(&eg1); add_eg(&eg_fore,&eg0,&eg1);
123: #ifdef TIMING
124: RECORD_TIME(3);
125: #endif
126: /**/
127: for ( i = 0; i < n; i++ ) {
128: AxBmodP( dftprod[i], unsigned int, dft1[i], dft2[i], P, Pinv );
129: }
130: #ifdef TIMING
131: RECORD_TIME(4);
132: #endif
133: /**/
134: ninv = P - (unsigned int)(((int)P) >> log2n);
135: get_eg(&eg0);
136: C_DFT_BACK( dftprod, n, 1, log2n, powa, Prod, 1, 0, d + 1, ninv, P, Pinv, true_wk );
137: get_eg(&eg1); add_eg(&eg_back,&eg0,&eg1);
138: }
139:
140: void MNpol_square_DFT( d1, C1, Prod, a, low0s, P, Pinv, wk )
141: unsigned int d1, low0s;
142: unsigned int C1[], Prod[], a, P, wk[];
143: double Pinv;
144: /*
145: * The amount of space of wk[] must be >= (11/2)*2^{\lceil \log_2(d1+d2+1) \rceil}.
146: */
147: {
148: int i, d, n, log2n, halfn;
149: unsigned int *dft1, *dft2, *dftprod, *powa, *true_wk, ninv;
150: struct oEGT eg0,eg1;
151:
152: d = d1 + d1;
153: /**/
154: for ( n = 2, log2n = 1; n <= d; n <<= 1) log2n++;
155: halfn = n >> 1;
156: dft1 = wk, dft2 = &wk[n];
157: dftprod = &dft2[n];
158: powa = &dftprod[n];
159: true_wk = &powa[halfn]; /* 2*n areas */
160: /**/
161: C_PREP_ALPHA( a, low0s, log2n, halfn -1, powa, P, Pinv );
162: #ifdef TIMING
163: RECORD_TIME(1);
164: #endif
165: /**/
166: get_eg(&eg0);
167: C_DFT_FORE( C1, d1 + 1, 1, log2n, powa, dft1, 1, P, Pinv, true_wk );
168: get_eg(&eg1); add_eg(&eg_fore,&eg0,&eg1);
169: /**/
170: for ( i = 0; i < n; i++ ) {
171: AxBmodP( dftprod[i], unsigned int, dft1[i], dft1[i], P, Pinv );
172: }
173: #ifdef TIMING
174: RECORD_TIME(4);
175: #endif
176: /**/
177: ninv = P - (unsigned int)(((int)P) >> log2n);
178: get_eg(&eg0);
179: C_DFT_BACK( dftprod, n, 1, log2n, powa, Prod, 1, 0, d + 1, ninv, P, Pinv, true_wk );
180: get_eg(&eg1); add_eg(&eg_back,&eg0,&eg1);
181: }
182:
183: /****************** TEST CODE ****************/
184:
185: /*
186: #define TEST
187: */
188: #ifdef TEST
189:
190: #define RETI 100
191:
192: #include <stdio.h>
193:
194: #define Prime 65537 /* 2^16 + 1 */
195: #define Low0bits 16
196: #define PrimRoot 3
197:
198: #define PInv (((double)1.0)/((double)Prime))
199:
200: #define MaxNPnts (1<<Low0bits)
201:
202: static unsigned int Pol1[MaxNPnts], Pol2[MaxNPnts], PolProdDFT[MaxNPnts], PolProd[MaxNPnts], wk[6*MaxNPnts];
203:
204: static void generate_pol(), dump_pol(), compare_vals();
205:
206: #ifdef TIMING
207: static void time_spent();
208:
209: #define pr_timing(pTiming) fprintf(stderr, "user:%3d.%06d + sys:%3d.%06d = %3d.%06d", \
210: (pTiming)->user.tv_sec,(pTiming)->user.tv_usec, \
211: (pTiming)->sys.tv_sec,(pTiming)->sys.tv_usec, \
212: (pTiming)->user.tv_sec+(pTiming)->sys.tv_sec, \
213: (pTiming)->user.tv_usec+(pTiming)->sys.tv_usec)
214:
215: #define Pr_timing(i) pr_timing(&time_duration[i])
216: #endif /* TIMING */
217: void main( argc, argv )
218: int argc;
219: char **argv;
220: {
221: int i, ac = argc - 1, d1, d2, d;
222:
223: int ii, *error;
224:
225: char **av = &argv[1];
226: unsigned int mnP = (unsigned int)Prime, mnProot = (unsigned int)PrimRoot;
227: double Pinv = PInv;
228:
229: for ( ; ac >= 2; ac -= 2, av += 2 ) {
230: d = (d1 = atoi( av[0] )) + (d2 = atoi( av[1] ));
231: if ( d1 <= 0 || d2 <= 0 || d >= MaxNPnts ) {
232: fprintf( stderr, "*** invalid degrees %s & %s. (sum < %d)\n",
233: av[0], av[1], MaxNPnts );
234: continue;
235: }
236: generate_pol( d1, Pol1 );
237: /* dump_pol( d1, Pol1 ); */
238: generate_pol( d2, Pol2 );
239: /* dump_pol( d2, Pol2 ); */
240: /**/
241:
242:
243:
244: for ( ii=0; ii < RETI; ii++ ){
245: #ifdef TIMING
246: RECORD_TIME(0);
247: #endif
248: MNpol_product_DFT( d1, Pol1, d2, Pol2, PolProdDFT, mnProot, Low0bits, mnP, Pinv, wk );
249: #ifdef TIMING
250: RECORD_TIME(5);
251: time_spent( 0, 1, 0 );
252: time_spent( 1, 2, 1 );
253: time_spent( 2, 3, 2 );
254: time_spent( 3, 4, 3 );
255: time_spent( 4, 5, 4 );
256: time_spent( 0, 5, 5 );
257: #endif
258: }
259:
260:
261: fprintf( stderr, "DFT mul done ( %d times )\n", RETI );
262:
263:
264:
265:
266:
267: for ( ii=0; ii < RETI; ii++ ){
268: #ifdef TIMING
269: RECORD_TIME(10);
270: #endif
271: MNpol_product( d1, Pol1, d2, Pol2, PolProd, mnP, Pinv );
272: #ifdef TIMING
273: RECORD_TIME(11);
274: time_spent( 10, 11, 10 );
275: #endif
276: }
277:
278:
279:
280: fprintf( stderr, "mul by classical alg. done ( %d times )\n", RETI );
281: /**/
282:
283:
284: /* PolProdDFT[20] = 0.0; */
285:
286:
287: compare_vals( PolProdDFT, PolProd, d+1 , error);
288: if ( *error == 0)
289: fprintf( stderr, "******* Result is OK *******\n" );
290: else
291: fprintf( stderr, "******* Result is NG *******\n" );
292:
293: /**/
294: #ifdef TIMING
295: fprintf( stderr, "DFT mul prep: " ); Pr_timing( 0 );
296: fprintf( stderr, "\nDFT transform 1: " ); Pr_timing( 1 );
297: fprintf( stderr, "\nDFT transform 2: " ); Pr_timing( 2 );
298: fprintf( stderr, "\nDFT mult: " ); Pr_timing( 3 );
299: fprintf( stderr, "\nDFT inv trans: " ); Pr_timing( 4 );
300: fprintf( stderr, "\nDFT pol-mult Tot: " ); Pr_timing( 5 );
301: /**/
302: fprintf( stderr, "\nClassical mult: " ); Pr_timing( 10 );
303: fprintf( stderr, "\n" );
304: #endif /* TIMING */
305: }
306: }
307:
308: static void generate_pol( d, cf )
309: int d;
310: unsigned int cf[];
311: {
312: unsigned int mnP = (unsigned int)Prime;
313: int i, c;
314:
315: for ( i = 0; i < d; i++ ) cf[i] = random() % Prime;
316: while ( (c = random() % Prime) == 0 ) ;
317: cf[d] = c;
318: }
319:
320: static void dump_pol( d, cf )
321: int d;
322: unsigned int cf[];
323: {
324: int i;
325:
326: printf( "Pol of degree %d over Z/(%d) :\n", d, Prime );
327: for ( i = 0; i <= d; i++ ) {
328: if ( (i%5) == 0 ) printf( "\t" );
329: printf( "%10d,", (int)cf[i] );
330: if ( (i%5) == 4 ) printf( "\n" );
331: }
332: if ( (i%5) != 0 ) printf( "\n" );
333: }
334:
335: void compare_vals( v1, v2, n , error)
336: unsigned int v1[], v2[];
337: int n;
338: int *error;
339: {
340: int i, tmp;
341:
342: *error = 0;
343:
344: for ( i = 0; i < n; i ++) {
345: tmp = (int)v1[i] - (int)v2[i];
346: tmp = abs(tmp);
347: if ( tmp > *error ) *error = tmp;
348: }
349:
350:
351:
352: /*
353: int i, j;
354:
355:
356:
357: for ( i = 0; i < n; i += 5 ) {
358: printf( " %6d:%10d", i, (int)v1[i] );
359: for ( j = 1; j < 5 && i + j < n; j++ ) printf( ",%10d", (int)v1[i+j] );
360: printf( "\n" );
361: if ( (int)v1[i] == (int)v2[i] && (int)v1[i+1] == (int)v2[i+1]
362: && (int)v1[i+2] == (int)v2[i+2] && (int)v1[i+3] == (int)v2[i+3]
363: && (int)v1[i+4] == (int)v2[i+4] ) continue;
364: printf( " " );
365: for ( j = 0; j < 5 && i + j < n; j++ )
366: if ( (int)v1[i+j] == (int)v2[i+j] ) printf( " " );
367: else printf( "%10d,", (int)v2[i+j] );
368: printf( "\n" );
369: }
370: */
371:
372:
373:
374: }
375:
376: #ifdef TIMING
377:
378: static void time_spent( s, e, n )
379: int s, e, n;
380: {
381: struct rusage *rus = &ru_time[s], *rue = &ru_time[e];
382: struct timeval *tms, *tme;
383: long sec, usec;
384: struct timing *pt = &time_duration[n];
385:
386: tms = &rus->ru_utime, tme = &rue->ru_utime;
387: sec = tme->tv_sec - tms->tv_sec,
388: usec = tme->tv_usec - tms->tv_usec;
389: if ( usec < 0 ) usec += 1000000, sec--;
390: /* pt->user.tv_sec = sec, pt->user.tv_usec = usec; */
391: pt->user.tv_sec += sec, pt->user.tv_usec += usec;
392:
393:
394:
395: if ( pt->user.tv_usec > 999999 ){
396: pt->user.tv_usec -= 1000000;
397: pt->user.tv_sec++;
398: }
399:
400:
401:
402: /**/
403: tms = &rus->ru_stime, tme = &rue->ru_stime;
404: sec = tme->tv_sec - tms->tv_sec,
405: usec = tme->tv_usec - tms->tv_usec;
406: if ( usec < 0 ) usec += 1000000, sec--;
407: /* pt->sys.tv_sec = sec, pt->sys.tv_usec = usec; */
408: pt->sys.tv_sec += sec, pt->sys.tv_usec += usec;
409:
410:
411:
412: if ( pt->sys.tv_usec > 999999 ){
413: pt->sys.tv_usec -= 1000000;
414: pt->sys.tv_sec++;
415: }
416:
417:
418: }
419: #endif /* TIMING */
420:
421: #endif /* TEST */
422:
423: void MNpol_product( d1, C1, d2, C2, Prod, P, Pinv )
424: unsigned int d1, d2;
425: unsigned int C1[], C2[], Prod[], P;
426: double Pinv;
427: {
428: unsigned int i, j;
429: unsigned int c;
430:
431: c = C1[0];
432: AxBmodP( Prod[0], unsigned int, c, C2[0], P, Pinv );
433: for ( i = 1; i <= d2; i++ ) {
434: AxBmodPnostrchk( Prod[i], unsigned int, c, C2[i], P, Pinv );
435: }
436: c = C1[d1];
437: if ( d1 > d2 ) {
438: for ( i = d2 + 1; i < d1; i++ ) Prod[i] = (unsigned int)0;
439: for ( i = 0; i < d2; i++ ) {
440: AxBmodPnostrchk( Prod[i+d1], unsigned int, c, C2[i], P, Pinv );
441: }
442: } else {
443: j = d2 - d1;
444: for ( i = 0; i <= j; i++ ) {
445: AxBplusCmodPnostrchk( Prod[i+d1], unsigned int, c, C2[i], Prod[i+d1], P, Pinv );
446: }
447: for ( ; i < d2; i++ ) {
448: AxBmodPnostrchk( Prod[i+d1], unsigned int, c, C2[i], P, Pinv );
449: }
450: }
451: AxBmodP( Prod[d1+d2], unsigned int, c, C2[d2], P, Pinv );
452: /**/
453: for ( j = 1; j < d1 - 1; j++ ) {
454: c = C1[j];
455: if ( c == (unsigned int)0 ) {
456: Prod[j] = AstrictmodP( Prod[j], P );
457: continue;
458: }
459: AxBplusCmodP( Prod[j], unsigned int, c, C2[0], Prod[j], P, Pinv );
460: for ( i = 1; i <= d2; i++ ) {
461: AxBplusCmodPnostrchk( Prod[i+j], unsigned int, c, C2[i], Prod[i+j], P, Pinv );
462: }
463: }
464: c = C1[j = d1-1];
465: for ( i = 0; i <= d2; i++ ) {
466: AxBplusCmodP( Prod[i+j], unsigned int, c, C2[i], Prod[i+j], P, Pinv );
467: }
468: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to polmul.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / fft