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