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