![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gf.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 gf.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.3 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.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gf.c,v 1.6 2001/06/25 01:35:20 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: struct resf_dlist {
54: int ib,id;
55: };
56:
57: int resf_degtest(int,int *,int,struct resf_dlist *);
58: void uhensel(P,NODE,int,int,NODE *);
59: void resf_hensel(int,P,int,P *,ML *);
60: void resf_dtest(P,ML,int,int *,int *,DCP *);
61: void resf_dtest_special(P,ML,int,int *,int *,DCP *);
62: void dtest_special(P,ML,P *);
63: void hensel_special(int,int,P,P *,ML *);
64:
65: void nullspace(UM **,UM,int,int,int *);
66: void nullspace_lm(LM **,int,int *);
67: void nullspace_gf2n(GF2N **,int,int *);
68: void nullspace_gfpn(GFPN **,int,int *);
1.5 noro 69: void nullspace_gfs(GFS **,int,int *);
1.1 noro 70: void null_to_sol(int **,int *,int,int,UM *);
71:
72: void showgfmat(UM **,int);
73: void pwr_mod(P,P,V,P,int,N,P *);
74: void rem_mod(P,P,V,P,int,P *);
75:
76: void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
1.6 noro 77: void Psfuhensel();
1.1 noro 78:
79: void Pnullspace_ff();
80:
81: void Psolve_linear_equation_gf2n();
82: void Plinear_form_to_vect(),Pvect_to_linear_form();
83:
84: void solve_linear_equation_gf2n(GF2N **,int,int,int *);
85: void linear_form_to_array(P,VL,int,Num *);
86: void array_to_linear_form(Num *,VL,int,P *);
1.6 noro 87: void sfuhensel(P,V,V,NODE *);
1.1 noro 88:
89: extern int current_ff;
90:
91: struct ftab gf_tab[] = {
92: {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
93: {"nullspace",Pnullspace,3},
94: {"nullspace_ff",Pnullspace_ff,1},
95: /* {"gcda_mod",Pgcda_mod,5}, */
96: {"ftest",Pftest,2},
97: {"resfmain",Presfmain,4},
98: {"pwr_mod",Ppwr_mod,6},
99: {"uhensel",Puhensel,4},
1.6 noro 100: {"sfuhensel",Psfuhensel,3},
1.1 noro 101: {0,0,0},
102: };
103:
104: int resf_degtest(k,in,nb,dlist)
105: int k;
106: int *in;
107: int nb;
108: struct resf_dlist *dlist;
109: {
110: int i,d0,d;
111: int dl_i;
112: struct resf_dlist *t;
113:
114: if ( k < nb )
115: return 0;
116: if ( nb == 1 )
117: return 1;
118: d0 = 0; d = 0; dl_i = 0;
119: for ( i = 0; i < k; i++ ) {
120: t = &dlist[in[i]];
121: if ( t->ib > dl_i + 1 )
122: return 0;
123: else if ( t->ib == dl_i )
124: d += t->id;
125: else if ( !d || (dl_i && d0 != d) )
126: return 0;
127: else {
128: d0 = d; dl_i++; d = t->id;
129: }
130: }
131: if ( dl_i != nb - 1 || d0 != d )
132: return 0;
133: else
134: return 1;
135: }
136:
137: void Puhensel(arg,rp)
138: NODE arg;
139: LIST *rp;
140: {
141: P f;
142: NODE mfl,r;
143: int mod,bound;
144:
145: f = (P)ARG0(arg);
146: mfl = BDY((LIST)ARG1(arg));
147: mod = QTOS((Q)ARG2(arg));
148: bound = QTOS((Q)ARG3(arg));
149: uhensel(f,mfl,mod,bound,&r);
150: MKLIST(*rp,r);
151: }
152:
153: void uhensel(f,mfl,mod,bound,rp)
154: P f;
155: NODE mfl;
156: int mod,bound;
157: NODE *rp;
158: {
159: ML blist,clist,rlist;
160: LUM fl;
161: int nf,i;
162: P s;
163: V v;
164: NODE t,top;
165:
166: nf = length(mfl);
167: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
168: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
169: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
170: ptoum(mod,(P)BDY(t),blist->c[i]);
171: }
172: gcdgen(f,blist,&clist);
173: blist->bound = clist->bound = bound;
174: W_LUMALLOC((int)UDEG(f),bound,fl);
175: ptolum(mod,bound,f,fl);
176: henmain(fl,blist,clist,&rlist);
177: v = VR(f);
178: for ( i = nf-1, top = 0; i >= 0; i-- ) {
179: lumtop(v,mod,bound,rlist->c[i],&s);
180: MKNODE(t,s,top); top = t;
1.6 noro 181: }
182: *rp = top;
183: }
184:
185: void Psfuhensel(arg,rp)
186: NODE arg;
187: LIST *rp;
188: {
189: P f;
190: V x,y;
191: NODE r;
192:
193: f = (P)ARG0(arg);
194: x = VR((P)ARG1(arg));
195: y = VR((P)ARG2(arg));
196: sfuhensel(f,x,y,&r);
197: MKLIST(*rp,r);
198: }
199:
200: void sfuhensel(f,x,y,rp)
201: P f;
202: V x,y;
203: NODE *rp;
204: {
205: ML lift;
206: int i;
1.7 ! noro 207: P s;
1.6 noro 208: NODE t,top;
209:
210: sfhensel(5,f,x,&lift);
211: for ( i = lift->n-1, top = 0; i >= 0; i-- ) {
1.7 ! noro 212: sfbmtop(CO,lift->bound,lift->c[i],x,y,&s);
! 213: MKNODE(t,s,top); top = t;
1.1 noro 214: }
215: *rp = top;
216: }
217:
218: void Presfmain(arg,rp)
219: NODE arg;
220: LIST *rp;
221: {
222: P f;
223: int mod,n,nb,i,j,k;
224: int *nf,*md;
225: P *mf;
226: NODE mfl,mdl,t,s,u;
227: ML list;
228: DCP dc;
229: int sflag;
230:
231: f = (P)ARG0(arg);
232: mod = QTOS((Q)ARG1(arg));
233: mfl = BDY((LIST)ARG2(arg));
234: mdl = BDY((LIST)ARG3(arg));
235: for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
236: for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
237: if ( n == nb ) {
238: /* f must be irreducible! */
239: NEWDC(dc);
240: dc->c = f; dc->d = ONE;
241: } else {
242: nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
243: for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
244: i++, t = NEXT(t), u = NEXT(u) ) {
245: if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
246: sflag = 0;
247: for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
248: mf[k++] = (P)BDY(s);
249: nf[i] = j; md[i] = QTOS((Q)BDY(u));
250: }
251: resf_hensel(mod,f,n,mf,&list);
252: if ( sflag )
253: resf_dtest_special(f,list,nb,nf,md,&dc);
254: else
255: resf_dtest(f,list,nb,nf,md,&dc);
256: }
257: dcptolist(dc,rp);
258: }
259:
260: void resf_hensel(mod,f,nf,mfl,listp)
261: int mod;
262: P f;
263: int nf;
264: P *mfl;
265: ML *listp;
266: {
267: register int i,j;
1.2 noro 268: int q,n,bound,inv,lc;
1.1 noro 269: int *p;
270: int **pp;
271: ML blist,clist,bqlist,cqlist,rlist;
272: UM *b;
273: LUM fl,tl;
274: LUM *l;
1.2 noro 275: UM w;
1.1 noro 276:
1.2 noro 277: w = W_UMALLOC(UDEG(f));
1.1 noro 278: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
1.2 noro 279:
280: /* c[0] must have lc(f) */
281: blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
282: ptoum(mod,mfl[0],w);
283: inv = invm(w->c[UDEG(mfl[0])],mod);
284: lc = rem(NM((Q)LC(f)),mod);
285: if ( SGN((Q)LC(f)) < 0 )
286: lc = (mod-lc)%mod;
287: lc = dmar(inv,lc,0,mod);
288: if ( lc == 1 )
289: copyum(w,blist->c[0]);
290: else
291: mulsum(mod,w,lc,blist->c[0]);
292:
293: /* c[i] (i=1,...,nf-1) must be monic */
294: for ( i = 1; i < nf; i++ ) {
1.1 noro 295: blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
1.2 noro 296: ptoum(mod,mfl[i],w);
297: inv = invm(w->c[UDEG(mfl[i])],mod);
298: if ( inv == 1 )
299: copyum(w,blist->c[i]);
300: else
301: mulsum(mod,w,inv,blist->c[i]);
1.1 noro 302: }
1.2 noro 303:
1.1 noro 304: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
305: n = bqlist->n; q = bqlist->mod;
306: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
307: if ( bound == 1 ) {
308: *listp = rlist = MLALLOC(n);
309: rlist->n = n; rlist->mod = q; rlist->bound = bound;
310: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
311: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
312: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
313: pp[j][0] = p[j];
314: }
315: } else {
316: W_LUMALLOC((int)UDEG(f),bound,fl);
317: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
318: }
319: }
320:
321: void resf_dtest(f,list,nb,nfl,mdl,dcp)
322: P f;
323: ML list;
324: int nb;
325: int *nfl,*mdl;
326: DCP *dcp;
327: {
328: int n,np,bound,q;
329: int i,j,k;
330: int *win;
331: P g,factor,cofactor;
332: Q csum,csumt;
333: DCP dcf,dcf0;
334: LUM *c;
335: ML wlist;
336: struct resf_dlist *dlist;
337: int z;
338:
339: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
340: win = W_ALLOC(np+1);
341: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
342: wlist = W_MLALLOC(np); wlist->n = list->n;
343: wlist->mod = list->mod; wlist->bound = list->bound;
344: c = (LUM *)COEF(wlist);
345: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
346: dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
347: for ( i = 0, j = 0; j < nb; j++ )
348: for ( k = 0; k < nfl[j]; k++, i++ ) {
349: dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
350: }
351: k = nb;
352: for ( i = 0; i < nb; i++ )
353: win[i] = i+1;
354: for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
355: #if 0
356: if ( !(z++ % 10000) )
357: fputc('.',stderr);
358: #endif
359: if ( resf_degtest(k,win,nb,dlist) &&
360: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
361: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
362: g = cofactor;
363: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
364: for ( i = 0; i < k - 1; i++ )
365: for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
366: c[j-i-1] = c[j];
367: dlist[j-i-1] = dlist[j];
368: }
369: for ( j = win[k-1] + 1; j <= np; j++ ) {
370: c[j-k] = c[j];
371: dlist[j-k] = dlist[j];
372: }
373: if ( ( np -= k ) < k )
374: break;
375: if ( np - win[0] + 1 < k )
376: if ( ++k > np )
377: break;
378: else
379: for ( i = 0; i < k; i++ )
380: win[i] = i + 1;
381: else
382: for ( i = 1; i < k; i++ )
383: win[i] = win[0] + i;
384: } else if ( !ncombi(1,np,k,win) )
385: if ( k == np )
386: break;
387: else
388: for ( i = 0, ++k; i < k; i++ )
389: win[i] = i + 1;
390: }
391: NEXTDC(dcf0,dcf); COEF(dcf) = g;
392: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
393: }
394:
395: void resf_dtest_special(f,list,nb,nfl,mdl,dcp)
396: P f;
397: ML list;
398: int nb;
399: int *nfl,*mdl;
400: DCP *dcp;
401: {
402: int n,np,bound,q;
403: int i,j;
404: int *win;
405: P g,factor,cofactor;
406: Q csum,csumt;
407: DCP dcf,dcf0;
408: LUM *c;
409: ML wlist;
410: int max;
411:
412: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
413: win = W_ALLOC(np+1);
414: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
415: wlist = W_MLALLOC(np); wlist->n = list->n;
416: wlist->mod = list->mod; wlist->bound = list->bound;
417: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
418: max = 1<<nb;
419: for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
420: #if 0
421: if ( !(i % 1000) )
422: fprintf(stderr,"i=%d\n",i);
423: #endif
424: for ( j = 0; j < nb; j++ )
425: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
426: if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
427: #if 0
428: fprintf(stderr,"success : i=%d\n",i);
429: #endif
430: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
431: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
432: NEXT(dcf) = 0;*dcp = dcf0;
433: return;
434: }
435: }
436: NEXTDC(dcf0,dcf); COEF(dcf) = g;
437: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
438: }
439:
440: #if 0
441: void Pftest(arg,rp)
442: NODE arg;
443: P *rp;
444: {
445: ML list;
446: DCP dc;
447: P p;
448: P *mfl;
449:
450: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
451: hensel_special(4,1,p,mfl,&list);
452: dtest_special(p,list,rp);
453: }
454:
455: void dtest_special(f,list,pr)
456: P f;
457: ML list;
458: P *pr;
459: {
460: int n,np,bound,q;
461: int i,j,k;
462: int *win;
463: P g,factor,cofactor;
464: Q csum,csumt;
465: DCP dc,dcf,dcf0;
466: LUM *c;
467: ML wlist;
468:
469: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
470: win = W_ALLOC(np+1);
471: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
472: wlist = W_MLALLOC(np); wlist->n = list->n;
473: wlist->mod = list->mod; wlist->bound = list->bound;
474: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
475: for ( g = f, i = 130000; i < (1<<20); i++ ) {
476: #if 0
477: if ( !(i % 1000) )
478: fprintf(stderr,"i=%d\n",i);
479: #endif
480: for ( j = 0; j < 20; j++ )
481: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
482: if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
483: #if 0
484: fprintf(stderr,"success : i=%d\n",i);
485: #endif
486: *pr = factor; return;
487: }
488: }
489: }
490:
491: void hensel_special(index,count,f,mfl,listp)
492: int index,count;
493: P f;
494: P *mfl;
495: ML *listp;
496: {
497: register int i,j;
498: int q,n,t,d,r,u,br,tmp,bound;
499: int *c,*p,*m,*w;
500: int **pp;
501: DCP dc;
502: ML blist,clist,bqlist,cqlist,rlist;
503: UM *b;
504: LUM fl,tl;
505: LUM *l;
506:
507: blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
508: for ( i = 0; i < 40; i++ ) {
509: blist->c[i] = (pointer)UMALLOC(6);
510: ptoum(11,mfl[i],blist->c[i]);
511: }
512: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
513: n = bqlist->n; q = bqlist->mod;
514: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
515: if ( bound == 1 ) {
516: *listp = rlist = MLALLOC(n);
517: rlist->n = n; rlist->mod = q; rlist->bound = bound;
518: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
519: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
520: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
521: pp[j][0] = p[j];
522: }
523: } else {
524: W_LUMALLOC(UDEG(f),bound,fl);
525: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
526: }
527: }
528: #endif
529:
530: #if 0
531: void Pftest(arg,rp)
532: NODE arg;
533: P *rp;
534: {
535: ML list;
536: DCP dc;
537: P p;
538: P *mfl;
539:
540: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
541: hensel_special(2,1,p,mfl,&list);
542: dtest_special(p,list,rp);
543: }
544:
545: void dtest_special(f,list,pr)
546: P f;
547: ML list;
548: P *pr;
549: {
550: int n,np,bound,q;
551: int i,j,k,t,b0;
552: int *win;
553: P g,factor,cofactor;
554: Q csum,csumt;
555: DCP dc,dcf,dcf0;
556: LUM *c;
557: ML wlist;
558: static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
559:
560: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
561: win = W_ALLOC(np+1);
562: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
563: wlist = W_MLALLOC(np); wlist->n = list->n;
564: wlist->mod = list->mod; wlist->bound = list->bound;
565: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
566: for ( g = f, i = 0; i < (1<<23); i++ ) {
567: #if 0
568: if ( !(i % 1000) )
569: fprintf(stderr,"i=%d\n",i);
570: #endif
571: t = i>>20; b0 = nbits[t];
572: if ( !b0 )
573: continue;
574: for ( j = 1; j < 6; j++ ) {
575: t = (i>>(20-4*j))&0xf;
576: if ( nbits[t] != b0 )
577: break;
578: }
579: if ( j != 6 )
580: continue;
581: for ( j = k = 0; j < 24; j++ )
582: if ( i & (1<<(23-j)) )
583: win[k++] = j;
584: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
585: #if 0
586: fprintf(stderr,"success : i=%d\n",i);
587: #endif
588: *pr = factor; return;
589: }
590: }
591: *pr = f;
592: }
593:
594: void hensel_special(index,count,f,mfl,listp)
595: int index,count;
596: P f;
597: P *mfl;
598: ML *listp;
599: {
600: register int i,j;
601: int q,n,t,d,r,u,br,tmp,bound;
602: int *c,*p,*m,*w;
603: int **pp;
604: DCP dc;
605: ML blist,clist,bqlist,cqlist,rlist;
606: UM *b;
607: LUM fl,tl;
608: LUM *l;
609:
610: blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
611: for ( i = 0; i < 24; i++ ) {
612: blist->c[i] = (pointer)UMALLOC(7);
613: ptoum(5,mfl[i],blist->c[i]);
614: }
615: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
616: n = bqlist->n; q = bqlist->mod;
617: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
618: if ( bound == 1 ) {
619: *listp = rlist = MLALLOC(n);
620: rlist->n = n; rlist->mod = q; rlist->bound = bound;
621: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
622: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
623: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
624: pp[j][0] = p[j];
625: }
626: } else {
627: W_LUMALLOC(UDEG(f),bound,fl);
628: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
629: }
630: }
631: #endif
632:
633: void Pftest(arg,rp)
634: NODE arg;
635: P *rp;
636: {
637: ML list;
638: P p;
639: P *mfl;
640:
641: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
642: hensel_special(5,1,p,mfl,&list);
643: dtest_special(p,list,rp);
644: }
645:
646: int nbits(a)
647: int a;
648: {
649: int i,s;
650:
651: for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
652: if ( a & 1 ) s++;
653: return s;
654: }
655:
656: void dtest_special(f,list,pr)
657: P f;
658: ML list;
659: P *pr;
660: {
661: int n,np,bound,q;
662: int i,j,k,b0;
663: int *win;
664: P g,factor,cofactor;
665: Q csum,csumt;
666: LUM *c;
667: ML wlist;
668:
669: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
670: win = W_ALLOC(np+1);
671: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
672: wlist = W_MLALLOC(np); wlist->n = list->n;
673: wlist->mod = list->mod; wlist->bound = list->bound;
674: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
675: for ( g = f, i = 0; i < (1<<19); i++ ) {
676: #if 0
677: if ( !(i % 10000) )
678: fprintf(stderr,"i=%d\n",i);
679: #endif
680: b0 = nbits(i>>10);
681: if ( !b0 || (nbits(i&0x3ff) != b0) )
682: continue;
683: for ( j = k = 0; j < 20; j++ )
684: if ( i & (1<<(19-j)) )
685: win[k++] = j;
686: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
687: #if 0
688: fprintf(stderr,"success : i=%d\n",i);
689: #endif
690: *pr = factor; return;
691: }
692: }
693: *pr = f;
694: }
695:
696: void hensel_special(index,count,f,mfl,listp)
697: int index,count;
698: P f;
699: P *mfl;
700: ML *listp;
701: {
702: register int i,j;
703: int q,n,bound;
704: int *p;
705: int **pp;
706: ML blist,clist,bqlist,cqlist,rlist;
707: UM *b;
708: LUM fl,tl;
709: LUM *l;
710:
711: blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
712: for ( i = 0; i < 20; i++ ) {
713: blist->c[i] = (pointer)UMALLOC(10);
714: ptoum(11,mfl[i],blist->c[i]);
715: }
716: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
717: n = bqlist->n; q = bqlist->mod;
718: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
719: if ( bound == 1 ) {
720: *listp = rlist = MLALLOC(n);
721: rlist->n = n; rlist->mod = q; rlist->bound = bound;
722: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
723: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
724: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
725: pp[j][0] = p[j];
726: }
727: } else {
728: W_LUMALLOC((int)UDEG(f),bound,fl);
729: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
730: }
731: }
732:
733: void Pnullspace(arg,rp)
734: NODE arg;
735: LIST *rp;
736: {
737: int i,j,n,mod;
738: MAT mat,r;
739: VECT u;
740: V v;
741: P p,z;
742: Q q;
743: UM **w;
744: UM mp;
745: P *t;
746: UM *s;
747: int *ind;
748: NODE n0,n1;
749:
750: mat = (MAT)ARG0(arg);
751: p = (P)ARG1(arg);
752: v = VR(p);
753: mod = QTOS((Q)ARG2(arg));
754: n = mat->row;
755: w = (UM **)almat_pointer(n,n);
756: for ( i = 0; i < n; i++ )
757: for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
758: ptomp(mod,t[j],&z);
759: s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
760: mptoum(z,s[j]);
761: }
762: mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
763: ind = (int *)ALLOCA(n*sizeof(int));
764: nullspace(w,mp,mod,n,ind);
765: MKMAT(r,n,n);
766: for ( i = 0; i < n; i++ )
767: for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
768: umtop(v,s[j],&t[j]);
769: MKVECT(u,n);
770: for ( i = 0; i < n; i++ ) {
771: STOQ(ind[i],q); u->body[i] = (pointer)q;
772: }
773: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
774: }
775:
776: void nullspace(mat,p,mod,n,ind)
777: UM **mat;
778: UM p;
779: int mod,n;
780: int *ind;
781: {
782: int i,j,l,s,d;
783: UM q,w,w1,h,inv;
784: UM *t,*u;
785:
786: d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
787: w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
788: bzero(ind,n*sizeof(int));
789: ind[0] = 0;
790: for ( i = j = 0; j < n; i++, j++ ) {
791: for ( ; j < n; j++ ) {
792: for ( l = i; l < n; l++ )
793: if ( DEG(mat[l][j])>=0 )
794: break;
795: if ( l < n ) {
796: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
797: } else
798: ind[j] = 1;
799: }
800: if ( j == n )
801: break;
802: invum(mod,p,mat[i][j],inv);
803: for ( s = j, t = mat[i]; s < n; s++ ) {
804: mulum(mod,t[s],inv,w);
805: DEG(w) = divum(mod,w,p,q);
806: cpyum(w,t[s]);
807: }
808: for ( l = 0; l < n; l++ ) {
809: if ( l == i )
810: continue;
811: u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
812: for ( s = j; s < n; s++ ) {
813: mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
814: DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
815: }
816: }
817: }
818: }
819:
820: void Pnullspace_ff(arg,rp)
821: NODE arg;
822: LIST *rp;
823: {
824: int i,j,n;
825: Q mod;
826: MAT mat,r;
827: VECT u;
828: Q q;
829: Obj **w;
830: Obj *t;
831: Obj *s;
832: int *ind;
833: NODE n0,n1;
834:
835: mat = (MAT)ARG0(arg);
836: n = mat->row;
837: w = (Obj **)almat_pointer(n,n);
838: for ( i = 0; i < n; i++ )
839: for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
840: s[j] = t[j];
841: ind = (int *)ALLOCA(n*sizeof(int));
842: switch ( current_ff ) {
843: case FF_GFP:
844: nullspace_lm((LM **)w,n,ind); break;
845: case FF_GF2N:
846: nullspace_gf2n((GF2N **)w,n,ind); break;
847: case FF_GFPN:
848: nullspace_gfpn((GFPN **)w,n,ind); break;
1.5 noro 849: case FF_GFS:
850: nullspace_gfs((GFS **)w,n,ind); break;
1.1 noro 851: default:
852: error("nullspace_ff : current_ff is not set");
853: }
854: MKMAT(r,n,n);
855: for ( i = 0; i < n; i++ )
856: for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
857: t[j] = s[j];
858: MKVECT(u,n);
859: for ( i = 0; i < n; i++ ) {
860: STOQ(ind[i],q); u->body[i] = (pointer)q;
861: }
862: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
863: }
864:
865: void nullspace_lm(mat,n,ind)
866: LM **mat;
867: int n;
868: int *ind;
869: {
870: int i,j,l,s;
871: Q q,mod;
872: N lm;
873: LM w,w1,h,inv;
874: LM *t,*u;
875:
876: getmod_lm(&lm); NTOQ(lm,1,mod);
877:
878: bzero(ind,n*sizeof(int));
879: ind[0] = 0;
880: for ( i = j = 0; j < n; i++, j++ ) {
881: for ( ; j < n; j++ ) {
882: for ( l = i; l < n; l++ ) {
883: simplm(mat[l][j],&w); mat[l][j] = w;
884: if ( mat[l][j] )
885: break;
886: }
887: if ( l < n ) {
888: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
889: } else
890: ind[j] = 1;
891: }
892: if ( j == n )
893: break;
894: NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);
895: for ( s = j, t = mat[i]; s < n; s++ ) {
896: mullm(t[s],inv,&w); t[s] = w;
897: }
898: for ( l = 0; l < n; l++ ) {
899: if ( l == i )
900: continue;
901: u = mat[l]; chsgnlm(u[j],&h);
902: for ( s = j; s < n; s++ ) {
903: mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
904: }
905: }
906: }
907: }
908:
909: void nullspace_gf2n(mat,n,ind)
910: GF2N **mat;
911: int n;
912: int *ind;
913: {
914: int i,j,l,s;
915: GF2N w,w1,h,inv;
916: GF2N *t,*u;
917: extern gf2n_lazy;
918:
919: bzero(ind,n*sizeof(int));
920: ind[0] = 0;
921: for ( i = j = 0; j < n; i++, j++ ) {
922: for ( ; j < n; j++ ) {
923: for ( l = i; l < n; l++ ) {
924: simpgf2n(mat[l][j],&w); mat[l][j] = w;
925: if ( mat[l][j] )
926: break;
927: }
928: if ( l < n ) {
929: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
930: } else
931: ind[j] = 1;
932: }
933: if ( j == n )
934: break;
935: invgf2n(mat[i][j],&inv);
936: for ( s = j, t = mat[i]; s < n; s++ ) {
937: mulgf2n(t[s],inv,&w); t[s] = w;
938: }
939: for ( l = 0; l < n; l++ ) {
940: if ( l == i )
941: continue;
942: u = mat[l]; h = u[j];
943: for ( s = j; s < n; s++ ) {
944: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
945: }
946: }
947: }
948: }
949:
950: void nullspace_gfpn(mat,n,ind)
951: GFPN **mat;
952: int n;
953: int *ind;
954: {
955: int i,j,l,s;
956: GFPN w,w1,h,inv;
957: GFPN *t,*u;
958:
959: bzero(ind,n*sizeof(int));
960: ind[0] = 0;
961: for ( i = j = 0; j < n; i++, j++ ) {
962: for ( ; j < n; j++ ) {
963: for ( l = i; l < n; l++ ) {
964: simpgfpn(mat[l][j],&w); mat[l][j] = w;
965: if ( mat[l][j] )
966: break;
967: }
968: if ( l < n ) {
969: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
970: } else
971: ind[j] = 1;
972: }
973: if ( j == n )
974: break;
975: divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
976: for ( s = j, t = mat[i]; s < n; s++ ) {
977: mulgfpn(t[s],inv,&w); t[s] = w;
978: }
979: for ( l = 0; l < n; l++ ) {
980: if ( l == i )
981: continue;
982: u = mat[l]; chsgngfpn(u[j],&h);
983: for ( s = j; s < n; s++ ) {
984: mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
985: }
986: }
987: }
988: }
1.5 noro 989:
990: void nullspace_gfs(mat,n,ind)
991: GFS **mat;
992: int n;
993: int *ind;
994: {
995: int i,j,l,s;
996: GFS w,w1,h,inv;
997: GFS *t,*u;
998: GFS one;
999:
1000: bzero(ind,n*sizeof(int));
1001: ind[0] = 0;
1002: mqtogfs(ONEM,&one);
1003:
1004: for ( i = j = 0; j < n; i++, j++ ) {
1005: for ( ; j < n; j++ ) {
1006: for ( l = i; l < n; l++ )
1007: if ( mat[l][j] )
1008: break;
1009: if ( l < n ) {
1010: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1011: } else
1012: ind[j] = 1;
1013: }
1014: if ( j == n )
1015: break;
1016: divgfs(one,mat[i][j],&inv);
1017: for ( s = j, t = mat[i]; s < n; s++ ) {
1018: mulgfs(t[s],inv,&w); t[s] = w;
1019: }
1020: for ( l = 0; l < n; l++ ) {
1021: if ( l == i )
1022: continue;
1023: u = mat[l];
1024: chsgngfs(u[j],&h);
1025: for ( s = j; s < n; s++ ) {
1026: mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
1027: }
1028: }
1029: }
1030: }
1031:
1.1 noro 1032: /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
1033:
1034: void linear_form_to_array(p,vl,m,array)
1035: P p;
1036: VL vl;
1037: int m;
1038: Num *array;
1039: {
1040: int i;
1041: DCP dc;
1042:
1043: bzero((char *)array,(m+1)*sizeof(Num *));
1044: for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
1045: if ( ID(p) == O_N )
1046: break;
1047: else if ( VR(p) == vl->v ) {
1048: dc = DC(p);
1049: array[i] = (Num)COEF(dc);
1050: dc = NEXT(dc);
1051: p = dc ? COEF(dc) : 0;
1052: }
1053: }
1054: array[m] = (Num)p;
1055: }
1056:
1057: void array_to_linear_form(array,vl,m,r)
1058: Num *array;
1059: VL vl;
1060: int m;
1061: P *r;
1062: {
1063: P t;
1064: DCP dc0,dc1;
1065:
1066: if ( !m )
1067: *r = (P)array[0];
1068: else {
1069: array_to_linear_form(array+1,NEXT(vl),m-1,&t);
1070: if ( !array[0] )
1071: *r = t;
1072: else {
1073: NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
1074: if ( !t )
1075: NEXT(dc0) = 0;
1076: else {
1077: NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
1078: NEXT(dc1) = 0;
1079: NEXT(dc0) = dc1;
1080: }
1081: MKP(vl->v,dc0,*r);
1082: }
1083: }
1084: }
1085:
1086: void Psolve_linear_equation_gf2n(arg,rp)
1087: NODE arg;
1088: LIST *rp;
1089: {
1090: NODE eqs,tn;
1091: VL vars,tvl;
1092: int i,j,n,m,dim,codim;
1093: GF2N **w;
1094: int *ind;
1095: NODE n0,n1;
1096:
1097: get_vars(ARG0(arg),&vars);
1098: eqs = BDY((LIST)ARG0(arg));
1099: for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
1100: for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
1101: w = (GF2N **)almat_pointer(n,m+1);
1102: for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
1103: linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
1104: ind = (int *)ALLOCA(m*sizeof(int));
1105: solve_linear_equation_gf2n(w,n,m,ind);
1106: for ( j = 0, dim = 0; j < m; j++ )
1107: if ( ind[j] )
1108: dim++;
1109: codim = m-dim;
1110: for ( i = codim; i < n; i++ )
1111: if ( w[i][m] ) {
1112: MKLIST(*rp,0); return;
1113: }
1114: for ( i = 0, n0 = 0; i < codim; i++ ) {
1115: NEXTNODE(n0,n1);
1116: array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
1117: }
1118: if ( n0 )
1119: NEXT(n1) = 0;
1120: MKLIST(*rp,n0);
1121: }
1122:
1123: void solve_linear_equation_gf2n(mat,n,m,ind)
1124: GF2N **mat;
1125: int n;
1126: int *ind;
1127: {
1128: int i,j,l,s;
1129: GF2N w,w1,h,inv;
1130: GF2N *t,*u;
1131: extern gf2n_lazy;
1132:
1133: bzero(ind,m*sizeof(int));
1134: ind[0] = 0;
1135: for ( i = j = 0; j < m; i++, j++ ) {
1136: for ( ; j < m; j++ ) {
1137: for ( l = i; l < n; l++ ) {
1138: simpgf2n(mat[l][j],&w); mat[l][j] = w;
1139: if ( mat[l][j] )
1140: break;
1141: }
1142: if ( l < n ) {
1143: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1144: } else
1145: ind[j] = 1;
1146: }
1147: if ( j == m )
1148: break;
1149: invgf2n(mat[i][j],&inv);
1150: for ( s = j, t = mat[i]; s <= m; s++ ) {
1151: mulgf2n(t[s],inv,&w); t[s] = w;
1152: }
1153: for ( l = 0; l < n; l++ ) {
1154: if ( l == i )
1155: continue;
1156: u = mat[l]; h = u[j];
1157: for ( s = j; s <= m; s++ ) {
1158: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
1159: }
1160: }
1161: }
1162: }
1163:
1164: /*
1165: void null_to_sol(mat,ind,mod,n,r)
1166: int **mat;
1167: int *ind;
1168: int mod,n;
1169: UM *r;
1170: {
1171: int i,j,k,l;
1172: int *c;
1173: UM w;
1174:
1175: for ( i = 0, l = 0; i < n; i++ ) {
1176: if ( !ind[i] )
1177: continue;
1178: w = UMALLOC(n);
1179: for ( j = k = 0, c = COEF(w); j < n; j++ )
1180: if ( ind[j] )
1181: c[j] = 0;
1182: else
1183: c[j] = mat[k++][i];
1184: c[i] = mod-1;
1185: for ( j = n; j >= 0; j-- )
1186: if ( c[j] )
1187: break;
1188: DEG(w) = j;
1189: r[l++] = w;
1190: }
1191: }
1192: */
1193:
1194: void showgfmat(mat,n)
1195: UM **mat;
1196: int n;
1197: {
1198: int i,j,k;
1199: int *c;
1200: UM p;
1201:
1202: for ( i = 0; i < n; i++ ) {
1203: for ( j = 0; j < n; j++ ) {
1204: p = mat[i][j];
1205: if ( DEG(p) < 0 )
1206: fprintf(asir_out,"0");
1207: else
1208: for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
1209: if ( c[k] )
1210: fprintf(asir_out,"+%d",c[k]);
1211: if ( k > 1 )
1212: fprintf(asir_out,"a^%d",k);
1213: else if ( k == 1 )
1214: fprintf(asir_out,"a",k);
1215: }
1216: fprintf(asir_out," ");
1217: }
1218: fprintf(asir_out,"\n");
1219: }
1220: }
1221:
1222: #if 0
1223: void Pgcda_mod(arg,rp)
1224: NODE arg;
1225: P *rp;
1226: {
1227: p1 = (P)ARG0(arg);
1228: p2 = (P)ARG1(arg);
1229: v = VR((P)ARG2(arg));
1230: d = (P)ARG3(arg);
1231: m = QTOS((Q)ARG4(arg));
1232: reordvar(CO,v,&vl);
1233: reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
1234: reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
1235: if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
1236: *rp = ONE; return;
1237: }
1238: if ( deg(v,m1) >= deg(v,m2) ) {
1239: t = m1; m1 = m2; m2 = t;
1240: }
1241: while ( 1 ) {
1242: inva_mod(COEF(DC(m2)),d,m,&inv);
1243: NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
1244: d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
1245: mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
1246: mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
1247: }
1248: }
1249: #endif
1250:
1251: void Ppwr_mod(arg,rp)
1252: NODE arg;
1253: P *rp;
1254: {
1255: P p,a,d,r;
1256: int m;
1257: Q q;
1258: N n;
1259:
1260: m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;
1261: ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
1262: ptomp(m,(P)ARG3(arg),&d);
1263: pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
1264: mptop(r,rp);
1265: }
1266:
1267: void pwr_mod(p,a,v,d,m,n,rp)
1268: P p,a,d;
1269: V v;
1270: int m;
1271: N n;
1272: P *rp;
1273: {
1274: int b;
1275: P t,s,r;
1276: N n1;
1277:
1278: if ( !n )
1279: *rp = (P)ONEM;
1280: else if ( UNIN(n) )
1281: *rp = p;
1282: else {
1283: b = divin(n,2,&n1);
1284: pwr_mod(p,a,v,d,m,n1,&t);
1285: mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
1286: if ( b ) {
1287: mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
1288: } else
1289: *rp = r;
1290: }
1291: }
1292:
1293: void rem_mod(p,a,v,d,m,rp)
1294: P p,a,d;
1295: V v;
1296: int m;
1297: P *rp;
1298: {
1299: P tmp,r1,r2;
1300:
1301: divsrmp(CO,m,p,d,&tmp,&r1);
1302: divsrmp(CO,m,r1,a,&tmp,&r2);
1303: divsrmp(CO,m,r2,d,&tmp,rp);
1304: }