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