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