![[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.24 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/poly.c,v 1.23 2011/07/20 03:19:11 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: {
1.23 noro 448: NODE n,tn,p,opt;
449: char *key;
450: Obj value;
451: int overwrite=0;
1.1 noro 452: LIST l;
453: VL vl,tvl,svl;
454: P t;
455: int i,j;
456: V *va;
457: V v;
458:
1.23 noro 459: if ( current_option ) {
460: for ( opt = current_option; opt; opt = NEXT(opt) ) {
461: p = BDY((LIST)BDY(opt));
462: key = BDY((STRING)BDY(p));
463: value = (Obj)BDY(NEXT(p));
464: if ( !strcmp(key,"overwrite") && value ) {
465: overwrite = value ? 1 : 0;
466: break;
467: }
468: }
469: }
470:
1.1 noro 471: if ( argc(arg) ) {
472: asir_assert(ARG0(arg),O_LIST,"ord");
473: for ( vl = 0, i = 0, n = BDY((LIST)ARG0(arg));
474: n; n = NEXT(n), i++ ) {
475: if ( !vl ) {
476: NEWVL(vl); tvl = vl;
477: } else {
478: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
479: }
480: if ( !(t = (P)BDY(n)) || (OID(t) != O_P) )
481: error("ord : invalid argument");
482: VR(tvl) = VR(t);
483: }
1.23 noro 484: if ( !overwrite ) {
485: va = (V *)ALLOCA(i*sizeof(V));
486: for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
487: va[j] = VR(svl);
488: for ( svl = CO; svl; svl = NEXT(svl) ) {
489: v = VR(svl);
490: for ( j = 0; j < i; j++ )
491: if ( v == va[j] )
492: break;
493: if ( j == i ) {
494: if ( !vl ) {
495: NEWVL(vl); tvl = vl;
496: } else {
497: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
498: }
499: VR(tvl) = v;
1.1 noro 500: }
501: }
1.24 ! noro 502: } else {
! 503: for ( svl = vl; svl; svl = NEXT(svl) ) {
! 504: if ( svl->v->attr == (pointer)V_PF )
! 505: ((PFINS)svl->v->priv)->pf->ins = 0;
! 506: }
1.1 noro 507: }
508: if ( vl )
509: NEXT(tvl) = 0;
510: CO = vl;
511: }
512: for ( n = 0, vl = CO; vl; vl = NEXT(vl) ) {
513: NEXTNODE(n,tn); MKV(VR(vl),t); BDY(tn) = (pointer)t;
514: }
515: NEXT(tn) = 0; MKLIST(l,n); *listp = l;
516: }
517:
1.22 noro 518: void Premove_vars(NODE arg,LIST *listp)
519: {
520: NODE l,nd,tnd;
521: V *v,*va;
522: int n,na,i,j;
523: VL vl,vl1;
524: P t;
525: LIST list;
526:
527: asir_assert(ARG0(arg),O_LIST,"remove_vars");
528: l = BDY((LIST)ARG0(arg)); n = length(l);
529: v = (V *)ALLOCA(n*sizeof(V));
530: for ( i = 0; i < n; i++, l = NEXT(l) )
531: if ( !(t = (P)BDY(l)) || (OID(t) != O_P) )
532: error("ord : invalid argument");
533: else v[i] = VR(t);
534:
535: for ( na = 0, vl = CO; vl; vl = NEXT(vl), na++ );
536: va = (V *)ALLOCA(na*sizeof(V));
537: for ( i = 0, vl = CO; i < na; i++, vl = NEXT(vl) ) va[i] = VR(vl);
538: for ( i = 0; i < na; i++ )
539: for ( j = 0; j < n; j++ ) if ( va[i] == v[j] ) va[i] = 0;
540: for ( vl = 0, i = na-1; i >= 0; i-- )
541: if ( va[i] ) {
542: NEWVL(vl1); VR(vl1) = va[i]; NEXT(vl1) = vl; vl = vl1;
543: }
544: CO = vl;
545: for ( nd = 0, vl = CO; vl; vl = NEXT(vl) ) {
546: NEXTNODE(nd,tnd); MKV(VR(vl),t); BDY(tnd) = (pointer)t;
547: }
548: if ( nd ) NEXT(tnd) = 0;
549: MKLIST(list,nd); *listp = list;
550: }
551:
1.16 noro 552: void Pcoef0(NODE arg,Obj *rp)
1.1 noro 553: {
554: Obj t,n;
555: P s;
556: DCP dc;
557: int id;
558: V v;
559: VL vl;
560:
561: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
562: *rp = 0;
563: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
564: *rp = 0;
565: else if ( id == O_N )
566: if ( !n )
567: *rp = t;
568: else
569: *rp = 0;
570: else {
571: if ( argc(arg) == 3 ) {
572: if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
573: reordvar(CO,v,&vl); reorderp(vl,CO,(P)t,&s);
574: } else
575: s = (P)t;
576: if ( VR(s) != v ) {
577: if ( n )
578: *rp = 0;
579: else
580: *rp = t;
581: return;
582: }
583: } else
584: s = (P)t;
585: for ( dc = DC(s); dc && cmpq(DEG(dc),(Q)n); dc = NEXT(dc) );
586: if ( dc )
587: *rp = (Obj)COEF(dc);
588: else
589: *rp = 0;
590: }
591: }
592:
1.16 noro 593: void Pcoef(NODE arg,Obj *rp)
1.1 noro 594: {
595: Obj t,n;
596: P s;
597: DCP dc;
598: int id;
599: V v;
600:
601: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
602: *rp = 0;
603: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
604: *rp = 0;
605: else if ( id == O_N ) {
606: if ( !n )
607: *rp = t;
608: else
609: *rp = 0;
610: } else {
611: if ( argc(arg) == 3 ) {
612: if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
613: getcoef(CO,(P)t,v,(Q)n,(P *)rp); return;
614: } else
615: s = (P)t;
616: if ( VR(s) != v ) {
617: if ( n )
618: *rp = 0;
619: else
620: *rp = t;
621: return;
622: }
623: } else
624: s = (P)t;
625: for ( dc = DC(s); dc && cmpq(DEG(dc),(Q)n); dc = NEXT(dc) );
626: if ( dc )
627: *rp = (Obj)COEF(dc);
628: else
629: *rp = 0;
630: }
631: }
632:
1.16 noro 633: void Pcoef_gf2n(NODE arg,Obj *rp)
1.1 noro 634: {
635: Obj t,n;
636: int id,d;
637: UP2 up2;
638:
639: if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
640: *rp = 0;
641: else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
642: *rp = 0;
643: else if ( id == O_N && NID((Num)t) == N_GF2N ) {
644: d = QTOS((Q)n);
645: up2 = ((GF2N)t)->body;
646: if ( d > degup2(up2) )
647: *rp = 0;
648: else
649: *rp = (Obj)(up2->b[d/BSH]&(((unsigned long)1)<<(d%BSH))?ONE:0);
650: } else
651: *rp = 0;
652: }
653:
1.16 noro 654: void Pdeg(NODE arg,Q *rp)
1.1 noro 655: {
656: Obj t,v;
657: int d;
658:
659: #if 0
660: if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
661: *rp = 0;
662: else if ( !(v = (Obj)ARG1(arg)) || (VR((P)v) != VR((P)t)) )
663: *rp = 0;
664: else
665: *rp = (Obj)DEG(DC((P)t));
666: #endif
667: if ( !(t = (Obj)ARG0(arg)) )
668: STOQ(-1,*rp);
669: else if ( OID(t) != O_P ) {
670: if ( OID(t) == O_N && NID(t) == N_GF2N
671: && (v=(Obj)ARG1(arg)) && OID(v)== O_N && NID(v) == N_GF2N ) {
672: d = degup2(((GF2N)t)->body);
673: STOQ(d,*rp);
674: } else
675: *rp = 0;
676: } else
677: degp(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
678: }
679:
1.16 noro 680: void Pmindeg(NODE arg,Q *rp)
1.1 noro 681: {
682: Obj t;
683:
684: if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
685: *rp = 0;
686: else
687: getmindeg(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
688: }
689:
1.16 noro 690: void Psetmod(NODE arg,Q *rp)
1.1 noro 691: {
692: if ( arg ) {
693: asir_assert(ARG0(arg),O_N,"setmod");
694: current_mod = QTOS((Q)ARG0(arg));
695: }
696: STOQ(current_mod,*rp);
697: }
698:
1.16 noro 699: void Psparsemod_gf2n(NODE arg,Q *rp)
1.1 noro 700: {
701: int id;
702:
703: if ( arg && current_mod_gf2n )
704: current_mod_gf2n->id = ARG0(arg)?1:0;
705: if ( !current_mod_gf2n )
706: id = -1;
707: else
708: id = current_mod_gf2n->id;
709: STOQ(id,*rp);
710: }
711:
1.16 noro 712: void Pmultest_gf2n(NODE arg,GF2N *rp)
1.1 noro 713: {
714: GF2N a,b,c;
715: int i;
716:
717: a = (GF2N)ARG0(arg); b = (GF2N)ARG0(arg);
718: for ( i = 0; i < 10000; i++ )
719: mulgf2n(a,b,&c);
720: *rp = c;
721: }
722:
1.16 noro 723: void Psquaretest_gf2n(NODE arg,GF2N *rp)
1.1 noro 724: {
725: GF2N a,c;
726: int i;
727:
728: a = (GF2N)ARG0(arg);
729: for ( i = 0; i < 10000; i++ )
730: squaregf2n(a,&c);
731: *rp = c;
732: }
733:
1.16 noro 734: void Pinvtest_gf2n(NODE arg,GF2N *rp)
1.1 noro 735: {
736: GF2N a,c;
737: int i;
738:
739: a = (GF2N)ARG0(arg);
740: for ( i = 0; i < 10000; i++ )
741: invgf2n(a,&c);
742: *rp = c;
743: }
744:
1.16 noro 745: void Pbininv_gf2n(NODE arg,GF2N *rp)
1.1 noro 746: {
747: UP2 a,inv;
748: int n;
749:
750: a = ((GF2N)ARG0(arg))->body;
751: n = QTOS((Q)ARG1(arg));
752: type1_bin_invup2(a,n,&inv);
753: MKGF2N(inv,*rp);
754: }
755:
1.16 noro 756: void Prinvtest_gf2n(Real *rp)
1.1 noro 757: {
758: GF2N *a;
759: GF2N c;
760: double t0,t1,r;
761: int i;
762: double get_clock();
763:
764: a = (GF2N *)ALLOCA(1000*sizeof(GF2N));
765: for ( i = 0; i < 1000; i++ ) {
766: randomgf2n(&a[i]);
767: }
768: t0 = get_clock();
769: for ( i = 0; i < 1000; i++ )
770: invgf2n(a[i],&c);
771: t1 = get_clock();
772: r = (t1-t0)/1000;
773: MKReal(r,*rp);
774: }
775:
1.16 noro 776: void Pfind_root_gf2n(NODE arg,GF2N *rp)
1.1 noro 777: {
778:
779: #if 0
780: UP p;
781:
782: ptoup((P)ARG0(arg),&p);
783: find_root_gf2n(p,rp);
784: #else
785: UP2 p;
786:
787: ptoup2((P)ARG0(arg),&p);
788: find_root_up2(p,rp);
789: #endif
790: }
791:
1.16 noro 792: void Pis_irred_gf2(NODE arg,Q *rp)
1.1 noro 793: {
794: UP2 t;
795:
796: ptoup2(ARG0(arg),&t);
797: *rp = irredcheckup2(t) ? ONE : 0;
798: }
799:
1.16 noro 800: void Pis_irred_ddd_gf2(NODE arg,Q *rp)
1.1 noro 801: {
802: UP2 t;
803: int ret;
804:
805: ptoup2(ARG0(arg),&t);
806: ret = irredcheck_dddup2(t);
807: STOQ(ret,*rp);
808: }
809:
1.16 noro 810: void Psetmod_ff(NODE arg,Obj *rp)
1.1 noro 811: {
812: int ac;
1.7 noro 813: int d;
1.1 noro 814: Obj mod,defpoly;
815: N n;
816: UP up;
817: UP2 up2;
1.14 noro 818: UM dp;
1.6 noro 819: Q q,r;
1.15 noro 820: P p,p1,y;
1.1 noro 821: NODE n0,n1;
822: LIST list;
823:
824: ac = argc(arg);
825: if ( ac == 1 ) {
826: mod = (Obj)ARG0(arg);
827: if ( !mod )
828: error("setmod_ff : invalid argument");
829: switch ( OID(mod) ) {
830: case O_N:
1.7 noro 831: current_ff = FF_GFP;
832: setmod_lm(NM((Q)mod));
1.6 noro 833: break;
1.1 noro 834: case O_P:
835: current_ff = FF_GF2N;
836: setmod_gf2n((P)mod); break;
837: default:
838: error("setmod_ff : invalid argument");
839: }
840: } else if ( ac == 2 ) {
1.7 noro 841: if ( OID(ARG0(arg)) == O_N ) {
842: /* small finite field; primitive root representation */
843: current_ff = FF_GFS;
844: setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
1.6 noro 845: } else {
1.7 noro 846: mod = (Obj)ARG1(arg);
1.6 noro 847: current_ff = FF_GFPN;
848: defpoly = (Obj)ARG0(arg);
849: if ( !mod || !defpoly )
850: error("setmod_ff : invalid argument");
851: setmod_lm(NM((Q)mod));
852: setmod_gfpn((P)defpoly);
853: }
1.14 noro 854: } else if ( ac == 3 ) {
855: /* finite extension of a small finite field */
856: current_ff = FF_GFS;
857: setmod_sf(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)));
858: d = QTOS((Q)ARG2(arg));
859: generate_defpoly_sfum(d,&dp);
1.15 noro 860: setmod_gfsn(dp);
861: current_ff = FF_GFSN;
1.1 noro 862: }
863: switch ( current_ff ) {
864: case FF_GFP:
865: getmod_lm(&n); NTOQ(n,1,q); *rp = (Obj)q; break;
866: case FF_GF2N:
867: getmod_gf2n(&up2); up2top(up2,&p); *rp = (Obj)p; break;
868: case FF_GFPN:
869: getmod_lm(&n); NTOQ(n,1,q);
870: getmod_gfpn(&up); uptop(up,&p);
871: MKNODE(n1,q,0); MKNODE(n0,p,n1);
872: MKLIST(list,n0);
873: *rp = (Obj)list; break;
1.6 noro 874: case FF_GFS:
1.15 noro 875: case FF_GFSN:
1.7 noro 876: STOQ(current_gfs_p,q);
1.15 noro 877: if ( current_gfs_ext )
1.7 noro 878: enc_to_p(current_gfs_p,current_gfs_iton[1],
879: VR(current_gfs_ext),&p);
1.15 noro 880: else {
1.11 noro 881: if ( current_gfs_p == 2 )
882: r = ONE;
1.17 noro 883: else if ( !current_gfs_ntoi )
884: r = 0;
1.11 noro 885: else
886: STOQ(current_gfs_iton[1],r);
1.15 noro 887: p = (P)r;
888: }
889: switch ( current_ff ) {
890: case FF_GFS:
891: n0 = mknode(3,q,current_gfs_ext,p);
892: break;
893: case FF_GFSN:
894: getmod_gfsn(&dp);
895: makevar("y",&y);
896: sfumtop(VR(y),dp,&p1);
897: n0 = mknode(4,q,current_gfs_ext,p,p1);
898: break;
1.7 noro 899: }
1.6 noro 900: MKLIST(list,n0);
901: *rp = (Obj)list; break;
1.1 noro 902: default:
903: *rp = 0; break;
904: }
905: }
906:
1.16 noro 907: void Pextdeg_ff(Q *rp)
1.1 noro 908: {
909: int d;
910: UP2 up2;
911: UP up;
1.15 noro 912: UM dp;
1.1 noro 913:
914: switch ( current_ff ) {
915: case FF_GFP:
916: *rp = ONE; break;
917: case FF_GF2N:
918: getmod_gf2n(&up2); d = degup2(up2); STOQ(d,*rp); break;
919: case FF_GFPN:
920: getmod_gfpn(&up); STOQ(up->d,*rp); break;
1.10 noro 921: case FF_GFS:
922: if ( !current_gfs_ext )
923: *rp = ONE;
924: else
925: *rp = DEG(DC(current_gfs_ext));
926: break;
1.15 noro 927: case FF_GFSN:
928: getmod_gfsn(&dp);
929: STOQ(DEG(dp),*rp);
930: break;
1.1 noro 931: default:
932: error("extdeg_ff : current_ff is not set");
933: }
934: }
935:
1.16 noro 936: void Pcharacteristic_ff(Q *rp)
1.1 noro 937: {
938: N lm;
939:
940: switch ( current_ff ) {
941: case FF_GFP:
942: case FF_GFPN:
943: getmod_lm(&lm); NTOQ(lm,1,*rp); break;
944: case FF_GF2N:
945: STOQ(2,*rp); break;
1.10 noro 946: case FF_GFS:
1.15 noro 947: case FF_GFSN:
1.10 noro 948: STOQ(current_gfs_p,*rp); break;
1.1 noro 949: default:
950: error("characteristic_ff : current_ff is not set");
951: }
952: }
953:
1.16 noro 954: void Pfield_type_ff(Q *rp)
1.1 noro 955: {
956: STOQ(current_ff,*rp);
957: }
958:
1.16 noro 959: void Pfield_order_ff(Q *rp)
1.1 noro 960: {
961: N n;
962:
963: field_order_ff(&n);
964: NTOQ(n,1,*rp);
965: }
966:
1.16 noro 967: void Prandom_ff(Obj *rp)
1.1 noro 968: {
969: LM l;
970: GF2N g;
971: GFPN p;
1.7 noro 972: GFS s;
1.15 noro 973: GFSN spn;
1.1 noro 974:
975: switch ( current_ff ) {
976: case FF_GFP:
977: random_lm(&l); *rp = (Obj)l; break;
978: case FF_GF2N:
979: randomgf2n(&g); *rp = (Obj)g; break;
980: case FF_GFPN:
981: randomgfpn(&p); *rp = (Obj)p; break;
1.7 noro 982: case FF_GFS:
983: randomgfs(&s); *rp = (Obj)s; break;
1.15 noro 984: case FF_GFSN:
985: randomgfsn(&spn); *rp = (Obj)spn; break;
1.1 noro 986: default:
987: error("random_ff : current_ff is not set");
988: }
989: }
990:
1.16 noro 991: void Psimp_ff(NODE arg,Obj *rp)
992: {
1.1 noro 993: simp_ff((Obj)ARG0(arg),rp);
994: }
995:
1.16 noro 996: void getcoef(VL vl,P p,V v,Q d,P *r)
1.1 noro 997: {
998: P s,t,u,a,b,x;
999: DCP dc;
1000: V w;
1001:
1002: if ( !p )
1003: *r = 0;
1004: else if ( NUM(p) )
1005: *r = d ? 0 : p;
1006: else if ( (w=VR(p)) == v ) {
1007: for ( dc = DC(p); dc && cmpq(DEG(dc),d); dc = NEXT(dc) );
1008: *r = dc ? COEF(dc) : 0;
1009: } else {
1010: MKV(w,x);
1011: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1012: getcoef(vl,COEF(dc),v,d,&t);
1013: if ( t ) {
1014: pwrp(vl,x,DEG(dc),&u); mulp(vl,t,u,&a);
1015: addp(vl,s,a,&b); s = b;
1016: }
1017: }
1018: *r = s;
1019: }
1020: }
1021:
1.16 noro 1022: void Pdeglist(NODE arg,LIST *rp)
1.1 noro 1023: {
1024: NODE d;
1025:
1026: getdeglist((P)ARG0(arg),VR((P)ARG1(arg)),&d);
1027: MKLIST(*rp,d);
1028: }
1029:
1.16 noro 1030: void Pch_mv(NODE arg,P *rp)
1.1 noro 1031: {
1032: change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
1033: }
1034:
1.16 noro 1035: void Pre_mv(NODE arg,P *rp)
1.1 noro 1036: {
1037: restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
1038: }
1039:
1.16 noro 1040: void change_mvar(VL vl,P p,V v,P *r)
1.1 noro 1041: {
1042: Q d;
1043: DCP dc,dc0;
1044: NODE dl;
1045:
1046: if ( !p || NUM(p) || (VR(p) == v) )
1047: *r = p;
1048: else {
1049: getdeglist(p,v,&dl);
1050: for ( dc0 = 0; dl; dl = NEXT(dl) ) {
1051: NEXTDC(dc0,dc); DEG(dc) = d = (Q)BDY(dl);
1052: getcoef(vl,p,v,d,&COEF(dc));
1053: }
1054: NEXT(dc) = 0; MKP(v,dc0,*r);
1055: }
1056: }
1057:
1.16 noro 1058: void restore_mvar(VL vl,P p,V v,P *r)
1.1 noro 1059: {
1060: P s,u,a,b,x;
1061: DCP dc;
1062:
1063: if ( !p || NUM(p) || (VR(p) != v) )
1064: *r = p;
1065: else {
1066: MKV(v,x);
1067: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
1068: pwrp(vl,x,DEG(dc),&u); mulp(vl,COEF(dc),u,&a);
1069: addp(vl,s,a,&b); s = b;
1070: }
1071: *r = s;
1072: }
1073: }
1074:
1.16 noro 1075: void getdeglist(P p,V v,NODE *d)
1.1 noro 1076: {
1077: NODE n,n0,d0,d1,d2;
1078: DCP dc;
1079:
1080: if ( !p || NUM(p) ) {
1081: MKNODE(n,0,0); *d = n;
1082: } else if ( VR(p) == v ) {
1083: for ( n0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
1084: NEXTNODE(n0,n); BDY(n) = (pointer)DEG(dc);
1085: }
1086: NEXT(n) = 0; *d = n0;
1087: } else {
1088: for ( dc = DC(p), d0 = 0; dc; dc = NEXT(dc) ) {
1089: getdeglist(COEF(dc),v,&d1); mergedeglist(d0,d1,&d2); d0 = d2;
1090: }
1091: *d = d0;
1092: }
1093: }
1.16 noro 1094:
1095: void Pmergelist(NODE arg,LIST *rp)
1.1 noro 1096: {
1097: NODE n;
1098:
1099: asir_assert(ARG0(arg),O_LIST,"mergelist");
1100: asir_assert(ARG1(arg),O_LIST,"mergelist");
1101: mergedeglist(BDY((LIST)ARG0(arg)),BDY((LIST)ARG1(arg)),&n);
1102: MKLIST(*rp,n);
1103: }
1104:
1.16 noro 1105: void mergedeglist(NODE d0,NODE d1,NODE *dr)
1.1 noro 1106: {
1107: NODE t0,t,dt;
1108: Q d;
1109: int c;
1110:
1111: if ( !d0 )
1112: *dr = d1;
1113: else {
1114: while ( d1 ) {
1115: dt = d1; d1 = NEXT(d1); d = (Q)BDY(dt);
1116: for ( t0 = 0, t = d0; t; t0 = t, t = NEXT(t) ) {
1117: c = cmpq(d,(Q)BDY(t));
1118: if ( !c )
1119: break;
1120: else if ( c > 0 ) {
1121: if ( !t0 ) {
1122: NEXT(dt) = d0; d0 = dt;
1123: } else {
1124: NEXT(t0) = dt; NEXT(dt) = t;
1125: }
1126: break;
1127: }
1128: }
1129: if ( !t ) {
1130: NEXT(t0) = dt; *dr = d0; return;
1131: }
1132: }
1133: *dr = d0;
1134: }
1135: }
1136:
1.16 noro 1137: void Pptomp(NODE arg,P *rp)
1.1 noro 1138: {
1139: ptomp(QTOS((Q)ARG1(arg)),(P)ARG0(arg),rp);
1140: }
1141:
1.16 noro 1142: void Pmptop(NODE arg,P *rp)
1.1 noro 1143: {
1144: mptop((P)ARG0(arg),rp);
1145: }
1146:
1.16 noro 1147: void Pptolmp(NODE arg,P *rp)
1.1 noro 1148: {
1149: ptolmp((P)ARG0(arg),rp);
1150: }
1151:
1.16 noro 1152: void Plmptop(NODE arg,P *rp)
1.1 noro 1153: {
1154: lmptop((P)ARG0(arg),rp);
1.11 noro 1155: }
1156:
1.16 noro 1157: void Psf_galois_action(NODE arg,P *rp)
1.11 noro 1158: {
1159: sf_galois_action(ARG0(arg),ARG1(arg),rp);
1.12 noro 1160: }
1161:
1162: /*
1163: sf_embed(F,B,PM)
1164: F : an element of GF(pn)
1165: B : the image of the primitive root of GF(pn)
1166: PM : order of GF(pm)
1167: */
1168:
1.16 noro 1169: void Psf_embed(NODE arg,P *rp)
1.12 noro 1170: {
1171: int k,pm;
1172:
1173: /* GF(pn)={0,1,a,a^2,...}->GF(pm)={0,1,b,b^2,...}; a->b^k */
1174: k = CONT((GFS)ARG1(arg));
1175: pm = QTOS((Q)ARG2(arg));
1176: sf_embed((P)ARG0(arg),k,pm,rp);
1.13 noro 1177: }
1178:
1.16 noro 1179: void Psf_log(NODE arg,Q *rp)
1.13 noro 1180: {
1181: int k;
1182:
1183: if ( !ARG0(arg) )
1184: error("sf_log : invalid armument");
1185: k = CONT((GFS)ARG0(arg));
1186: STOQ(k,*rp);
1.12 noro 1187: }
1188:
1.16 noro 1189: void Psf_find_root(NODE arg,GFS *rp)
1.12 noro 1190: {
1191: P p;
1192: Obj t;
1193: int d;
1194: UM u;
1195: int *root;
1196:
1197: p = (P)ARG0(arg);
1198: simp_ff((Obj)p,&t); p = (P)t;
1199: d = getdeg(VR(p),p);
1200: u = W_UMALLOC(d);
1201: ptosfum(p,u);
1202: root = (int *)ALLOCA(d*sizeof(int));
1203: find_rootsf(u,root);
1204: MKGFS(IFTOF(root[0]),*rp);
1205: }
1206:
1.16 noro 1207: void Psf_minipoly(NODE arg,P *rp)
1.12 noro 1208: {
1209: Obj t;
1210: P p1,p2;
1211: int d1,d2;
1212: UM up1,up2,m;
1213:
1214: p1 = (P)ARG0(arg); simp_ff((Obj)p1,&t); p1 = (P)t;
1215: p2 = (P)ARG1(arg); simp_ff((Obj)p2,&t); p2 = (P)t;
1216: d1 = getdeg(VR(p1),p1); up1 = W_UMALLOC(d1); ptosfum(p1,up1);
1217: d2 = getdeg(VR(p2),p2); up2 = W_UMALLOC(d2); ptosfum(p2,up2);
1218: m = W_UMALLOC(d2);
1219: minipolysf(up1,up2,m);
1220: sfumtop(VR(p2),m,&p1);
1221: sfptop(p1,rp);
1.11 noro 1222: }
1223:
1.16 noro 1224: void Pptosfp(NODE arg,P *rp)
1.11 noro 1225: {
1226: ptosfp(ARG0(arg),rp);
1.6 noro 1227: }
1228:
1.16 noro 1229: void Psfptop(NODE arg,P *rp)
1.6 noro 1230: {
1231: sfptop((P)ARG0(arg),rp);
1.18 noro 1232: }
1233:
1234: void Psfptopsfp(NODE arg,P *rp)
1235: {
1236: sfptopsfp((P)ARG0(arg),VR((P)ARG1(arg)),rp);
1.1 noro 1237: }
1238:
1.16 noro 1239: void Pptogf2n(NODE arg,GF2N *rp)
1.1 noro 1240: {
1241: ptogf2n((Obj)ARG0(arg),rp);
1242: }
1243:
1.16 noro 1244: void Pgf2ntop(NODE arg,P *rp)
1.1 noro 1245: {
1246: if ( argc(arg) == 2 )
1247: up2_var = VR((P)ARG1(arg));
1248: gf2ntop((GF2N)ARG0(arg),rp);
1249: }
1250:
1.16 noro 1251: void Pgf2ntovect(NODE arg,VECT *rp)
1.1 noro 1252: {
1253: gf2ntovect((GF2N)ARG0(arg),rp);
1254: }
1255:
1.16 noro 1256: void Pptogfpn(NODE arg,GFPN *rp)
1.1 noro 1257: {
1258: ptogfpn((Obj)ARG0(arg),rp);
1259: }
1260:
1.16 noro 1261: void Pgfpntop(NODE arg,P *rp)
1.1 noro 1262: {
1263: if ( argc(arg) == 2 )
1264: up_var = VR((P)ARG1(arg));
1265: gfpntop((GFPN)ARG0(arg),rp);
1266: }
1267:
1.16 noro 1268: void Pureverse(NODE arg,P *rp)
1.1 noro 1269: {
1270: UP p,r;
1271:
1272: ptoup((P)ARG0(arg),&p);
1.3 noro 1273: if ( argc(arg) == 1 )
1274: reverseup(p,p->d,&r);
1275: else
1276: reverseup(p,QTOS((Q)ARG1(arg)),&r);
1.1 noro 1277: uptop(r,rp);
1278: }
1279:
1.16 noro 1280: void Putrunc(NODE arg,P *rp)
1.1 noro 1281: {
1282: UP p,r;
1283:
1284: ptoup((P)ARG0(arg),&p);
1285: truncup(p,QTOS((Q)ARG1(arg))+1,&r);
1286: uptop(r,rp);
1287: }
1288:
1.16 noro 1289: void Pudecomp(NODE arg,LIST *rp)
1.1 noro 1290: {
1291: P u,l;
1292: UP p,up,low;
1293: NODE n0,n1;
1294:
1295: ptoup((P)ARG0(arg),&p);
1296: decompup(p,QTOS((Q)ARG1(arg))+1,&low,&up);
1297: uptop(low,&l);
1298: uptop(up,&u);
1299: MKNODE(n1,u,0); MKNODE(n0,l,n1);
1300: MKLIST(*rp,n0);
1301: }
1302:
1.16 noro 1303: void Purembymul(NODE arg,P *rp)
1.1 noro 1304: {
1305: UP p1,p2,r;
1306:
1307: if ( !ARG0(arg) || !ARG1(arg) )
1308: *rp = 0;
1309: else {
1310: ptoup((P)ARG0(arg),&p1);
1311: ptoup((P)ARG1(arg),&p2);
1312: rembymulup(p1,p2,&r);
1313: uptop(r,rp);
1314: }
1315: }
1316:
1317: /*
1318: * d1 = deg(p1), d2 = deg(p2)
1319: * d1 <= 2*d2-1
1320: * p2*inv = 1 mod x^d2
1321: */
1322:
1.16 noro 1323: void Purembymul_precomp(NODE arg,P *rp)
1.1 noro 1324: {
1325: UP p1,p2,inv,r;
1326:
1327: if ( !ARG0(arg) || !ARG1(arg) )
1328: *rp = 0;
1329: else {
1330: ptoup((P)ARG0(arg),&p1);
1331: ptoup((P)ARG1(arg),&p2);
1332: ptoup((P)ARG2(arg),&inv);
1333: if ( p1->d >= 2*p2->d ) {
1334: error("urembymul_precomp : degree of 1st arg is too large");
1335: /* fprintf(stderr,"urembymul_precomp : degree of 1st arg is too large"); */
1336: remup(p1,p2,&r);
1337: } else
1338: hybrid_rembymulup_special(current_ff,p1,p2,inv,&r);
1339: uptop(r,rp);
1340: }
1341: }
1342:
1.16 noro 1343: void Puinvmod(NODE arg,P *rp)
1.1 noro 1344: {
1345: UP p,r;
1346:
1347: ptoup((P)ARG0(arg),&p);
1348: invmodup(p,QTOS((Q)ARG1(arg)),&r);
1349: uptop(r,rp);
1350: }
1351:
1.16 noro 1352: void Purevinvmod(NODE arg,P *rp)
1.1 noro 1353: {
1354: UP p,pr,r;
1355:
1356: ptoup((P)ARG0(arg),&p);
1357: reverseup(p,p->d,&pr);
1358: invmodup(pr,QTOS((Q)ARG1(arg)),&r);
1359: uptop(r,rp);
1360: }
1361:
1.16 noro 1362: void Ppwrmod_ff(NODE arg,P *rp)
1.1 noro 1363: {
1364: UP p1,p2;
1365:
1366: ptoup((P)ARG0(arg),&p1);
1367: switch ( current_ff ) {
1368: case FF_GFP:
1369: hybrid_powermodup(p1,&p2); break;
1370: case FF_GF2N:
1371: powermodup_gf2n(p1,&p2); break;
1372: case FF_GFPN:
1.10 noro 1373: case FF_GFS:
1.15 noro 1374: case FF_GFSN:
1.1 noro 1375: powermodup(p1,&p2); break;
1376: default:
1377: error("pwrmod_ff : current_ff is not set");
1378: }
1379: uptop(p2,rp);
1380: }
1381:
1.16 noro 1382: void Pgeneric_pwrmod_ff(NODE arg,P *rp)
1.1 noro 1383: {
1384: UP g,f,r;
1385:
1386: ptoup((P)ARG0(arg),&g);
1387: ptoup((P)ARG1(arg),&f);
1388: switch ( current_ff ) {
1389: case FF_GFP:
1390: hybrid_generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1391: case FF_GF2N:
1392: generic_powermodup_gf2n(g,f,(Q)ARG2(arg),&r); break;
1393: case FF_GFPN:
1.10 noro 1394: case FF_GFS:
1.15 noro 1395: case FF_GFSN:
1.1 noro 1396: generic_powermodup(g,f,(Q)ARG2(arg),&r); break;
1397: default:
1398: error("generic_pwrmod_ff : current_ff is not set");
1399: }
1400: uptop(r,rp);
1401: }
1402:
1.16 noro 1403: void Ppwrtab_ff(NODE arg,VECT *rp)
1.1 noro 1404: {
1405: UP f,xp;
1406: UP *tab;
1407: VECT r;
1408: int i,d;
1409:
1410: ptoup((P)ARG0(arg),&f);
1411: ptoup((P)ARG1(arg),&xp);
1412: d = f->d;
1413:
1414: tab = (UP *)ALLOCA(d*sizeof(UP));
1415: switch ( current_ff ) {
1416: case FF_GFP:
1417: hybrid_powertabup(f,xp,tab); break;
1418: case FF_GF2N:
1419: powertabup_gf2n(f,xp,tab); break;
1420: case FF_GFPN:
1.10 noro 1421: case FF_GFS:
1.15 noro 1422: case FF_GFSN:
1.1 noro 1423: powertabup(f,xp,tab); break;
1424: default:
1425: error("pwrtab_ff : current_ff is not set");
1426: }
1427: MKVECT(r,d); *rp = r;
1428: for ( i = 0; i < d; i++ )
1429: uptop(tab[i],(P *)&BDY(r)[i]);
1430: }
1431:
1.16 noro 1432: void Pkpwrmod_lm(NODE arg,P *rp)
1.1 noro 1433: {
1434: UP p1,p2;
1435:
1436: ptoup((P)ARG0(arg),&p1);
1437: powermodup(p1,&p2);
1438: uptop(p2,rp);
1439: }
1440:
1.16 noro 1441: void Pkgeneric_pwrmod_lm(NODE arg,P *rp)
1.1 noro 1442: {
1443: UP g,f,r;
1444:
1445: ptoup((P)ARG0(arg),&g);
1446: ptoup((P)ARG1(arg),&f);
1447: generic_powermodup(g,f,(Q)ARG2(arg),&r);
1448: uptop(r,rp);
1449: }
1450:
1.16 noro 1451: void Pkpwrtab_lm(NODE arg,VECT *rp)
1.1 noro 1452: {
1453: UP f,xp;
1454: UP *tab;
1455: VECT r;
1456: int i,d;
1457:
1458: ptoup((P)ARG0(arg),&f);
1459: ptoup((P)ARG1(arg),&xp);
1460: d = f->d;
1461:
1462: tab = (UP *)ALLOCA(d*sizeof(UP));
1463: powertabup(f,xp,tab);
1464: MKVECT(r,d); *rp = r;
1465: for ( i = 0; i < d; i++ )
1466: uptop(tab[i],(P *)&BDY(r)[i]);
1467: }
1468:
1.16 noro 1469: void Plazy_lm(NODE arg,Q *rp)
1.1 noro 1470: {
1471: lm_lazy = QTOS((Q)ARG0(arg));
1472: *rp = 0;
1473: }
1474:
1.16 noro 1475: void Pkmul(NODE arg,P *rp)
1.1 noro 1476: {
1477: P n1,n2;
1478:
1479: n1 = (P)ARG0(arg); n2 = (P)ARG1(arg);
1480: asir_assert(n1,O_P,"kmul");
1481: asir_assert(n2,O_P,"kmul");
1482: kmulp(CO,n1,n2,rp);
1483: }
1484:
1.16 noro 1485: void Pksquare(NODE arg,P *rp)
1.1 noro 1486: {
1487: P n1;
1488:
1489: n1 = (P)ARG0(arg);
1490: asir_assert(n1,O_P,"ksquare");
1491: ksquarep(CO,n1,rp);
1492: }
1493:
1.16 noro 1494: void Pktmul(NODE arg,P *rp)
1.1 noro 1495: {
1496: UP p1,p2,r;
1497:
1498: ptoup((P)ARG0(arg),&p1);
1499: ptoup((P)ARG1(arg),&p2);
1500: tkmulup(p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1501: uptop(r,rp);
1502: }
1503:
1.16 noro 1504: void Pumul(NODE arg,P *rp)
1.1 noro 1505: {
1506: P a1,a2;
1507: UP p1,p2,r;
1508:
1509: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1510: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1511: mulp(CO,a1,a2,rp);
1512: else {
1.16 noro 1513: if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
1.1 noro 1514: error("umul : invalid argument");
1515: ptoup(a1,&p1);
1516: ptoup(a2,&p2);
1517: hybrid_mulup(0,p1,p2,&r);
1518: uptop(r,rp);
1519: }
1520: }
1521:
1.16 noro 1522: void Pusquare(NODE arg,P *rp)
1.1 noro 1523: {
1.16 noro 1524: UP p1,r;
1.1 noro 1525:
1526: ptoup((P)ARG0(arg),&p1);
1527: hybrid_squareup(0,p1,&r);
1528: uptop(r,rp);
1529: }
1530:
1.16 noro 1531: void Putmul(NODE arg,P *rp)
1.1 noro 1532: {
1533: UP p1,p2,r;
1534:
1535: ptoup((P)ARG0(arg),&p1);
1536: ptoup((P)ARG1(arg),&p2);
1537: hybrid_tmulup(0,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1538: uptop(r,rp);
1539: }
1540:
1.16 noro 1541: void Pumul_ff(NODE arg,Obj *rp)
1.1 noro 1542: {
1543: P a1,a2;
1544: UP p1,p2,r;
1545: P p;
1546:
1547: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1548: ptoup(a1,&p1);
1549: ptoup(a2,&p2);
1550: hybrid_mulup(current_ff,p1,p2,&r);
1551: uptop(r,&p);
1552: simp_ff((Obj)p,rp);
1553: }
1554:
1.16 noro 1555: void Pusquare_ff(NODE arg,Obj *rp)
1.1 noro 1556: {
1.16 noro 1557: UP p1,r;
1.1 noro 1558: P p;
1559:
1560: ptoup((P)ARG0(arg),&p1);
1561: hybrid_squareup(current_ff,p1,&r);
1562: uptop(r,&p);
1563: simp_ff((Obj)p,rp);
1564: }
1565:
1.16 noro 1566: void Putmul_ff(NODE arg,Obj *rp)
1.1 noro 1567: {
1568: UP p1,p2,r;
1569: P p;
1570:
1571: ptoup((P)ARG0(arg),&p1);
1572: ptoup((P)ARG1(arg),&p2);
1573: hybrid_tmulup(current_ff,p1,p2,QTOS((Q)ARG2(arg))+1,&r);
1574: uptop(r,&p);
1575: simp_ff((Obj)p,rp);
1576: }
1577:
1.16 noro 1578: void Phfmul_lm(NODE arg,P *rp)
1.1 noro 1579: {
1580: UP p1,p2,r;
1581: UP hi,lo,mid,t,s,p10,p11,p20,p21;
1582: unsigned int m,d;
1583: int i;
1584:
1585: ptoup((P)ARG0(arg),&p1);
1586: ptoup((P)ARG1(arg),&p2);
1587: d = MAX(p1->d,p2->d);
1588: for ( m = 1; m < d; m <<= 1 );
1589: if ( m > d )
1590: m >>= 1;
1591: if ( d-m < 10000 ) {
1592: decompup(p1,m,&p10,&p11); /* p1 = p11*x^m+p10 */
1593: decompup(p2,m,&p20,&p21); /* p2 = p21*x^m+p20 */
1594: fft_mulup_lm(p10,p20,&lo);
1595: kmulup(p11,p21,&hi);
1596: kmulup(p11,p20,&t); kmulup(p10,p21,&s); addup(t,s,&mid);
1597: r = UPALLOC(2*d);
1598: bzero((char *)COEF(r),(2*d+1)*sizeof(Num));
1599: if ( lo )
1600: bcopy((char *)COEF(lo),(char *)COEF(r),(DEG(lo)+1)*sizeof(Num));
1601: if ( hi )
1602: bcopy((char *)COEF(hi),(char *)(COEF(r)+2*m),(DEG(hi)+1)*sizeof(Num));
1603: for ( i = 2*d; i >= 0 && !COEF(r)[i]; i-- );
1604: if ( i < 0 )
1605: r = 0;
1606: else {
1607: DEG(r) = i;
1608: t = UPALLOC(DEG(mid)+m); DEG(t) = DEG(mid)+m;
1609: bzero((char *)COEF(t),(DEG(mid)+m+1)*sizeof(Num));
1610: bcopy((char *)COEF(mid),(char *)(COEF(t)+m),(DEG(mid)+1)*sizeof(Num));
1611: addup(t,r,&s);
1612: r = s;
1613: }
1614: } else
1615: fft_mulup_lm(p1,p2,&r);
1616: uptop(r,rp);
1617: }
1618:
1.16 noro 1619: void Pfmultest(NODE arg,LIST *rp)
1.1 noro 1620: {
1621: P p1,p2,r;
1622: int d1,d2;
1623: UM w1,w2,wr;
1624: unsigned int *f1,*f2,*fr,*w;
1625: int index,mod,root,d,maxint,i;
1626: int cond;
1627: Q prime;
1628: NODE n0,n1;
1629:
1630: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); index = QTOS((Q)ARG2(arg));
1631: FFT_primes(index,&mod,&root,&d);
1632: maxint = 1<<d;
1633: d1 = UDEG(p1); d2 = UDEG(p2);
1634: if ( maxint < d1+d2+1 )
1635: *rp = 0;
1636: else {
1637: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1638: wr = W_UMALLOC(d1+d2);
1639: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1640: f1 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1641: f2 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1642: fr = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
1643: w = (unsigned int *)ALLOCA(12*maxint*sizeof(unsigned int));
1644:
1645: for ( i = 0; i <= d1; i++ )
1646: f1[i] = (unsigned int)w1->c[i];
1647: for ( i = 0; i <= d2; i++ )
1648: f2[i] = (unsigned int)w2->c[i];
1649:
1650: cond = FFT_pol_product(d1,f1,d2,f2,fr,index,w);
1651: if ( cond )
1652: error("fmultest : ???");
1653: wr->d = d1+d2;
1654: for ( i = 0; i <= d1+d2; i++ )
1655: wr->c[i] = (unsigned int)fr[i];
1656: umtop(VR(p1),wr,&r);
1657: STOQ(mod,prime);
1658: MKNODE(n1,prime,0);
1659: MKNODE(n0,r,n1);
1660: MKLIST(*rp,n0);
1661: }
1662: }
1663:
1.16 noro 1664: void Pkmulum(NODE arg,P *rp)
1.1 noro 1665: {
1666: P p1,p2;
1667: int d1,d2,mod;
1668: UM w1,w2,wr;
1669:
1670: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = QTOS((Q)ARG2(arg));
1671: d1 = UDEG(p1); d2 = UDEG(p2);
1672: w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
1673: wr = W_UMALLOC(d1+d2);
1674: ptoum(mod,p1,w1); ptoum(mod,p2,w2);
1675: kmulum(mod,w1,w2,wr);
1676: umtop(VR(p1),wr,rp);
1677: }
1678:
1.16 noro 1679: void Pksquareum(NODE arg,P *rp)
1.1 noro 1680: {
1681: P p1;
1682: int d1,mod;
1683: UM w1,wr;
1684:
1685: p1 = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
1686: d1 = UDEG(p1);
1687: w1 = W_UMALLOC(d1);
1688: wr = W_UMALLOC(2*d1);
1689: ptoum(mod,p1,w1);
1690: kmulum(mod,w1,w1,wr);
1691: umtop(VR(p1),wr,rp);
1692: }
1693:
1.16 noro 1694: void Ptracemod_gf2n(NODE arg,P *rp)
1.1 noro 1695: {
1696: UP g,f,r;
1697:
1698: ptoup((P)ARG0(arg),&g);
1699: ptoup((P)ARG1(arg),&f);
1700: tracemodup_gf2n(g,f,(Q)ARG2(arg),&r);
1701: uptop(r,rp);
1.2 noro 1702: }
1703:
1.16 noro 1704: void Pumul_specialmod(NODE arg,P *rp)
1.2 noro 1705: {
1706: P a1,a2;
1707: UP p1,p2,r;
1708: int i,nmod;
1709: int *modind;
1710: NODE t,n;
1711:
1712: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1713: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1714: mulp(CO,a1,a2,rp);
1715: else {
1.16 noro 1716: if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
1.2 noro 1717: error("umul_specialmod : invalid argument");
1718: ptoup(a1,&p1);
1719: ptoup(a2,&p2);
1720: n = BDY((LIST)ARG2(arg));
1721: nmod = length(n);
1722: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1723: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1724: modind[i] = QTOS((Q)BDY(t));
1725: fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
1726: uptop(r,rp);
1727: }
1728: }
1729:
1.16 noro 1730: void Pusquare_specialmod(NODE arg,P *rp)
1.2 noro 1731: {
1732: P a1;
1733: UP p1,r;
1734: int i,nmod;
1735: int *modind;
1736: NODE t,n;
1737:
1738: a1 = (P)ARG0(arg);
1739: if ( !a1 || NUM(a1) )
1740: mulp(CO,a1,a1,rp);
1741: else {
1.16 noro 1742: if ( !uzpcheck((Obj)a1) )
1.2 noro 1743: error("usquare_specialmod : invalid argument");
1744: ptoup(a1,&p1);
1745: n = BDY((LIST)ARG1(arg));
1746: nmod = length(n);
1747: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1748: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1749: modind[i] = QTOS((Q)BDY(t));
1750: fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
1751: uptop(r,rp);
1752: }
1753: }
1754:
1.16 noro 1755: void Putmul_specialmod(NODE arg,P *rp)
1.2 noro 1756: {
1757: P a1,a2;
1758: UP p1,p2,r;
1759: int i,nmod;
1760: int *modind;
1761: NODE t,n;
1762:
1763: a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
1764: if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
1765: mulp(CO,a1,a2,rp);
1766: else {
1.16 noro 1767: if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
1.2 noro 1768: error("utmul_specialmod : invalid argument");
1769: ptoup(a1,&p1);
1770: ptoup(a2,&p2);
1771: n = BDY((LIST)ARG3(arg));
1772: nmod = length(n);
1773: modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
1774: for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
1775: modind[i] = QTOS((Q)BDY(t));
1776: fft_mulup_specialmod_main(p1,p2,QTOS((Q)ARG2(arg))+1,modind,nmod,&r);
1777: uptop(r,rp);
1778: }
1.1 noro 1779: }