![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to D.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to D.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/D.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void dtest(f,list,hint,dcp)
5: P f;
6: ML list;
7: int hint;
8: DCP *dcp;
9: {
10: int n,np,bound,q;
11: int i,j,k;
12: int *win;
13: P g,factor,cofactor;
14: Q csum,csumt;
15: DCP dcf,dcf0;
16: LUM *c;
17: ML wlist;
18: int z;
19:
20: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
21: win = W_ALLOC(np+1);
22: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
23: wlist = W_MLALLOC(np); wlist->n = list->n;
24: wlist->mod = list->mod; wlist->bound = list->bound;
25: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
26: for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; ) {
27: #if 0
28: if ( !(++z % 10000) )
29: fprintf(stderr,"z=%d\n",z);
30: #endif
31: if ( degtest(k,win,wlist,hint) &&
32: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
33: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
34: g = cofactor;
35: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
36: for ( i = 0; i < k - 1; i++ )
37: for ( j = win[i] + 1; j < win[i + 1]; j++ )
38: c[j-i-1] = c[j];
39: for ( j = win[k-1] + 1; j <= np; j++ )
40: c[j-k] = c[j];
41: if ( ( np -= k ) < k )
42: break;
43: if ( np - win[0] + 1 < k )
44: if ( ++k > np )
45: break;
46: else
47: for ( i = 0; i < k; i++ )
48: win[i] = i + 1;
49: else
50: for ( i = 1; i < k; i++ )
51: win[i] = win[0] + i;
52: } else if ( !ncombi(1,np,k,win) )
53: if ( k == np )
54: break;
55: else
56: for ( i = 0, ++k; i < k; i++ )
57: win[i] = i + 1;
58: }
59: NEXTDC(dcf0,dcf); COEF(dcf) = g;
60: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
61: }
62:
63: void dtestsql(f,list,dc,dcp)
64: P f;
65: ML list;
66: struct oDUM *dc;
67: DCP *dcp;
68: {
69: int j,n,m,b;
70: P t,s,fq,fr;
71: P *true;
72: Q tq;
73: LUM *c;
74: DCP dcr,dcr0;
75:
76: n = list->n; m = list->mod; b = list->bound; c = (LUM *)list->c;
77: true = (P*)ALLOCA(n*sizeof(P));
78: for ( j = 0; j < n; j++ ) {
79: dtestsq(m,b,f,c[j],&t);
80: if ( t )
81: true[j] = t;
82: else {
83: *dcp = 0;
84: return;
85: }
86: }
87: for ( t = f, j = 0; j < n; j++ ) {
88: STOQ(dc[j].n,tq); pwrp(CO,true[j],tq,&s); udivpz(t,s,&fq,&fr);
89: if ( fq && !fr )
90: t = fq;
91: else {
92: *dcp = 0;
93: return;
94: }
95: }
96: for ( j = 0, dcr = dcr0 = 0; j < n; j++ ) {
97: NEXTDC(dcr0,dcr); STOQ(dc[j].n,DEG(dcr)); COEF(dcr) = true[j];
98: }
99: NEXT(dcr) = 0; *dcp = dcr0;
100: }
101:
102: void dtestsq(q,bound,f,fl,g)
103: int q,bound;
104: P f;
105: LUM fl;
106: P *g;
107: {
108: P lcf,t,fq,fr,s;
109: struct oML list;
110: int in = 0;
111:
112: list.n = 1;
113: list.mod = q;
114: list.bound = bound;
115: list.c[0] = (pointer)fl;
116:
117: mullumarray(f,&list,1,&in,&t); mulp(CO,f,COEF(DC(f)),&lcf);
118: udivpz(lcf,t,&fq,&fr);
119: if( fq && !fr )
120: ptozp(t,1,(Q *)&s,g);
121: else
122: *g = 0;
123: }
124:
125: void dtestroot(m,b,f,fl,dc,dcp)
126: int m,b;
127: P f;
128: LUM fl;
129: struct oDUM *dc;
130: DCP *dcp;
131: {
132: P t,s,u;
133: DCP dcr;
134: Q q;
135:
136: dtestroot1(m,b,f,fl,&t);
137: if ( !t ) {
138: *dcp = 0;
139: return;
140: }
141: STOQ(dc[0].n,q); pwrp(CO,t,q,&s); subp(CO,s,f,&u);
142: if ( u )
143: *dcp = 0;
144: else {
145: NEWDC(dcr); STOQ(dc[0].n,DEG(dcr));
146: COEF(dcr) = t; NEXT(dcr) = 0; *dcp = dcr;
147: }
148: }
149:
150: void dtestroot1(q,bound,f,fl,g)
151: int q,bound;
152: P f;
153: LUM fl;
154: P *g;
155: {
156: P fq,fr,t;
157:
158: lumtop(VR(f),q,bound,fl,&t); udivpz(f,t,&fq,&fr);
159: *g = (fq && !fr) ? t : 0;
160: }
161:
162: int dtestmain(g,csum,list,k,in,fp,cofp)
163: P g;
164: Q csum;
165: ML list;
166: int k;
167: int *in;
168: P *fp,*cofp;
169: {
170: int mod;
171: P fmul,lcg;
172: Q csumg;
173: N nq,nr;
174: P fq,fr,t;
175:
176: if (!ctest(g,list,k,in))
177: return 0;
178: mod = list->mod;
179: mullumarray(g,list,k,in,&fmul); mulp(CO,g,COEF(DC(g)),&lcg);
180: if ( csum ) {
181: ucsump(fmul,&csumg);
182: if ( csumg ) {
183: divn(NM(csum),NM(csumg),&nq,&nr);
184: if ( nr )
185: return 0;
186: }
187: }
188: udivpz(lcg,fmul,&fq,&fr);
189: if ( fq && !fr ) {
190: ptozp(fq,1,(Q *)&t,cofp); ptozp(fmul,1,(Q *)&t,fp);
191: return 1;
192: } else
193: return 0;
194: }
195:
196: int degtest(k,win,list,hint)
197: int k;
198: int *win;
199: ML list;
200: int hint;
201: {
202: register int i,d;
203: LUM *c;
204:
205: if ( hint == 1 )
206: return 1;
207: for ( c = (LUM*)list->c, i = 0, d = 0; i < k; i++ )
208: d += DEG(c[win[i]]);
209: return !(d % hint);
210: }
211:
212: int ctest(g,list,k,in)
213: P g;
214: ML list;
215: int k;
216: int *in;
217: {
218: register int i;
219: int q,bound;
220: int *wm,*wm1,*tmpp;
221: DCP dc;
222: Q dvr;
223: N lcn,cstn,dndn,dmyn,rn;
224: LUM *l;
225:
226: for ( dc = DC(g); dc && DEG(dc); dc = NEXT(dc) );
227: if ( dc )
228: cstn = NM((Q)COEF(dc));
229: else
230: return 1;
231: q = list->mod; bound = list->bound;
232: ntobn(q,NM((Q)COEF(DC(g))),&lcn);;
233: W_CALLOC(bound+1,int,wm); W_CALLOC(bound+1,int,wm1);
234: for ( i = 0; i < PL(lcn); i++ )
235: wm[i] = BD(lcn)[i];
236: for ( i = 0, l = (LUM *)list->c; i < k; i++ ) {
237: mulpadic(q,bound,wm,COEF(l[in[i]])[0],wm1);
238: tmpp = wm; wm = wm1; wm1 = tmpp;
239: }
240: padictoq(q,bound,wm,&dvr);
241: kmuln(NM((Q)COEF(DC(g))),cstn,&dndn); divn(dndn,NM(dvr),&dmyn,&rn);
242: return rn ? 0 : 1;
243: }
244:
245: /*
246: int ncombi(n0,n,k,c)
247: int n0,n,k,*c;
248: {
249: register int i,tmp;
250:
251: if ( !k )
252: return 0;
253: if ( !ncombi(c[1],n,k-1,c+1) ) {
254: if ( c[0] + k > n )
255: return 0;
256: else {
257: for ( i = 0, tmp = c[0]; i < k; i++ )
258: c[i] = tmp + i + 1;
259: return 1;
260: }
261: } else
262: return 1;
263: }
264: */
265:
266: int ncombi(n0,n,k,c)
267: int n0,n,k,*c;
268: {
269: register int i,t;
270:
271: if ( !k )
272: return 0;
273: for ( i = k-1; i >= 0 && c[i] == n+i-(k-1); i-- );
274: if ( i < 0 )
275: return 0;
276: t = ++c[i++];
277: for ( t++ ; i < k; i++, t++ )
278: c[i] = t;
279: return 1;
280: }
281:
282: void nthrootn(number,n,root)
283: int n;
284: N number,*root;
285: {
286: N s,t,u,pn,base,n1,n2,q,r,gcd,num;
287: int sgn,index,p,i,tmp,tp,mlr,num0;
288:
289: for ( i = 0; !(n % 2); n /= 2, i++ );
290: for ( index = 0, num = number; ; index++ ) {
291: if ( n == 1 )
292: goto TAIL;
293: p = lprime[index];
294: if ( !p )
295: error("nthrootn : lprime[] exhausted.");
296: if ( !(num0 = rem(num,p)) )
297: continue;
298: STON(n,n1); STON(p-1,n2); gcdn(n1,n2,&gcd);
299: if ( !UNIN(gcd) )
300: continue;
301: tp = pwrm(p,num0,invm(n,p-1)); STON(tp,s);
302: mlr = invm(dmb(p,n,pwrm(p,tp,n-1),&tmp),p);
303: STON(p,base); STON(p,pn);
304: while ( 1 ) {
305: pwrn(s,n,&t); sgn = subn(num,t,&u);
306: if ( !u ) {
307: num = s;
308: break;
309: }
310: if ( sgn < 0 ) {
311: *root = 0;
312: return;
313: }
314: divn(u,base,&q,&r);
315: if ( r ) {
316: *root = 0;
317: return;
318: }
319: STON(dmb(p,mlr,rem(q,p),&tmp),t);
320: kmuln(t,base,&u); addn(u,s,&t); s = t;
321: kmuln(base,pn,&t); base = t;
322: }
323: TAIL :
324: for ( ; i; i-- ) {
325: sqrtn(num,&t);
326: if ( !t ) {
327: *root = 0;
328: return;
329: }
330: num = t;
331: }
332: *root = num;
333: return;
334: }
335: }
336:
337: void sqrtn(number,root)
338: N number,*root;
339: {
340: N a,s,r,q;
341: int sgn;
342:
343: for ( a = ONEN; ; ) {
344: divn(number,a,&q,&r); sgn = subn(q,a,&s);
345: if ( !s ) {
346: *root = !r ? a : 0;
347: return;
348: } else if ( UNIN(s) ) {
349: *root = 0;
350: return;
351: } else {
352: divin(s,2,&q);
353: if ( sgn > 0 )
354: addn(a,q,&r);
355: else
356: subn(a,q,&r);
357: a = r;
358: }
359: }
360: }
361:
362: void lumtop(v,mod,bound,f,g)
363: V v;
364: int mod;
365: int bound;
366: LUM f;
367: P *g;
368: {
369: DCP dc,dc0;
370: int **l;
371: int i;
372: Q q;
373:
374: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
375: padictoq(mod,bound,l[i],&q);
376: if ( q ) {
377: NEXTDC(dc0,dc);
378: if ( i )
379: STOQ(i,DEG(dc));
380: else
381: DEG(dc) = 0;
382: COEF(dc) = (P)q;
383: }
384: }
385: if ( !dc0 )
386: *g = 0;
387: else {
388: NEXT(dc) = 0; MKP(v,dc0,*g);
389: }
390: }
391:
392: void padictoq(mod,bound,p,qp)
393: int mod,bound,*p;
394: Q *qp;
395: {
396: register int h,i,t;
397: int br,sgn;
398: unsigned int *ptr;
399: N n,tn;
400: int *c;
401:
402: c = W_ALLOC(bound);
403: for ( i = 0; i < bound; i++ )
404: c[i] = p[i];
405: h = (mod%2?(mod-1)/2:mod/2); i = bound - 1;
406: while ( i >= 0 && c[i] == h ) i--;
407: if ( i == -1 || c[i] > h ) {
408: for (i = 0, br = 0; i < bound; i++ )
409: if ( ( t = -(c[i] + br) ) < 0 ) {
410: c[i] = t + mod; br = 1;
411: } else {
412: c[i] = 0; br = 0;
413: }
414: sgn = -1;
415: } else
416: sgn = 1;
417: for ( i = bound - 1; ( i >= 0 ) && ( c[i] == 0 ); i--);
418: if ( i == -1 )
419: *qp = 0;
420: else {
421: n = NALLOC(i+1); PL(n) = i+1;
422: for ( i = 0, ptr = BD(n); i < PL(n); i++ )
423: ptr[i] = c[i];
424: bnton(mod,n,&tn); NTOQ(tn,sgn,*qp);
425: }
426: }
427:
428: void mullumarray(f,list,k,in,g)
429: P f;
430: ML list;
431: int k;
432: int *in;
433: P *g;
434: {
435: int np,bound,q,n,i,u;
436: int *tmpp;
437: LUM lclum,wb0,wb1,tlum;
438: LUM *l;
439: N lc;
440:
441: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
442: W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
443: W_LUMALLOC(0,bound,lclum);
444: ntobn(q,NM((Q)COEF(DC(f))),&lc);
445: for ( i = 0, tmpp = COEF(lclum)[0], u = MIN(bound,PL(lc));
446: i < u; i++ )
447: tmpp[i] = BD(lc)[i];
448: l = (LUM *)list->c;
449: mullum(q,bound,lclum,l[in[0]],wb0);
450: for ( i = 1; i < k; i++ ) {
451: mullum(q,bound,l[in[i]],wb0,wb1);
452: tlum = wb0; wb0 = wb1; wb1 = tlum;
453: }
454: lumtop(VR(f),q,bound,wb0,g);
455: }
456:
457: void ucsump(f,s)
458: P f;
459: Q *s;
460: {
461: Q t,u;
462: DCP dc;
463:
464: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
465: addq((Q)COEF(dc),t,&u); t = u;
466: }
467: *s = t;
468: }
469: