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