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