Annotation of OpenXM_contrib2/asir2000/builtin/poly.c, Revision 1.10
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.10 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/poly.c,v 1.9 2001/05/02 09:03:52 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;
1.10 ! noro 825: case FF_GFS:
! 826: if ( !current_gfs_ext )
! 827: *rp = ONE;
! 828: else
! 829: *rp = DEG(DC(current_gfs_ext));
! 830: break;
1.1 noro 831: default:
832: error("extdeg_ff : current_ff is not set");
833: }
834: }
835:
836: void Pcharacteristic_ff(rp)
837: Q *rp;
838: {
839: N lm;
840:
841: switch ( current_ff ) {
842: case FF_GFP:
843: case FF_GFPN:
844: getmod_lm(&lm); NTOQ(lm,1,*rp); break;
845: case FF_GF2N:
846: STOQ(2,*rp); break;
1.10 ! noro 847: case FF_GFS:
! 848: STOQ(current_gfs_p,*rp); break;
1.1 noro 849: default:
850: error("characteristic_ff : current_ff is not set");
851: }
852: }
853:
854: void Pfield_type_ff(rp)
855: Q *rp;
856: {
857: STOQ(current_ff,*rp);
858: }
859:
860: void Pfield_order_ff(rp)
861: Q *rp;
862: {
863: N n;
864:
865: field_order_ff(&n);
866: NTOQ(n,1,*rp);
867: }
868:
869: void field_order_ff(order)
870: N *order;
871: {
872: UP2 up2;
873: UP up;
874: N m;
875: int d,w;
876:
877: switch ( current_ff ) {
878: case FF_GFP:
879: getmod_lm(order); break;
880: case FF_GF2N:
881: getmod_gf2n(&up2); d = degup2(up2);
882: w = (d>>5)+1;
883: *order = m = NALLOC(w);
884: PL(m)=w;
885: bzero((char *)BD(m),w*sizeof(unsigned int));
886: BD(m)[d>>5] |= 1<<(d&31);
887: break;
888: case FF_GFPN:
889: getmod_lm(&m);
890: getmod_gfpn(&up); pwrn(m,up->d,order); break;
1.10 ! noro 891: case FF_GFS:
! 892: STON(current_gfs_q,*order); break;
1.1 noro 893: default:
894: error("field_order_ff : current_ff is not set");
895: }
896: }
897:
898: void Prandom_ff(rp)
899: Obj *rp;
900: {
901: LM l;
902: GF2N g;
903: GFPN p;
1.7 noro 904: GFS s;
1.1 noro 905:
906: switch ( current_ff ) {
907: case FF_GFP:
908: random_lm(&l); *rp = (Obj)l; break;
909: case FF_GF2N:
910: randomgf2n(&g); *rp = (Obj)g; break;
911: case FF_GFPN:
912: randomgfpn(&p); *rp = (Obj)p; break;
1.7 noro 913: case FF_GFS:
914: randomgfs(&s); *rp = (Obj)s; break;
1.1 noro 915: default:
916: error("random_ff : current_ff is not set");
917: }
918: }
919:
920: void Psimp_ff(arg,rp)
921: NODE arg;
922: Obj *rp;
923: {
924: LM r;
925: GF2N rg;
926: extern lm_lazy;
927:
928: simp_ff((Obj)ARG0(arg),rp);
929: }
930:
931: void simp_ff(p,rp)
932: Obj p;
933: Obj *rp;
934: {
935: Num n;
936: LM r,s;
937: DCP dc,dcr0,dcr;
938: GF2N rg,sg;
939: GFPN rpn,spn;
1.6 noro 940: GFS rs;
1.1 noro 941: P t;
942: Obj obj;
943:
944: lm_lazy = 0;
945: if ( !p )
946: *rp = 0;
947: else if ( OID(p) == O_N ) {
948: switch ( current_ff ) {
949: case FF_GFP:
950: ptolmp((P)p,&t); simplm((LM)t,&s); *rp = (Obj)s;
951: break;
952: case FF_GF2N:
953: ptogf2n((Obj)p,&rg); simpgf2n((GF2N)rg,&sg); *rp = (Obj)sg;
954: break;
955: case FF_GFPN:
956: ntogfpn((Obj)p,&rpn); simpgfpn((GFPN)rpn,&spn); *rp = (Obj)spn;
957: break;
1.6 noro 958: case FF_GFS:
1.8 noro 959: if ( NID((Num)p) == N_GFS )
960: *rp = p;
961: else {
962: ptomp(current_gfs_p,(P)p,&t); mqtogfs(t,&rs);
963: *rp = (Obj)rs;
964: }
1.6 noro 965: break;
1.1 noro 966: default:
967: *rp = (Obj)p;
968: break;
969: }
970: } else if ( OID(p) == O_P ) {
971: for ( dc = DC((P)p), dcr0 = 0; dc; dc = NEXT(dc) ) {
972: simp_ff((Obj)COEF(dc),&obj);
973: if ( obj ) {
974: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)obj;
975: }
976: }
977: if ( !dcr0 )
978: *rp = 0;
979: else {
980: NEXT(dcr) = 0; MKP(VR((P)p),dcr0,t); *rp = (Obj)t;
981: }
982: } else
983: error("simp_ff : not implemented yet");
984: }
985:
986: void getcoef(vl,p,v,d,r)
987: VL vl;
988: P p;
989: V v;
990: Q d;
991: P *r;
992: {
993: P s,t,u,a,b,x;
994: DCP dc;
995: V w;
996:
997: if ( !p )
998: *r = 0;
999: else if ( NUM(p) )
1000: *r = d ? 0 : p;
1001: else if ( (w=VR(p)) == v ) {
1002: for ( dc = DC(p); dc && cmpq(DEG(dc),d); dc = NEXT(dc) );
1003: *r = dc ? COEF(dc) : 0;
1004: } else {
1005: MKV(w,x);
1006: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1007: getcoef(vl,COEF(dc),v,d,&t);
1008: if ( t ) {
1009: pwrp(vl,x,DEG(dc),&u); mulp(vl,t,u,&a);
1010: addp(vl,s,a,&b); s = b;
1011: }
1012: }
1013: *r = s;
1014: }
1015: }
1016:
1017: void Pdeglist(arg,rp)
1018: NODE arg;
1019: LIST *rp;
1020: {
1021: NODE d;
1022:
1023: getdeglist((P)ARG0(arg),VR((P)ARG1(arg)),&d);
1024: MKLIST(*rp,d);
1025: }
1026:
1027: void Pch_mv(arg,rp)
1028: NODE arg;
1029: P *rp;
1030: {
1031: change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
1032: }
1033:
1034: void Pre_mv(arg,rp)
1035: NODE arg;
1036: P *rp;
1037: {
1038: restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
1039: }
1040:
1041: void change_mvar(vl,p,v,r)
1042: VL vl;
1043: P p;
1044: V v;
1045: P *r;
1046: {
1047: Q d;
1048: DCP dc,dc0;
1049: NODE dl;
1050:
1051: if ( !p || NUM(p) || (VR(p) == v) )
1052: *r = p;
1053: else {
1054: getdeglist(p,v,&dl);
1055: for ( dc0 = 0; dl; dl = NEXT(dl) ) {
1056: NEXTDC(dc0,dc); DEG(dc) = d = (Q)BDY(dl);
1057: getcoef(vl,p,v,d,&COEF(dc));
1058: }
1059: NEXT(dc) = 0; MKP(v,dc0,*r);
1060: }
1061: }
1062:
1063: void restore_mvar(vl,p,v,r)
1064: VL vl;
1065: P p;
1066: V v;
1067: P *r;
1068: {
1069: P s,u,a,b,x;
1070: DCP dc;
1071:
1072: if ( !p || NUM(p) || (VR(p) != v) )
1073: *r = p;
1074: else {
1075: MKV(v,x);
1076: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1077: pwrp(vl,x,DEG(dc),&u); mulp(vl,COEF(dc),u,&a);
1078: addp(vl,s,a,&b); s = b;
1079: }
1080: *r = s;
1081: }
1082: }
1083:
1084: void getdeglist(p,v,d)
1085: P p;
1086: V v;
1087: NODE *d;
1088: {
1089: NODE n,n0,d0,d1,d2;
1090: DCP dc;
1091:
1092: if ( !p || NUM(p) ) {
1093: MKNODE(n,0,0); *d = n;
1094: } else if ( VR(p) == v ) {
1095: for ( n0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
1096: NEXTNODE(n0,n); BDY(n) = (pointer)DEG(dc);
1097: }
1098: NEXT(n) = 0; *d = n0;
1099: } else {
1100: for ( dc = DC(p), d0 = 0; dc; dc = NEXT(dc) ) {
1101: getdeglist(COEF(dc),v,&d1); mergedeglist(d0,d1,&d2); d0 = d2;
1102: }
1103: *d = d0;
1104: }
1105: }
1106: void Pmergelist(arg,rp)
1107: NODE arg;
1108: LIST *rp;
1109: {
1110: NODE n;
1111:
1112: asir_assert(ARG0(arg),O_LIST,"mergelist");
1113: asir_assert(ARG1(arg),O_LIST,"mergelist");
1114: mergedeglist(BDY((LIST)ARG0(arg)),BDY((LIST)ARG1(arg)),&n);
1115: MKLIST(*rp,n);
1116: }
1117:
1118: void mergedeglist(d0,d1,dr)
1119: NODE d0,d1,*dr;
1120: {
1121: NODE t0,t,dt;
1122: Q d;
1123: int c;
1124:
1125: if ( !d0 )
1126: *dr = d1;
1127: else {
1128: while ( d1 ) {
1129: dt = d1; d1 = NEXT(d1); d = (Q)BDY(dt);
1130: for ( t0 = 0, t = d0; t; t0 = t, t = NEXT(t) ) {
1131: c = cmpq(d,(Q)BDY(t));
1132: if ( !c )
1133: break;
1134: else if ( c > 0 ) {
1135: if ( !t0 ) {
1136: NEXT(dt) = d0; d0 = dt;
1137: } else {
1138: NEXT(t0) = dt; NEXT(dt) = t;
1139: }
1140: break;
1141: }
1142: }
1143: if ( !t ) {
1144: NEXT(t0) = dt; *dr = d0; return;
1145: }
1146: }
1147: *dr = d0;
1148: }
1149: }
1150:
1151: void Pptomp(arg,rp)
1152: NODE arg;
1153: P *rp;
1154: {
1155: ptomp(QTOS((Q)ARG1(arg)),(P)ARG0(arg),rp);
1156: }
1157:
1158: void Pmptop(arg,rp)
1159: NODE arg;
1160: P *rp;
1161: {
1162: mptop((P)ARG0(arg),rp);
1163: }
1164:
1165: void Pptolmp(arg,rp)
1166: NODE arg;
1167: P *rp;
1168: {
1169: ptolmp((P)ARG0(arg),rp);
1170: }
1171:
1172: void Plmptop(arg,rp)
1173: NODE arg;
1174: P *rp;
1175: {
1176: lmptop((P)ARG0(arg),rp);
1.6 noro 1177: }
1178:
1179: void Psfptop(arg,rp)
1180: NODE arg;
1181: P *rp;
1182: {
1183: sfptop((P)ARG0(arg),rp);
1.1 noro 1184: }
1185:
1186: void Pptogf2n(arg,rp)
1187: NODE arg;
1188: GF2N *rp;
1189: {
1190: ptogf2n((Obj)ARG0(arg),rp);
1191: }
1192:
1193: void Pgf2ntop(arg,rp)
1194: NODE arg;
1195: P *rp;
1196: {
1197: extern V up2_var;
1198:
1199: if ( argc(arg) == 2 )
1200: up2_var = VR((P)ARG1(arg));
1201: gf2ntop((GF2N)ARG0(arg),rp);
1202: }
1203:
1204: void Pgf2ntovect(arg,rp)
1205: NODE arg;
1206: VECT *rp;
1207: {
1208: gf2ntovect((GF2N)ARG0(arg),rp);
1209: }
1210:
1211: void Pptogfpn(arg,rp)
1212: NODE arg;
1213: GF2N *rp;
1214: {
1215: ptogfpn((Obj)ARG0(arg),rp);
1216: }
1217:
1218: void Pgfpntop(arg,rp)
1219: NODE arg;
1220: P *rp;
1221: {
1222: extern V up_var;
1223:
1224: if ( argc(arg) == 2 )
1225: up_var = VR((P)ARG1(arg));
1226: gfpntop((GFPN)ARG0(arg),rp);
1227: }
1228:
1229: void Pureverse(arg,rp)
1230: NODE arg;
1231: P *rp;
1232: {
1233: UP p,r;
1234:
1235: ptoup((P)ARG0(arg),&p);
1.3 noro 1236: if ( argc(arg) == 1 )
1237: reverseup(p,p->d,&r);
1238: else
1239: reverseup(p,QTOS((Q)ARG1(arg)),&r);
1.1 noro 1240: uptop(r,rp);
1241: }
1242:
1243: void Putrunc(arg,rp)
1244: NODE arg;
1245: P *rp;
1246: {
1247: UP p,r;
1248:
1249: ptoup((P)ARG0(arg),&p);
1250: truncup(p,QTOS((Q)ARG1(arg))+1,&r);
1251: uptop(r,rp);
1252: }
1253:
1254: void Pudecomp(arg,rp)
1255: NODE arg;
1256: LIST *rp;
1257: {
1258: P u,l;
1259: UP p,up,low;
1260: NODE n0,n1;
1261:
1262: ptoup((P)ARG0(arg),&p);
1263: decompup(p,QTOS((Q)ARG1(arg))+1,&low,&up);
1264: uptop(low,&l);
1265: uptop(up,&u);
1266: MKNODE(n1,u,0); MKNODE(n0,l,n1);
1267: MKLIST(*rp,n0);
1268: }
1269:
1270: void Purembymul(arg,rp)
1271: NODE arg;
1272: P *rp;
1273: {
1274: UP p1,p2,r;
1275:
1276: if ( !ARG0(arg) || !ARG1(arg) )
1277: *rp = 0;
1278: else {
1279: ptoup((P)ARG0(arg),&p1);
1280: ptoup((P)ARG1(arg),&p2);
1281: rembymulup(p1,p2,&r);
1282: uptop(r,rp);
1283: }
1284: }
1285:
1286: /*
1287: * d1 = deg(p1), d2 = deg(p2)
1288: * d1 <= 2*d2-1
1289: * p2*inv = 1 mod x^d2
1290: */
1291:
1292: void Purembymul_precomp(arg,rp)
1293: NODE arg;
1294: P *rp;
1295: {
1296: UP p1,p2,inv,r;
1297:
1298: if ( !ARG0(arg) || !ARG1(arg) )
1299: *rp = 0;
1300: else {
1301: ptoup((P)ARG0(arg),&p1);
1302: ptoup((P)ARG1(arg),&p2);
1303: ptoup((P)ARG2(arg),&inv);
1304: if ( p1->d >= 2*p2->d ) {
1305: error("urembymul_precomp : degree of 1st arg is too large");
1306: /* fprintf(stderr,"urembymul_precomp : degree of 1st arg is too large"); */
1307: remup(p1,p2,&r);
1308: } else
1309: hybrid_rembymulup_special(current_ff,p1,p2,inv,&r);
1310: uptop(r,rp);
1311: }
1312: }
1313:
1314: void Puinvmod(arg,rp)
1315: NODE arg;
1316: P *rp;
1317: {
1318: UP p,r;
1319:
1320: ptoup((P)ARG0(arg),&p);
1321: invmodup(p,QTOS((Q)ARG1(arg)),&r);
1322: uptop(r,rp);
1323: }
1324:
1325: void Purevinvmod(arg,rp)
1326: NODE arg;
1327: P *rp;
1328: {
1329: UP p,pr,r;
1330:
1331: ptoup((P)ARG0(arg),&p);
1332: reverseup(p,p->d,&pr);
1333: invmodup(pr,QTOS((Q)ARG1(arg)),&r);
1334: uptop(r,rp);
1335: }
1336:
1337: void Ppwrmod_ff(arg,rp)
1338: NODE arg;
1339: P *rp;
1340: {
1341: UP p1,p2;
1342:
1343: ptoup((P)ARG0(arg),&p1);
1344: switch ( current_ff ) {
1345: case FF_GFP:
1346: hybrid_powermodup(p1,&p2); break;
1347: case FF_GF2N:
1348: powermodup_gf2n(p1,&p2); break;
1349: case FF_GFPN:
1.10 ! noro 1350: case FF_GFS:
1.1 noro 1351: powermodup(p1,&p2); break;
1352: default:
1353: error("pwrmod_ff : current_ff is not set");
1354: }
1355: uptop(p2,rp);
1356: }
1357:
1358: void Pgeneric_pwrmod_ff(arg,rp)
1359: NODE arg;
1360: P *rp;
1361: {
1362: UP g,f,r;
1363:
1364: ptoup((P)ARG0(arg),&g);
1365: ptoup((P)ARG1(arg),&f);
1366: switch ( current_ff ) {
1367: case FF_GFP:
1368: hybrid_generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1369: case FF_GF2N:
1370: generic_powermodup_gf2n(g,f,(Q)ARG2(arg),&r); break;
1371: case FF_GFPN:
1.10 ! noro 1372: case FF_GFS:
1.1 noro 1373: generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1374: default:
1375: error("generic_pwrmod_ff : current_ff is not set");
1376: }
1377: uptop(r,rp);
1378: }
1379:
1380: void Ppwrtab_ff(arg,rp)
1381: NODE arg;
1382: VECT *rp;
1383: {
1384: UP f,xp;
1385: UP *tab;
1386: VECT r;
1387: int i,d;
1388:
1389: ptoup((P)ARG0(arg),&f);
1390: ptoup((P)ARG1(arg),&xp);
1391: d = f->d;
1392:
1393: tab = (UP *)ALLOCA(d*sizeof(UP));
1394: switch ( current_ff ) {
1395: case FF_GFP:
1396: hybrid_powertabup(f,xp,tab); break;
1397: case FF_GF2N:
1398: powertabup_gf2n(f,xp,tab); break;
1399: case FF_GFPN:
1.10 ! noro 1400: case FF_GFS:
1.1 noro 1401: powertabup(f,xp,tab); break;
1402: default:
1403: error("pwrtab_ff : current_ff is not set");
1404: }
1405: MKVECT(r,d); *rp = r;
1406: for ( i = 0; i < d; i++ )
1407: uptop(tab[i],(P *)&BDY(r)[i]);
1408: }
1409:
1410: void Pkpwrmod_lm(arg,rp)
1411: NODE arg;
1412: P *rp;
1413: {
1414: UP p1,p2;
1415:
1416: ptoup((P)ARG0(arg),&p1);
1417: powermodup(p1,&p2);
1418: uptop(p2,rp);
1419: }
1420:
1421: void Pkgeneric_pwrmod_lm(arg,rp)
1422: NODE arg;
1423: P *rp;
1424: {
1425: UP g,f,r;
1426:
1427: ptoup((P)ARG0(arg),&g);
1428: ptoup((P)ARG1(arg),&f);
1429: generic_powermodup(g,f,(Q)ARG2(arg),&r);
1430: uptop(r,rp);
1431: }
1432:
1433: void Pkpwrtab_lm(arg,rp)
1434: NODE arg;
1435: VECT *rp;
1436: {
1437: UP f,xp;
1438: UP *tab;
1439: VECT r;
1440: int i,d;
1441:
1442: ptoup((P)ARG0(arg),&f);
1443: ptoup((P)ARG1(arg),&xp);
1444: d = f->d;
1445:
1446: tab = (UP *)ALLOCA(d*sizeof(UP));
1447: powertabup(f,xp,tab);
1448: MKVECT(r,d); *rp = r;
1449: for ( i = 0; i < d; i++ )
1450: uptop(tab[i],(P *)&BDY(r)[i]);
1451: }
1452:
1453: void Plazy_lm(arg,rp)
1454: NODE arg;
1455: Q *rp;
1456: {
1457: lm_lazy = QTOS((Q)ARG0(arg));
1458: *rp = 0;
1459: }
1460:
1461: void Pkmul(arg,rp)
1462: NODE arg;
1463: P *rp;
1464: {
1465: P n1,n2;
1466:
1467: n1 = (P)ARG0(arg); n2 = (P)ARG1(arg);
1468: asir_assert(n1,O_P,"kmul");
1469: asir_assert(n2,O_P,"kmul");
1470: kmulp(CO,n1,n2,rp);
1471: }
1472:
1473: void Pksquare(arg,rp)
1474: NODE arg;
1475: P *rp;
1476: {
1477: P n1;
1478:
1479: n1 = (P)ARG0(arg);
1480: asir_assert(n1,O_P,"ksquare");
1481: ksquarep(CO,n1,rp);
1482: }
1483:
1484: void Pktmul(arg,rp)
1485: NODE arg;
1486: P *rp;
1487: {
1488: UP p1,p2,r;
1489:
1490: ptoup((P)ARG0(arg),&p1);
1491: ptoup((P)ARG1(arg),&p2);
1492: tkmulup(p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1493: uptop(r,rp);
1494: }
1495:
1496: void Pumul(arg,rp)
1497: NODE arg;
1498: P *rp;
1499: {
1500: P a1,a2;
1501: UP p1,p2,r;
1502:
1503: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1504: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1505: mulp(CO,a1,a2,rp);
1506: else {
1507: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1508: error("umul : invalid argument");
1509: ptoup(a1,&p1);
1510: ptoup(a2,&p2);
1511: hybrid_mulup(0,p1,p2,&r);
1512: uptop(r,rp);
1513: }
1514: }
1515:
1516: void Pusquare(arg,rp)
1517: NODE arg;
1518: P *rp;
1519: {
1520: UP p1,p2,r;
1521:
1522: ptoup((P)ARG0(arg),&p1);
1523: hybrid_squareup(0,p1,&r);
1524: uptop(r,rp);
1525: }
1526:
1527: void Putmul(arg,rp)
1528: NODE arg;
1529: P *rp;
1530: {
1531: UP p1,p2,r;
1532:
1533: ptoup((P)ARG0(arg),&p1);
1534: ptoup((P)ARG1(arg),&p2);
1535: hybrid_tmulup(0,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1536: uptop(r,rp);
1537: }
1538:
1539: void Pumul_ff(arg,rp)
1540: NODE arg;
1541: Obj *rp;
1542: {
1543: P a1,a2;
1544: UP p1,p2,r;
1545: P p;
1546:
1547: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1548: ptoup(a1,&p1);
1549: ptoup(a2,&p2);
1550: hybrid_mulup(current_ff,p1,p2,&r);
1551: uptop(r,&p);
1552: simp_ff((Obj)p,rp);
1553: }
1554:
1555: void Pusquare_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: hybrid_squareup(current_ff,p1,&r);
1564: uptop(r,&p);
1565: simp_ff((Obj)p,rp);
1566: }
1567:
1568: void Putmul_ff(arg,rp)
1569: NODE arg;
1570: Obj *rp;
1571: {
1572: UP p1,p2,r;
1573: P p;
1574:
1575: ptoup((P)ARG0(arg),&p1);
1576: ptoup((P)ARG1(arg),&p2);
1577: hybrid_tmulup(current_ff,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1578: uptop(r,&p);
1579: simp_ff((Obj)p,rp);
1580: }
1581:
1582: void Phfmul_lm(arg,rp)
1583: NODE arg;
1584: P *rp;
1585: {
1586: UP p1,p2,r;
1587: UP hi,lo,mid,t,s,p10,p11,p20,p21;
1588: unsigned int m,d;
1589: int i;
1590:
1591: ptoup((P)ARG0(arg),&p1);
1592: ptoup((P)ARG1(arg),&p2);
1593: d = MAX(p1->d,p2->d);
1594: for ( m = 1; m < d; m <<= 1 );
1595: if ( m > d )
1596: m >>= 1;
1597: if ( d-m < 10000 ) {
1598: decompup(p1,m,&p10,&p11); /* p1 = p11*x^m+p10 */
1599: decompup(p2,m,&p20,&p21); /* p2 = p21*x^m+p20 */
1600: fft_mulup_lm(p10,p20,&lo);
1601: kmulup(p11,p21,&hi);
1602: kmulup(p11,p20,&t); kmulup(p10,p21,&s); addup(t,s,&mid);
1603: r = UPALLOC(2*d);
1604: bzero((char *)COEF(r),(2*d+1)*sizeof(Num));
1605: if ( lo )
1606: bcopy((char *)COEF(lo),(char *)COEF(r),(DEG(lo)+1)*sizeof(Num));
1607: if ( hi )
1608: bcopy((char *)COEF(hi),(char *)(COEF(r)+2*m),(DEG(hi)+1)*sizeof(Num));
1609: for ( i = 2*d; i >= 0 && !COEF(r)[i]; i-- );
1610: if ( i < 0 )
1611: r = 0;
1612: else {
1613: DEG(r) = i;
1614: t = UPALLOC(DEG(mid)+m); DEG(t) = DEG(mid)+m;
1615: bzero((char *)COEF(t),(DEG(mid)+m+1)*sizeof(Num));
1616: bcopy((char *)COEF(mid),(char *)(COEF(t)+m),(DEG(mid)+1)*sizeof(Num));
1617: addup(t,r,&s);
1618: r = s;
1619: }
1620: } else
1621: fft_mulup_lm(p1,p2,&r);
1622: uptop(r,rp);
1623: }
1624:
1625: void Pfmultest(arg,rp)
1626: NODE arg;
1627: LIST *rp;
1628: {
1629: P p1,p2,r;
1630: int d1,d2;
1631: UM w1,w2,wr;
1632: unsigned int *f1,*f2,*fr,*w;
1633: int index,mod,root,d,maxint,i;
1634: int cond;
1635: Q prime;
1636: NODE n0,n1;
1637:
1638: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); index = QTOS((Q)ARG2(arg));
1639: FFT_primes(index,&mod,&root,&d);
1640: maxint = 1<<d;
1641: d1 = UDEG(p1); d2 = UDEG(p2);
1642: if ( maxint < d1+d2+1 )
1643: *rp = 0;
1644: else {
1645: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1646: wr = W_UMALLOC(d1+d2);
1647: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1648: f1 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1649: f2 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1650: fr = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1651: w = (unsigned int *)ALLOCA(12*maxint*sizeof(unsigned int));
1652:
1653: for ( i = 0; i <= d1; i++ )
1654: f1[i] = (unsigned int)w1->c[i];
1655: for ( i = 0; i <= d2; i++ )
1656: f2[i] = (unsigned int)w2->c[i];
1657:
1658: cond = FFT_pol_product(d1,f1,d2,f2,fr,index,w);
1659: if ( cond )
1660: error("fmultest : ???");
1661: wr->d = d1+d2;
1662: for ( i = 0; i <= d1+d2; i++ )
1663: wr->c[i] = (unsigned int)fr[i];
1664: umtop(VR(p1),wr,&r);
1665: STOQ(mod,prime);
1666: MKNODE(n1,prime,0);
1667: MKNODE(n0,r,n1);
1668: MKLIST(*rp,n0);
1669: }
1670: }
1671:
1672: void Pkmulum(arg,rp)
1673: NODE arg;
1674: P *rp;
1675: {
1676: P p1,p2;
1677: int d1,d2,mod;
1678: UM w1,w2,wr;
1679:
1680: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
1681: d1 = UDEG(p1); d2 = UDEG(p2);
1682: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1683: wr = W_UMALLOC(d1+d2);
1684: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1685: kmulum(mod,w1,w2,wr);
1686: umtop(VR(p1),wr,rp);
1687: }
1688:
1689: void Pksquareum(arg,rp)
1690: NODE arg;
1691: P *rp;
1692: {
1693: P p1;
1694: int d1,mod;
1695: UM w1,wr;
1696:
1697: p1 = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
1698: d1 = UDEG(p1);
1699: w1 = W_UMALLOC(d1);
1700: wr = W_UMALLOC(2*d1);
1701: ptoum(mod,p1,w1);
1702: kmulum(mod,w1,w1,wr);
1703: umtop(VR(p1),wr,rp);
1704: }
1705:
1706: void Ptracemod_gf2n(arg,rp)
1707: NODE arg;
1708: P *rp;
1709: {
1710: UP g,f,r;
1711:
1712: ptoup((P)ARG0(arg),&g);
1713: ptoup((P)ARG1(arg),&f);
1714: tracemodup_gf2n(g,f,(Q)ARG2(arg),&r);
1715: uptop(r,rp);
1.2 noro 1716: }
1717:
1718: void Pumul_specialmod(arg,rp)
1719: NODE arg;
1720: P *rp;
1721: {
1722: P a1,a2;
1723: UP p1,p2,r;
1724: int i,nmod;
1725: int *modind;
1726: NODE t,n;
1727:
1728: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1729: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1730: mulp(CO,a1,a2,rp);
1731: else {
1732: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1733: error("umul_specialmod : invalid argument");
1734: ptoup(a1,&p1);
1735: ptoup(a2,&p2);
1736: n = BDY((LIST)ARG2(arg));
1737: nmod = length(n);
1738: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1739: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1740: modind[i] = QTOS((Q)BDY(t));
1741: fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
1742: uptop(r,rp);
1743: }
1744: }
1745:
1746: void Pusquare_specialmod(arg,rp)
1747: NODE arg;
1748: P *rp;
1749: {
1750: P a1;
1751: UP p1,r;
1752: int i,nmod;
1753: int *modind;
1754: NODE t,n;
1755:
1756: a1 = (P)ARG0(arg);
1757: if ( !a1 || NUM(a1) )
1758: mulp(CO,a1,a1,rp);
1759: else {
1760: if ( !uzpcheck(a1) )
1761: error("usquare_specialmod : invalid argument");
1762: ptoup(a1,&p1);
1763: n = BDY((LIST)ARG1(arg));
1764: nmod = length(n);
1765: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1766: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1767: modind[i] = QTOS((Q)BDY(t));
1768: fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
1769: uptop(r,rp);
1770: }
1771: }
1772:
1773: void Putmul_specialmod(arg,rp)
1774: NODE arg;
1775: P *rp;
1776: {
1777: P a1,a2;
1778: UP p1,p2,r;
1779: int i,nmod;
1780: int *modind;
1781: NODE t,n;
1782:
1783: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1784: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1785: mulp(CO,a1,a2,rp);
1786: else {
1787: if ( !uzpcheck(a1) || !uzpcheck(a2) || VR(a1) != VR(a2) )
1788: error("utmul_specialmod : invalid argument");
1789: ptoup(a1,&p1);
1790: ptoup(a2,&p2);
1791: n = BDY((LIST)ARG3(arg));
1792: nmod = length(n);
1793: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1794: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1795: modind[i] = QTOS((Q)BDY(t));
1796: fft_mulup_specialmod_main(p1,p2,QTOS((Q)ARG2(arg))+1,modind,nmod,&r);
1797: uptop(r,rp);
1798: }
1.1 noro 1799: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>