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