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