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