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