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