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