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