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