![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to poly.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to poly.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.3 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/poly.c,v 1.2 1999/12/27 04:16:30 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4: #include "base.h"
5:
6: void Pranp();
7:
8: void Pumul(),Pumul_ff(),Pusquare(),Pusquare_ff(),Putmul(),Putmul_ff();
9: void Pkmul(),Pksquare(),Pktmul();
10: void Pord(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Psetmod();
11: void Pcoef_gf2n();
12: void getcoef(), getdeglist(), mergedeglist(), change_mvar(), restore_mvar();
13:
1.2 noro 14: void Pp_mag(),Pmaxblen();
1.1 noro 15: void Pmergelist(), Pch_mv(), Pre_mv(), Pdeglist();
16: void Pptomp(),Pmptop();
17: void Pptolmp(),Plmptop();
18: void Pptogf2n(),Pgf2ntop(),Pgf2ntovect();
19: void Pptogfpn(),Pgfpntop();
20: void Pfind_root_gf2n();
1.2 noro 21: void Pumul_specialmod(),Pusquare_specialmod(),Putmul_specialmod();
1.1 noro 22:
23: void Pureverse(),Putrunc(),Pudecomp(),Purembymul(),Purembymul_precomp();
24: void Puinvmod(),Purevinvmod();
25: void Ppwrmod_ff(),Ppwrtab_ff(),Pgeneric_pwrmod_ff();
26: void Pkpwrmod_lm(),Pkpwrtab_lm(),Pkgeneric_pwrmod_lm();
27:
28: void Pkmulum();
29: void Pksquareum();
30:
31: void Pfmultest();
32: void Phfmul_lm();
33: void Plazy_lm();
34:
35: void Psetmod_ff();
36: void Psimp_ff();
37: void Pextdeg_ff();
38: void Pcharacteristic_ff();
39: void Pfield_type_ff();
40: void Pfield_order_ff();
41: void Prandom_ff();
42:
43: void Ptracemod_gf2n();
44: void Psparsemod_gf2n();
45: void Pmultest_gf2n();
46: void Psquaretest_gf2n();
47: void Pinvtest_gf2n();
48: void Pbininv_gf2n();
49: void Prinvtest_gf2n();
50: void Pis_irred_gf2();
51: void Pis_irred_ddd_gf2();
1.2 noro 52: void Pget_next_fft_prime();
1.3 ! noro 53: void Puadj_coef();
1.1 noro 54: void simp_ff(Obj,Obj *);
55: void ranp(int,UP *);
56: void field_order_ff(N *);
57:
58: extern int current_mod;
59: extern GEN_UP2 current_mod_gf2n;
60: extern int lm_lazy;
61:
62: int current_ff;
63:
64: struct ftab poly_tab[] = {
1.3 ! noro 65: {"uadj_coef",Puadj_coef,3},
1.1 noro 66: {"ranp",Pranp,2},
67: {"p_mag",Pp_mag,1},
1.2 noro 68: {"maxblen",Pmaxblen,1},
1.1 noro 69: {"ord",Pord,-1},
70: {"coef0",Pcoef0,-3},
71: {"coef",Pcoef,-3},
72: {"coef_gf2n",Pcoef_gf2n,2},
73: {"deg",Pdeg,2},
74: {"mindeg",Pmindeg,2},
75: {"setmod",Psetmod,-1},
76:
77: {"sparsemod_gf2n",Psparsemod_gf2n,-1},
78:
79: {"setmod_ff",Psetmod_ff,-2},
80: {"simp_ff",Psimp_ff,1},
81: {"extdeg_ff",Pextdeg_ff,0},
82: {"characteristic_ff",Pcharacteristic_ff,0},
83: {"field_type_ff",Pfield_type_ff,0},
84: {"field_order_ff",Pfield_order_ff,0},
85: {"random_ff",Prandom_ff,0},
86:
87: {"deglist",Pdeglist,2},
88: {"mergelist",Pmergelist,2},
89: {"ch_mv",Pch_mv,2},
90: {"re_mv",Pre_mv,2},
91:
92: {"ptomp",Pptomp,2},
93: {"mptop",Pmptop,1},
94:
95: {"ptolmp",Pptolmp,1},
96: {"lmptop",Plmptop,1},
97:
98: {"ptogf2n",Pptogf2n,1},
99: {"gf2ntop",Pgf2ntop,-2},
100: {"gf2ntovect",Pgf2ntovect,1},
101:
102: {"ptogfpn",Pptogfpn,1},
103: {"gfpntop",Pgfpntop,-2},
104:
105: {"kmul",Pkmul,2},
106: {"ksquare",Pksquare,1},
107: {"ktmul",Pktmul,3},
108:
109: {"umul",Pumul,2},
110: {"usquare",Pusquare,1},
111: {"ureverse_inv_as_power_series",Purevinvmod,2},
112: {"uinv_as_power_series",Puinvmod,2},
113:
1.2 noro 114: {"umul_specialmod",Pumul_specialmod,3},
115: {"usquare_specialmod",Pusquare_specialmod,2},
116: {"utmul_specialmod",Putmul_specialmod,4},
117:
1.1 noro 118: {"utmul",Putmul,3},
119: {"umul_ff",Pumul_ff,2},
120: {"usquare_ff",Pusquare_ff,1},
121: {"utmul_ff",Putmul_ff,3},
122:
123: /* for historical reason */
124: {"trunc",Putrunc,2},
125: {"decomp",Pudecomp,2},
126:
127: {"utrunc",Putrunc,2},
128: {"udecomp",Pudecomp,2},
1.3 ! noro 129: {"ureverse",Pureverse,-2},
1.1 noro 130: {"urembymul",Purembymul,2},
131: {"urembymul_precomp",Purembymul_precomp,3},
132:
133: {"lazy_lm",Plazy_lm,1},
134: {"lazy_ff",Plazy_lm,1},
135:
136: {"pwrmod_ff",Ppwrmod_ff,1},
137: {"generic_pwrmod_ff",Pgeneric_pwrmod_ff,3},
138: {"pwrtab_ff",Ppwrtab_ff,2},
139:
140: {"tracemod_gf2n",Ptracemod_gf2n,3},
141: {"b_find_root_gf2n",Pfind_root_gf2n,1},
142:
143: {"is_irred_gf2",Pis_irred_gf2,1},
144: {"is_irred_ddd_gf2",Pis_irred_ddd_gf2,1},
145:
146: {"kpwrmod_lm",Pkpwrmod_lm,1},
147: {"kgeneric_pwrmod_lm",Pkgeneric_pwrmod_lm,3},
148: {"kpwrtab_lm",Pkpwrtab_lm,2},
149:
150: {"kmulum",Pkmulum,3},
151: {"ksquareum",Pksquareum,2},
152:
153: {"fmultest",Pfmultest,3},
154: {"hfmul_lm",Phfmul_lm,2},
155:
156: {"multest_gf2n",Pmultest_gf2n,2},
157: {"squaretest_gf2n",Psquaretest_gf2n,1},
158: {"bininv_gf2n",Pbininv_gf2n,2},
159: {"invtest_gf2n",Pinvtest_gf2n,1},
160: {"rinvtest_gf2n",Prinvtest_gf2n,0},
1.2 noro 161: {"get_next_fft_prime",Pget_next_fft_prime,2},
1.1 noro 162: {0,0,0},
163: };
164:
165: extern V up_var;
1.3 ! noro 166:
! 167: /*
! 168: uadj_coef(F,M,M2)
! 169: if ( F is a non-negative integer )
! 170: return F > M2 ? F-M : M2;
! 171: else
! 172: F = CN*V^N+...+C0
! 173: return uadj_coef(CN,M,M2)*V^N+...+uadj_coef(C0,M,M2);
! 174: */
! 175:
! 176: void Puadj_coef(arg,rp)
! 177: NODE arg;
! 178: P *rp;
! 179: {
! 180: UP f,r;
! 181: N m,m2;
! 182:
! 183: ptoup((P)ARG0(arg),&f);
! 184: m = NM((Q)ARG1(arg));
! 185: m2 = NM((Q)ARG2(arg));
! 186: adj_coefup(f,m,m2,&r);
! 187: uptop(r,rp);
! 188: }
! 189:
1.2 noro 190: /*
191: get_next_fft_prime(StartIndex,Bits)
192: tries to find smallest Index >= StartIndex s.t.
193: 2^(Bits-1)|FFTprime[Index]-1
194: return [Index,Mod] or 0 (not exist)
195: */
196:
197: void Pget_next_fft_prime(arg,rp)
198: NODE arg;
199: LIST *rp;
200: {
201: unsigned int mod,d;
202: int start,bits,i;
203: NODE n;
204: Q q,ind;
205:
206: start = QTOS((Q)ARG0(arg));
207: bits = QTOS((Q)ARG1(arg));
208: for ( i = start; ; i++ ) {
209: get_fft_prime(i,&mod,&d);
210: if ( !mod ) {
211: *rp = 0; return;
212: }
213: if ( bits <= d ) {
214: UTOQ(mod,q);
215: UTOQ(i,ind);
216: n = mknode(2,ind,q);
217: MKLIST(*rp,n);
218: return;
219: }
220: }
221: }
1.1 noro 222:
223: void Pranp(arg,rp)
224: NODE arg;
225: P *rp;
226: {
227: int n;
228: UP c;
229:
230: n = QTOS((Q)ARG0(arg));
231: ranp(n,&c);
232: if ( c ) {
233: up_var = VR((P)ARG1(arg));
234: uptop(c,rp);
235: } else
236: *rp = 0;
237: }
238:
239: void ranp(n,nr)
240: int n;
241: UP *nr;
242: {
243: int i;
244: unsigned int r;
245: Q q;
246: UP c;
247:
248: *nr = c = UPALLOC(n);
249: for ( i = 0; i <= n; i++ ) {
250: r = random();
251: UTOQ(r,q);
252: c->c[i] = (Num)q;
253: }
254: for ( i = n; i >= 0 && !c->c[i]; i-- );
255: if ( i >= 0 )
256: c->d = i;
257: else
258: *nr = 0;
259: }
260:
1.2 noro 261: void Pmaxblen(arg,rp)
262: NODE arg;
263: Q *rp;
264: {
265: int l;
266: l = maxblenp(ARG0(arg));
267: STOQ(l,*rp);
268: }
269:
1.1 noro 270: void Pp_mag(arg,rp)
271: NODE arg;
272: Q *rp;
273: {
274: int l;
275: l = p_mag(ARG0(arg));
276: STOQ(l,*rp);
277: }
278:
279: void Pord(arg,listp)
280: NODE arg;
281: LIST *listp;
282: {
283: NODE n,tn;
284: LIST l;
285: VL vl,tvl,svl;
286: P t;
287: int i,j;
288: V *va;
289: V v;
290:
291: if ( argc(arg) ) {
292: asir_assert(ARG0(arg),O_LIST,"ord");
293: for ( vl = 0, i = 0, n = BDY((LIST)ARG0(arg));
294: n; n = NEXT(n), i++ ) {
295: if ( !vl ) {
296: NEWVL(vl); tvl = vl;
297: } else {
298: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
299: }
300: if ( !(t = (P)BDY(n)) || (OID(t) != O_P) )
301: error("ord : invalid argument");
302: VR(tvl) = VR(t);
303: }
304: va = (V *)ALLOCA(i*sizeof(V));
305: for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
306: va[j] = VR(svl);
307: for ( svl = CO; svl; svl = NEXT(svl) ) {
308: v = VR(svl);
309: for ( j = 0; j < i; j++ )
310: if ( v == va[j] )
311: break;
312: if ( j == i ) {
313: if ( !vl ) {
314: NEWVL(vl); tvl = vl;
315: } else {
316: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
317: }
318: VR(tvl) = v;
319: }
320: }
321: if ( vl )
322: NEXT(tvl) = 0;
323: CO = vl;
324: }
325: for ( n = 0, vl = CO; vl; vl = NEXT(vl) ) {
326: NEXTNODE(n,tn); MKV(VR(vl),t); BDY(tn) = (pointer)t;
327: }
328: NEXT(tn) = 0; MKLIST(l,n); *listp = l;
329: }
330:
331: void Pcoef0(arg,rp)
332: NODE arg;
333: Obj *rp;
334: {
335: Obj t,n;
336: P s;
337: DCP dc;
338: int id;
339: V v;
340: VL vl;
341:
342: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
343: *rp = 0;
344: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
345: *rp = 0;
346: else if ( id == O_N )
347: if ( !n )
348: *rp = t;
349: else
350: *rp = 0;
351: else {
352: if ( argc(arg) == 3 ) {
353: if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
354: reordvar(CO,v,&vl); reorderp(vl,CO,(P)t,&s);
355: } else
356: s = (P)t;
357: if ( VR(s) != v ) {
358: if ( n )
359: *rp = 0;
360: else
361: *rp = t;
362: return;
363: }
364: } else
365: s = (P)t;
366: for ( dc = DC(s); dc && cmpq(DEG(dc),(Q)n); dc = NEXT(dc) );
367: if ( dc )
368: *rp = (Obj)COEF(dc);
369: else
370: *rp = 0;
371: }
372: }
373:
374: void Pcoef(arg,rp)
375: NODE arg;
376: Obj *rp;
377: {
378: Obj t,n;
379: P s;
380: DCP dc;
381: int id;
382: V v;
383:
384: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
385: *rp = 0;
386: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
387: *rp = 0;
388: else if ( id == O_N ) {
389: if ( !n )
390: *rp = t;
391: else
392: *rp = 0;
393: } else {
394: if ( argc(arg) == 3 ) {
395: if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
396: getcoef(CO,(P)t,v,(Q)n,(P *)rp); return;
397: } else
398: s = (P)t;
399: if ( VR(s) != v ) {
400: if ( n )
401: *rp = 0;
402: else
403: *rp = t;
404: return;
405: }
406: } else
407: s = (P)t;
408: for ( dc = DC(s); dc && cmpq(DEG(dc),(Q)n); dc = NEXT(dc) );
409: if ( dc )
410: *rp = (Obj)COEF(dc);
411: else
412: *rp = 0;
413: }
414: }
415:
416: void Pcoef_gf2n(arg,rp)
417: NODE arg;
418: Obj *rp;
419: {
420: Obj t,n;
421: int id,d;
422: UP2 up2;
423:
424: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
425: *rp = 0;
426: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
427: *rp = 0;
428: else if ( id == O_N && NID((Num)t) == N_GF2N ) {
429: d = QTOS((Q)n);
430: up2 = ((GF2N)t)->body;
431: if ( d > degup2(up2) )
432: *rp = 0;
433: else
434: *rp = (Obj)(up2->b[d/BSH]&(((unsigned long)1)<<(d%BSH))?ONE:0);
435: } else
436: *rp = 0;
437: }
438:
439: void Pdeg(arg,rp)
440: NODE arg;
441: Q *rp;
442: {
443: Obj t,v;
444: int d;
445:
446: #if 0
447: if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
448: *rp = 0;
449: else if ( !(v = (Obj)ARG1(arg)) || (VR((P)v) != VR((P)t)) )
450: *rp = 0;
451: else
452: *rp = (Obj)DEG(DC((P)t));
453: #endif
454: if ( !(t = (Obj)ARG0(arg)) )
455: STOQ(-1,*rp);
456: else if ( OID(t) != O_P ) {
457: if ( OID(t) == O_N && NID(t) == N_GF2N
458: && (v=(Obj)ARG1(arg)) && OID(v)== O_N && NID(v) == N_GF2N ) {
459: d = degup2(((GF2N)t)->body);
460: STOQ(d,*rp);
461: } else
462: *rp = 0;
463: } else
464: degp(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
465: }
466:
467: void Pmindeg(arg,rp)
468: NODE arg;
469: Q *rp;
470: {
471: Obj t;
472:
473: if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
474: *rp = 0;
475: else
476: getmindeg(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
477: }
478:
479: void Psetmod(arg,rp)
480: NODE arg;
481: Q *rp;
482: {
483: if ( arg ) {
484: asir_assert(ARG0(arg),O_N,"setmod");
485: current_mod = QTOS((Q)ARG0(arg));
486: }
487: STOQ(current_mod,*rp);
488: }
489:
490: void Psparsemod_gf2n(arg,rp)
491: NODE arg;
492: Q *rp;
493: {
494: int id;
495:
496: if ( arg && current_mod_gf2n )
497: current_mod_gf2n->id = ARG0(arg)?1:0;
498: if ( !current_mod_gf2n )
499: id = -1;
500: else
501: id = current_mod_gf2n->id;
502: STOQ(id,*rp);
503: }
504:
505: void Pmultest_gf2n(arg,rp)
506: NODE arg;
507: GF2N *rp;
508: {
509: GF2N a,b,c;
510: int i;
511:
512: a = (GF2N)ARG0(arg); b = (GF2N)ARG0(arg);
513: for ( i = 0; i < 10000; i++ )
514: mulgf2n(a,b,&c);
515: *rp = c;
516: }
517:
518: void Psquaretest_gf2n(arg,rp)
519: NODE arg;
520: GF2N *rp;
521: {
522: GF2N a,c;
523: int i;
524:
525: a = (GF2N)ARG0(arg);
526: for ( i = 0; i < 10000; i++ )
527: squaregf2n(a,&c);
528: *rp = c;
529: }
530:
531: void Pinvtest_gf2n(arg,rp)
532: NODE arg;
533: GF2N *rp;
534: {
535: GF2N a,c;
536: int i;
537:
538: a = (GF2N)ARG0(arg);
539: for ( i = 0; i < 10000; i++ )
540: invgf2n(a,&c);
541: *rp = c;
542: }
543:
544: void Pbininv_gf2n(arg,rp)
545: NODE arg;
546: GF2N *rp;
547: {
548: UP2 a,inv;
549: int n;
550:
551: a = ((GF2N)ARG0(arg))->body;
552: n = QTOS((Q)ARG1(arg));
553: type1_bin_invup2(a,n,&inv);
554: MKGF2N(inv,*rp);
555: }
556:
557: void Prinvtest_gf2n(rp)
558: Real *rp;
559: {
560: GF2N *a;
561: GF2N c;
562: double t0,t1,r;
563: int i;
564: double get_clock();
565:
566: a = (GF2N *)ALLOCA(1000*sizeof(GF2N));
567: for ( i = 0; i < 1000; i++ ) {
568: randomgf2n(&a[i]);
569: }
570: t0 = get_clock();
571: for ( i = 0; i < 1000; i++ )
572: invgf2n(a[i],&c);
573: t1 = get_clock();
574: r = (t1-t0)/1000;
575: MKReal(r,*rp);
576: }
577:
578: void Pfind_root_gf2n(arg,rp)
579: NODE arg;
580: GF2N *rp;
581: {
582:
583: #if 0
584: UP p;
585:
586: ptoup((P)ARG0(arg),&p);
587: find_root_gf2n(p,rp);
588: #else
589: UP2 p;
590:
591: ptoup2((P)ARG0(arg),&p);
592: find_root_up2(p,rp);
593: #endif
594: }
595:
596: void Pis_irred_gf2(arg,rp)
597: NODE arg;
598: Q *rp;
599: {
600: UP2 t;
601:
602: ptoup2(ARG0(arg),&t);
603: *rp = irredcheckup2(t) ? ONE : 0;
604: }
605:
606: void Pis_irred_ddd_gf2(arg,rp)
607: NODE arg;
608: Q *rp;
609: {
610: UP2 t;
611: int ret;
612:
613: ptoup2(ARG0(arg),&t);
614: ret = irredcheck_dddup2(t);
615: STOQ(ret,*rp);
616: }
617:
618: void Psetmod_ff(arg,rp)
619: NODE arg;
620: Obj *rp;
621: {
622: int ac;
623: Obj mod,defpoly;
624: N n;
625: UP up;
626: UP2 up2;
627: Q q;
628: P p;
629: NODE n0,n1;
630: LIST list;
631:
632: ac = argc(arg);
633: if ( ac == 1 ) {
634: mod = (Obj)ARG0(arg);
635: if ( !mod )
636: error("setmod_ff : invalid argument");
637: switch ( OID(mod) ) {
638: case O_N:
639: current_ff = FF_GFP;
640: setmod_lm(NM((Q)mod)); break;
641: case O_P:
642: current_ff = FF_GF2N;
643: setmod_gf2n((P)mod); break;
644: default:
645: error("setmod_ff : invalid argument");
646: }
647: } else if ( ac == 2 ) {
648: current_ff = FF_GFPN;
649: defpoly = (Obj)ARG0(arg);
650: mod = (Obj)ARG1(arg);
651: if ( !mod || !defpoly )
652: error("setmod_ff : invalid argument");
653: setmod_lm(NM((Q)mod));
654: setmod_gfpn((P)defpoly);
655: }
656: switch ( current_ff ) {
657: case FF_GFP:
658: getmod_lm(&n); NTOQ(n,1,q); *rp = (Obj)q; break;
659: case FF_GF2N:
660: getmod_gf2n(&up2); up2top(up2,&p); *rp = (Obj)p; break;
661: case FF_GFPN:
662: getmod_lm(&n); NTOQ(n,1,q);
663: getmod_gfpn(&up); uptop(up,&p);
664: MKNODE(n1,q,0); MKNODE(n0,p,n1);
665: MKLIST(list,n0);
666: *rp = (Obj)list; break;
667: default:
668: *rp = 0; break;
669: }
670: }
671:
672: void Pextdeg_ff(rp)
673: Q *rp;
674: {
675: int d;
676: UP2 up2;
677: UP up;
678:
679: switch ( current_ff ) {
680: case FF_GFP:
681: *rp = ONE; break;
682: case FF_GF2N:
683: getmod_gf2n(&up2); d = degup2(up2); STOQ(d,*rp); break;
684: case FF_GFPN:
685: getmod_gfpn(&up); STOQ(up->d,*rp); break;
686: default:
687: error("extdeg_ff : current_ff is not set");
688: }
689: }
690:
691: void Pcharacteristic_ff(rp)
692: Q *rp;
693: {
694: N lm;
695:
696: switch ( current_ff ) {
697: case FF_GFP:
698: case FF_GFPN:
699: getmod_lm(&lm); NTOQ(lm,1,*rp); break;
700: case FF_GF2N:
701: STOQ(2,*rp); break;
702: default:
703: error("characteristic_ff : current_ff is not set");
704: }
705: }
706:
707: void Pfield_type_ff(rp)
708: Q *rp;
709: {
710: STOQ(current_ff,*rp);
711: }
712:
713: void Pfield_order_ff(rp)
714: Q *rp;
715: {
716: N n;
717:
718: field_order_ff(&n);
719: NTOQ(n,1,*rp);
720: }
721:
722: void field_order_ff(order)
723: N *order;
724: {
725: UP2 up2;
726: UP up;
727: N m;
728: int d,w;
729:
730: switch ( current_ff ) {
731: case FF_GFP:
732: getmod_lm(order); break;
733: case FF_GF2N:
734: getmod_gf2n(&up2); d = degup2(up2);
735: w = (d>>5)+1;
736: *order = m = NALLOC(w);
737: PL(m)=w;
738: bzero((char *)BD(m),w*sizeof(unsigned int));
739: BD(m)[d>>5] |= 1<<(d&31);
740: break;
741: case FF_GFPN:
742: getmod_lm(&m);
743: getmod_gfpn(&up); pwrn(m,up->d,order); break;
744: default:
745: error("field_order_ff : current_ff is not set");
746: }
747: }
748:
749: void Prandom_ff(rp)
750: Obj *rp;
751: {
752: LM l;
753: GF2N g;
754: GFPN p;
755:
756: switch ( current_ff ) {
757: case FF_GFP:
758: random_lm(&l); *rp = (Obj)l; break;
759: case FF_GF2N:
760: randomgf2n(&g); *rp = (Obj)g; break;
761: case FF_GFPN:
762: randomgfpn(&p); *rp = (Obj)p; break;
763: default:
764: error("random_ff : current_ff is not set");
765: }
766: }
767:
768: void Psimp_ff(arg,rp)
769: NODE arg;
770: Obj *rp;
771: {
772: LM r;
773: GF2N rg;
774: extern lm_lazy;
775:
776: simp_ff((Obj)ARG0(arg),rp);
777: }
778:
779: void simp_ff(p,rp)
780: Obj p;
781: Obj *rp;
782: {
783: Num n;
784: LM r,s;
785: DCP dc,dcr0,dcr;
786: GF2N rg,sg;
787: GFPN rpn,spn;
788: P t;
789: Obj obj;
790:
791: lm_lazy = 0;
792: if ( !p )
793: *rp = 0;
794: else if ( OID(p) == O_N ) {
795: switch ( current_ff ) {
796: case FF_GFP:
797: ptolmp((P)p,&t); simplm((LM)t,&s); *rp = (Obj)s;
798: break;
799: case FF_GF2N:
800: ptogf2n((Obj)p,&rg); simpgf2n((GF2N)rg,&sg); *rp = (Obj)sg;
801: break;
802: case FF_GFPN:
803: ntogfpn((Obj)p,&rpn); simpgfpn((GFPN)rpn,&spn); *rp = (Obj)spn;
804: break;
805: default:
806: *rp = (Obj)p;
807: break;
808: }
809: } else if ( OID(p) == O_P ) {
810: for ( dc = DC((P)p), dcr0 = 0; dc; dc = NEXT(dc) ) {
811: simp_ff((Obj)COEF(dc),&obj);
812: if ( obj ) {
813: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)obj;
814: }
815: }
816: if ( !dcr0 )
817: *rp = 0;
818: else {
819: NEXT(dcr) = 0; MKP(VR((P)p),dcr0,t); *rp = (Obj)t;
820: }
821: } else
822: error("simp_ff : not implemented yet");
823: }
824:
825: void getcoef(vl,p,v,d,r)
826: VL vl;
827: P p;
828: V v;
829: Q d;
830: P *r;
831: {
832: P s,t,u,a,b,x;
833: DCP dc;
834: V w;
835:
836: if ( !p )
837: *r = 0;
838: else if ( NUM(p) )
839: *r = d ? 0 : p;
840: else if ( (w=VR(p)) == v ) {
841: for ( dc = DC(p); dc && cmpq(DEG(dc),d); dc = NEXT(dc) );
842: *r = dc ? COEF(dc) : 0;
843: } else {
844: MKV(w,x);
845: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
846: getcoef(vl,COEF(dc),v,d,&t);
847: if ( t ) {
848: pwrp(vl,x,DEG(dc),&u); mulp(vl,t,u,&a);
849: addp(vl,s,a,&b); s = b;
850: }
851: }
852: *r = s;
853: }
854: }
855:
856: void Pdeglist(arg,rp)
857: NODE arg;
858: LIST *rp;
859: {
860: NODE d;
861:
862: getdeglist((P)ARG0(arg),VR((P)ARG1(arg)),&d);
863: MKLIST(*rp,d);
864: }
865:
866: void Pch_mv(arg,rp)
867: NODE arg;
868: P *rp;
869: {
870: change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
871: }
872:
873: void Pre_mv(arg,rp)
874: NODE arg;
875: P *rp;
876: {
877: restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
878: }
879:
880: void change_mvar(vl,p,v,r)
881: VL vl;
882: P p;
883: V v;
884: P *r;
885: {
886: Q d;
887: DCP dc,dc0;
888: NODE dl;
889:
890: if ( !p || NUM(p) || (VR(p) == v) )
891: *r = p;
892: else {
893: getdeglist(p,v,&dl);
894: for ( dc0 = 0; dl; dl = NEXT(dl) ) {
895: NEXTDC(dc0,dc); DEG(dc) = d = (Q)BDY(dl);
896: getcoef(vl,p,v,d,&COEF(dc));
897: }
898: NEXT(dc) = 0; MKP(v,dc0,*r);
899: }
900: }
901:
902: void restore_mvar(vl,p,v,r)
903: VL vl;
904: P p;
905: V v;
906: P *r;
907: {
908: P s,u,a,b,x;
909: DCP dc;
910:
911: if ( !p || NUM(p) || (VR(p) != v) )
912: *r = p;
913: else {
914: MKV(v,x);
915: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
916: pwrp(vl,x,DEG(dc),&u); mulp(vl,COEF(dc),u,&a);
917: addp(vl,s,a,&b); s = b;
918: }
919: *r = s;
920: }
921: }
922:
923: void getdeglist(p,v,d)
924: P p;
925: V v;
926: NODE *d;
927: {
928: NODE n,n0,d0,d1,d2;
929: DCP dc;
930:
931: if ( !p || NUM(p) ) {
932: MKNODE(n,0,0); *d = n;
933: } else if ( VR(p) == v ) {
934: for ( n0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
935: NEXTNODE(n0,n); BDY(n) = (pointer)DEG(dc);
936: }
937: NEXT(n) = 0; *d = n0;
938: } else {
939: for ( dc = DC(p), d0 = 0; dc; dc = NEXT(dc) ) {
940: getdeglist(COEF(dc),v,&d1); mergedeglist(d0,d1,&d2); d0 = d2;
941: }
942: *d = d0;
943: }
944: }
945: void Pmergelist(arg,rp)
946: NODE arg;
947: LIST *rp;
948: {
949: NODE n;
950:
951: asir_assert(ARG0(arg),O_LIST,"mergelist");
952: asir_assert(ARG1(arg),O_LIST,"mergelist");
953: mergedeglist(BDY((LIST)ARG0(arg)),BDY((LIST)ARG1(arg)),&n);
954: MKLIST(*rp,n);
955: }
956:
957: void mergedeglist(d0,d1,dr)
958: NODE d0,d1,*dr;
959: {
960: NODE t0,t,dt;
961: Q d;
962: int c;
963:
964: if ( !d0 )
965: *dr = d1;
966: else {
967: while ( d1 ) {
968: dt = d1; d1 = NEXT(d1); d = (Q)BDY(dt);
969: for ( t0 = 0, t = d0; t; t0 = t, t = NEXT(t) ) {
970: c = cmpq(d,(Q)BDY(t));
971: if ( !c )
972: break;
973: else if ( c > 0 ) {
974: if ( !t0 ) {
975: NEXT(dt) = d0; d0 = dt;
976: } else {
977: NEXT(t0) = dt; NEXT(dt) = t;
978: }
979: break;
980: }
981: }
982: if ( !t ) {
983: NEXT(t0) = dt; *dr = d0; return;
984: }
985: }
986: *dr = d0;
987: }
988: }
989:
990: void Pptomp(arg,rp)
991: NODE arg;
992: P *rp;
993: {
994: ptomp(QTOS((Q)ARG1(arg)),(P)ARG0(arg),rp);
995: }
996:
997: void Pmptop(arg,rp)
998: NODE arg;
999: P *rp;
1000: {
1001: mptop((P)ARG0(arg),rp);
1002: }
1003:
1004: void Pptolmp(arg,rp)
1005: NODE arg;
1006: P *rp;
1007: {
1008: ptolmp((P)ARG0(arg),rp);
1009: }
1010:
1011: void Plmptop(arg,rp)
1012: NODE arg;
1013: P *rp;
1014: {
1015: lmptop((P)ARG0(arg),rp);
1016: }
1017:
1018: void Pptogf2n(arg,rp)
1019: NODE arg;
1020: GF2N *rp;
1021: {
1022: ptogf2n((Obj)ARG0(arg),rp);
1023: }
1024:
1025: void Pgf2ntop(arg,rp)
1026: NODE arg;
1027: P *rp;
1028: {
1029: extern V up2_var;
1030:
1031: if ( argc(arg) == 2 )
1032: up2_var = VR((P)ARG1(arg));
1033: gf2ntop((GF2N)ARG0(arg),rp);
1034: }
1035:
1036: void Pgf2ntovect(arg,rp)
1037: NODE arg;
1038: VECT *rp;
1039: {
1040: gf2ntovect((GF2N)ARG0(arg),rp);
1041: }
1042:
1043: void Pptogfpn(arg,rp)
1044: NODE arg;
1045: GF2N *rp;
1046: {
1047: ptogfpn((Obj)ARG0(arg),rp);
1048: }
1049:
1050: void Pgfpntop(arg,rp)
1051: NODE arg;
1052: P *rp;
1053: {
1054: extern V up_var;
1055:
1056: if ( argc(arg) == 2 )
1057: up_var = VR((P)ARG1(arg));
1058: gfpntop((GFPN)ARG0(arg),rp);
1059: }
1060:
1061: void Pureverse(arg,rp)
1062: NODE arg;
1063: P *rp;
1064: {
1065: UP p,r;
1066:
1067: ptoup((P)ARG0(arg),&p);
1.3 ! noro 1068: if ( argc(arg) == 1 )
! 1069: reverseup(p,p->d,&r);
! 1070: else
! 1071: reverseup(p,QTOS((Q)ARG1(arg)),&r);
1.1 noro 1072: uptop(r,rp);
1073: }
1074:
1075: void Putrunc(arg,rp)
1076: NODE arg;
1077: P *rp;
1078: {
1079: UP p,r;
1080:
1081: ptoup((P)ARG0(arg),&p);
1082: truncup(p,QTOS((Q)ARG1(arg))+1,&r);
1083: uptop(r,rp);
1084: }
1085:
1086: void Pudecomp(arg,rp)
1087: NODE arg;
1088: LIST *rp;
1089: {
1090: P u,l;
1091: UP p,up,low;
1092: NODE n0,n1;
1093:
1094: ptoup((P)ARG0(arg),&p);
1095: decompup(p,QTOS((Q)ARG1(arg))+1,&low,&up);
1096: uptop(low,&l);
1097: uptop(up,&u);
1098: MKNODE(n1,u,0); MKNODE(n0,l,n1);
1099: MKLIST(*rp,n0);
1100: }
1101:
1102: void Purembymul(arg,rp)
1103: NODE arg;
1104: P *rp;
1105: {
1106: UP p1,p2,r;
1107:
1108: if ( !ARG0(arg) || !ARG1(arg) )
1109: *rp = 0;
1110: else {
1111: ptoup((P)ARG0(arg),&p1);
1112: ptoup((P)ARG1(arg),&p2);
1113: rembymulup(p1,p2,&r);
1114: uptop(r,rp);
1115: }
1116: }
1117:
1118: /*
1119: * d1 = deg(p1), d2 = deg(p2)
1120: * d1 <= 2*d2-1
1121: * p2*inv = 1 mod x^d2
1122: */
1123:
1124: void Purembymul_precomp(arg,rp)
1125: NODE arg;
1126: P *rp;
1127: {
1128: UP p1,p2,inv,r;
1129:
1130: if ( !ARG0(arg) || !ARG1(arg) )
1131: *rp = 0;
1132: else {
1133: ptoup((P)ARG0(arg),&p1);
1134: ptoup((P)ARG1(arg),&p2);
1135: ptoup((P)ARG2(arg),&inv);
1136: if ( p1->d >= 2*p2->d ) {
1137: error("urembymul_precomp : degree of 1st arg is too large");
1138: /* fprintf(stderr,"urembymul_precomp : degree of 1st arg is too large"); */
1139: remup(p1,p2,&r);
1140: } else
1141: hybrid_rembymulup_special(current_ff,p1,p2,inv,&r);
1142: uptop(r,rp);
1143: }
1144: }
1145:
1146: void Puinvmod(arg,rp)
1147: NODE arg;
1148: P *rp;
1149: {
1150: UP p,r;
1151:
1152: ptoup((P)ARG0(arg),&p);
1153: invmodup(p,QTOS((Q)ARG1(arg)),&r);
1154: uptop(r,rp);
1155: }
1156:
1157: void Purevinvmod(arg,rp)
1158: NODE arg;
1159: P *rp;
1160: {
1161: UP p,pr,r;
1162:
1163: ptoup((P)ARG0(arg),&p);
1164: reverseup(p,p->d,&pr);
1165: invmodup(pr,QTOS((Q)ARG1(arg)),&r);
1166: uptop(r,rp);
1167: }
1168:
1169: void Ppwrmod_ff(arg,rp)
1170: NODE arg;
1171: P *rp;
1172: {
1173: UP p1,p2;
1174:
1175: ptoup((P)ARG0(arg),&p1);
1176: switch ( current_ff ) {
1177: case FF_GFP:
1178: hybrid_powermodup(p1,&p2); break;
1179: case FF_GF2N:
1180: powermodup_gf2n(p1,&p2); break;
1181: case FF_GFPN:
1182: powermodup(p1,&p2); break;
1183: default:
1184: error("pwrmod_ff : current_ff is not set");
1185: }
1186: uptop(p2,rp);
1187: }
1188:
1189: void Pgeneric_pwrmod_ff(arg,rp)
1190: NODE arg;
1191: P *rp;
1192: {
1193: UP g,f,r;
1194:
1195: ptoup((P)ARG0(arg),&g);
1196: ptoup((P)ARG1(arg),&f);
1197: switch ( current_ff ) {
1198: case FF_GFP:
1199: hybrid_generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1200: case FF_GF2N:
1201: generic_powermodup_gf2n(g,f,(Q)ARG2(arg),&r); break;
1202: case FF_GFPN:
1203: generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1204: default:
1205: error("generic_pwrmod_ff : current_ff is not set");
1206: }
1207: uptop(r,rp);
1208: }
1209:
1210: void Ppwrtab_ff(arg,rp)
1211: NODE arg;
1212: VECT *rp;
1213: {
1214: UP f,xp;
1215: UP *tab;
1216: VECT r;
1217: int i,d;
1218:
1219: ptoup((P)ARG0(arg),&f);
1220: ptoup((P)ARG1(arg),&xp);
1221: d = f->d;
1222:
1223: tab = (UP *)ALLOCA(d*sizeof(UP));
1224: switch ( current_ff ) {
1225: case FF_GFP:
1226: hybrid_powertabup(f,xp,tab); break;
1227: case FF_GF2N:
1228: powertabup_gf2n(f,xp,tab); break;
1229: case FF_GFPN:
1230: powertabup(f,xp,tab); break;
1231: default:
1232: error("pwrtab_ff : current_ff is not set");
1233: }
1234: MKVECT(r,d); *rp = r;
1235: for ( i = 0; i < d; i++ )
1236: uptop(tab[i],(P *)&BDY(r)[i]);
1237: }
1238:
1239: void Pkpwrmod_lm(arg,rp)
1240: NODE arg;
1241: P *rp;
1242: {
1243: UP p1,p2;
1244:
1245: ptoup((P)ARG0(arg),&p1);
1246: powermodup(p1,&p2);
1247: uptop(p2,rp);
1248: }
1249:
1250: void Pkgeneric_pwrmod_lm(arg,rp)
1251: NODE arg;
1252: P *rp;
1253: {
1254: UP g,f,r;
1255:
1256: ptoup((P)ARG0(arg),&g);
1257: ptoup((P)ARG1(arg),&f);
1258: generic_powermodup(g,f,(Q)ARG2(arg),&r);
1259: uptop(r,rp);
1260: }
1261:
1262: void Pkpwrtab_lm(arg,rp)
1263: NODE arg;
1264: VECT *rp;
1265: {
1266: UP f,xp;
1267: UP *tab;
1268: VECT r;
1269: int i,d;
1270:
1271: ptoup((P)ARG0(arg),&f);
1272: ptoup((P)ARG1(arg),&xp);
1273: d = f->d;
1274:
1275: tab = (UP *)ALLOCA(d*sizeof(UP));
1276: powertabup(f,xp,tab);
1277: MKVECT(r,d); *rp = r;
1278: for ( i = 0; i < d; i++ )
1279: uptop(tab[i],(P *)&BDY(r)[i]);
1280: }
1281:
1282: void Plazy_lm(arg,rp)
1283: NODE arg;
1284: Q *rp;
1285: {
1286: lm_lazy = QTOS((Q)ARG0(arg));
1287: *rp = 0;
1288: }
1289:
1290: void Pkmul(arg,rp)
1291: NODE arg;
1292: P *rp;
1293: {
1294: P n1,n2;
1295:
1296: n1 = (P)ARG0(arg); n2 = (P)ARG1(arg);
1297: asir_assert(n1,O_P,"kmul");
1298: asir_assert(n2,O_P,"kmul");
1299: kmulp(CO,n1,n2,rp);
1300: }
1301:
1302: void Pksquare(arg,rp)
1303: NODE arg;
1304: P *rp;
1305: {
1306: P n1;
1307:
1308: n1 = (P)ARG0(arg);
1309: asir_assert(n1,O_P,"ksquare");
1310: ksquarep(CO,n1,rp);
1311: }
1312:
1313: void Pktmul(arg,rp)
1314: NODE arg;
1315: P *rp;
1316: {
1317: UP p1,p2,r;
1318:
1319: ptoup((P)ARG0(arg),&p1);
1320: ptoup((P)ARG1(arg),&p2);
1321: tkmulup(p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1322: uptop(r,rp);
1323: }
1324:
1325: void Pumul(arg,rp)
1326: NODE arg;
1327: P *rp;
1328: {
1329: P a1,a2;
1330: UP p1,p2,r;
1331:
1332: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1333: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1334: mulp(CO,a1,a2,rp);
1335: else {
1336: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1337: error("umul : invalid argument");
1338: ptoup(a1,&p1);
1339: ptoup(a2,&p2);
1340: hybrid_mulup(0,p1,p2,&r);
1341: uptop(r,rp);
1342: }
1343: }
1344:
1345: void Pusquare(arg,rp)
1346: NODE arg;
1347: P *rp;
1348: {
1349: UP p1,p2,r;
1350:
1351: ptoup((P)ARG0(arg),&p1);
1352: hybrid_squareup(0,p1,&r);
1353: uptop(r,rp);
1354: }
1355:
1356: void Putmul(arg,rp)
1357: NODE arg;
1358: P *rp;
1359: {
1360: UP p1,p2,r;
1361:
1362: ptoup((P)ARG0(arg),&p1);
1363: ptoup((P)ARG1(arg),&p2);
1364: hybrid_tmulup(0,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1365: uptop(r,rp);
1366: }
1367:
1368: void Pumul_ff(arg,rp)
1369: NODE arg;
1370: Obj *rp;
1371: {
1372: P a1,a2;
1373: UP p1,p2,r;
1374: P p;
1375:
1376: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1377: ptoup(a1,&p1);
1378: ptoup(a2,&p2);
1379: hybrid_mulup(current_ff,p1,p2,&r);
1380: uptop(r,&p);
1381: simp_ff((Obj)p,rp);
1382: }
1383:
1384: void Pusquare_ff(arg,rp)
1385: NODE arg;
1386: Obj *rp;
1387: {
1388: UP p1,p2,r;
1389: P p;
1390:
1391: ptoup((P)ARG0(arg),&p1);
1392: hybrid_squareup(current_ff,p1,&r);
1393: uptop(r,&p);
1394: simp_ff((Obj)p,rp);
1395: }
1396:
1397: void Putmul_ff(arg,rp)
1398: NODE arg;
1399: Obj *rp;
1400: {
1401: UP p1,p2,r;
1402: P p;
1403:
1404: ptoup((P)ARG0(arg),&p1);
1405: ptoup((P)ARG1(arg),&p2);
1406: hybrid_tmulup(current_ff,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1407: uptop(r,&p);
1408: simp_ff((Obj)p,rp);
1409: }
1410:
1411: void Phfmul_lm(arg,rp)
1412: NODE arg;
1413: P *rp;
1414: {
1415: UP p1,p2,r;
1416: UP hi,lo,mid,t,s,p10,p11,p20,p21;
1417: unsigned int m,d;
1418: int i;
1419:
1420: ptoup((P)ARG0(arg),&p1);
1421: ptoup((P)ARG1(arg),&p2);
1422: d = MAX(p1->d,p2->d);
1423: for ( m = 1; m < d; m <<= 1 );
1424: if ( m > d )
1425: m >>= 1;
1426: if ( d-m < 10000 ) {
1427: decompup(p1,m,&p10,&p11); /* p1 = p11*x^m+p10 */
1428: decompup(p2,m,&p20,&p21); /* p2 = p21*x^m+p20 */
1429: fft_mulup_lm(p10,p20,&lo);
1430: kmulup(p11,p21,&hi);
1431: kmulup(p11,p20,&t); kmulup(p10,p21,&s); addup(t,s,&mid);
1432: r = UPALLOC(2*d);
1433: bzero((char *)COEF(r),(2*d+1)*sizeof(Num));
1434: if ( lo )
1435: bcopy((char *)COEF(lo),(char *)COEF(r),(DEG(lo)+1)*sizeof(Num));
1436: if ( hi )
1437: bcopy((char *)COEF(hi),(char *)(COEF(r)+2*m),(DEG(hi)+1)*sizeof(Num));
1438: for ( i = 2*d; i >= 0 && !COEF(r)[i]; i-- );
1439: if ( i < 0 )
1440: r = 0;
1441: else {
1442: DEG(r) = i;
1443: t = UPALLOC(DEG(mid)+m); DEG(t) = DEG(mid)+m;
1444: bzero((char *)COEF(t),(DEG(mid)+m+1)*sizeof(Num));
1445: bcopy((char *)COEF(mid),(char *)(COEF(t)+m),(DEG(mid)+1)*sizeof(Num));
1446: addup(t,r,&s);
1447: r = s;
1448: }
1449: } else
1450: fft_mulup_lm(p1,p2,&r);
1451: uptop(r,rp);
1452: }
1453:
1454: void Pfmultest(arg,rp)
1455: NODE arg;
1456: LIST *rp;
1457: {
1458: P p1,p2,r;
1459: int d1,d2;
1460: UM w1,w2,wr;
1461: unsigned int *f1,*f2,*fr,*w;
1462: int index,mod,root,d,maxint,i;
1463: int cond;
1464: Q prime;
1465: NODE n0,n1;
1466:
1467: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); index = QTOS((Q)ARG2(arg));
1468: FFT_primes(index,&mod,&root,&d);
1469: maxint = 1<<d;
1470: d1 = UDEG(p1); d2 = UDEG(p2);
1471: if ( maxint < d1+d2+1 )
1472: *rp = 0;
1473: else {
1474: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1475: wr = W_UMALLOC(d1+d2);
1476: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1477: f1 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1478: f2 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1479: fr = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1480: w = (unsigned int *)ALLOCA(12*maxint*sizeof(unsigned int));
1481:
1482: for ( i = 0; i <= d1; i++ )
1483: f1[i] = (unsigned int)w1->c[i];
1484: for ( i = 0; i <= d2; i++ )
1485: f2[i] = (unsigned int)w2->c[i];
1486:
1487: cond = FFT_pol_product(d1,f1,d2,f2,fr,index,w);
1488: if ( cond )
1489: error("fmultest : ???");
1490: wr->d = d1+d2;
1491: for ( i = 0; i <= d1+d2; i++ )
1492: wr->c[i] = (unsigned int)fr[i];
1493: umtop(VR(p1),wr,&r);
1494: STOQ(mod,prime);
1495: MKNODE(n1,prime,0);
1496: MKNODE(n0,r,n1);
1497: MKLIST(*rp,n0);
1498: }
1499: }
1500:
1501: void Pkmulum(arg,rp)
1502: NODE arg;
1503: P *rp;
1504: {
1505: P p1,p2;
1506: int d1,d2,mod;
1507: UM w1,w2,wr;
1508:
1509: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
1510: d1 = UDEG(p1); d2 = UDEG(p2);
1511: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1512: wr = W_UMALLOC(d1+d2);
1513: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1514: kmulum(mod,w1,w2,wr);
1515: umtop(VR(p1),wr,rp);
1516: }
1517:
1518: void Pksquareum(arg,rp)
1519: NODE arg;
1520: P *rp;
1521: {
1522: P p1;
1523: int d1,mod;
1524: UM w1,wr;
1525:
1526: p1 = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
1527: d1 = UDEG(p1);
1528: w1 = W_UMALLOC(d1);
1529: wr = W_UMALLOC(2*d1);
1530: ptoum(mod,p1,w1);
1531: kmulum(mod,w1,w1,wr);
1532: umtop(VR(p1),wr,rp);
1533: }
1534:
1535: void Ptracemod_gf2n(arg,rp)
1536: NODE arg;
1537: P *rp;
1538: {
1539: UP g,f,r;
1540:
1541: ptoup((P)ARG0(arg),&g);
1542: ptoup((P)ARG1(arg),&f);
1543: tracemodup_gf2n(g,f,(Q)ARG2(arg),&r);
1544: uptop(r,rp);
1.2 noro 1545: }
1546:
1547: void Pumul_specialmod(arg,rp)
1548: NODE arg;
1549: P *rp;
1550: {
1551: P a1,a2;
1552: UP p1,p2,r;
1553: int i,nmod;
1554: int *modind;
1555: NODE t,n;
1556:
1557: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1558: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1559: mulp(CO,a1,a2,rp);
1560: else {
1561: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1562: error("umul_specialmod : invalid argument");
1563: ptoup(a1,&p1);
1564: ptoup(a2,&p2);
1565: n = BDY((LIST)ARG2(arg));
1566: nmod = length(n);
1567: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1568: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1569: modind[i] = QTOS((Q)BDY(t));
1570: fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
1571: uptop(r,rp);
1572: }
1573: }
1574:
1575: void Pusquare_specialmod(arg,rp)
1576: NODE arg;
1577: P *rp;
1578: {
1579: P a1;
1580: UP p1,r;
1581: int i,nmod;
1582: int *modind;
1583: NODE t,n;
1584:
1585: a1 = (P)ARG0(arg);
1586: if ( !a1 || NUM(a1) )
1587: mulp(CO,a1,a1,rp);
1588: else {
1589: if ( !uzpcheck(a1) )
1590: error("usquare_specialmod : invalid argument");
1591: ptoup(a1,&p1);
1592: n = BDY((LIST)ARG1(arg));
1593: nmod = length(n);
1594: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1595: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1596: modind[i] = QTOS((Q)BDY(t));
1597: fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
1598: uptop(r,rp);
1599: }
1600: }
1601:
1602: void Putmul_specialmod(arg,rp)
1603: NODE arg;
1604: P *rp;
1605: {
1606: P a1,a2;
1607: UP p1,p2,r;
1608: int i,nmod;
1609: int *modind;
1610: NODE t,n;
1611:
1612: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1613: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1614: mulp(CO,a1,a2,rp);
1615: else {
1616: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1617: error("utmul_specialmod : invalid argument");
1618: ptoup(a1,&p1);
1619: ptoup(a2,&p2);
1620: n = BDY((LIST)ARG3(arg));
1621: nmod = length(n);
1622: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1623: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1624: modind[i] = QTOS((Q)BDY(t));
1625: fft_mulup_specialmod_main(p1,p2,QTOS((Q)ARG2(arg))+1,modind,nmod,&r);
1626: uptop(r,rp);
1627: }
1.1 noro 1628: }