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