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