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