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