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