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