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