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