![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to A.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 A.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/A.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void pdiva(vl,p0,p1,p2,pr)
5: VL vl;
6: P p1,p2,p0;
7: P *pr;
8: {
9: Q m,d;
10: P t,q,r;
11:
12: if ( VR(p2) == VR(p0) ) {
13: *pr = p1;
14: return;
15: }
16: if ( VR(p1) == VR(p0) ) {
17: error("pdiva : ???");
18: return;
19: }
20: subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d);
21: pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t);
22: divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr);
23: }
24:
25: void rema(vl,p0,p1,p2,pr)
26: VL vl;
27: P p1,p2,p0,*pr;
28: {
29: P m,m1,m2,g,lc,x,t;
30: Q q;
31: V v;
32: int d1,d2;
33:
34: if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) )
35: *pr = 0;
36: else if ( NUM(p1) )
37: *pr = p1;
38: else if ( ( v = VR(p1) ) == VR(p0) )
39: remsdcp(vl,p1,p0,pr);
40: else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) )
41: *pr = p1;
42: else if ( d1 == d2 ) {
43: mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1);
44: subp(vl,m,m1,&m2);
45: remsdcp(vl,m2,p0,pr);
46: } else {
47: g = p1; lc = LC(p2); MKV(v,x);
48: while ( ( d1 = deg(v,g) ) >= d2 ) {
49: mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1);
50: STOQ(d1-d2,q); pwrp(vl,x,q,&t);
51: mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1);
52: remsdcp(vl,m1,p0,&g);
53: }
54: *pr = g;
55: }
56: }
57:
58: void gcda(vl,p0,p1,p2,pr)
59: VL vl;
60: P p1,p2,p0,*pr;
61: {
62: P g1,g2,r,t,c;
63:
64: remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2);
65: while ( 1 ) {
66: if ( NUM(g2) || (VR(g2) == VR(p0)) ) {
67: *pr = (P)ONE;
68: return;
69: }
70: pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r);
71: if ( !r ) {
72: pmonic(vl,p0,g2,pr);
73: return;
74: }
75: g1 = g2; g2 = r;
76: }
77: }
78:
79: void sqa(vl,p0,p,dcp)
80: VL vl;
81: P p,p0;
82: DCP *dcp;
83: {
84: V v;
85: int i,j,n;
86: P f,f1,f2,gcd,flat,q,r,t;
87: P *a;
88: int *b;
89: DCP dc,dc0;
90:
91: if ( ( v = VR(p) ) == VR(p0) ) {
92: NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE;
93: NEXT(dc) = 0;
94: } else {
95: n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b);
96: for ( i = 0, f = p; ;i++ ) {
97: if ( VR(f) != v )
98: break;
99: remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2);
100: while ( 1 ) {
101: pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r);
102: if ( !r ) {
103: pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r);
104: break;
105: }
106: f1 = f2; f2 = r;
107: }
108: if ( VR(gcd) != v ) {
109: a[i] = f; b[i] = 1;
110: break;
111: }
112: for ( j = 1, f = gcd; ; j++ ) {
113: pqra(vl,p0,f,flat,&q,&r);
114: if ( r )
115: break;
116: else
117: f = q;
118: }
119: a[i] = flat; b[i] = j;
120: }
121: for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) {
122: NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc));
123: if ( a[i+1] )
124: pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t);
125: else
126: COEF(dc) = a[i];
127: }
128: NEXT(dc) = 0; *dcp = dc0;
129: }
130: }
131:
132: void pqra(vl, p0, p1, p2, p, pr)
133: VL vl;
134: P p0,p1,p2;
135: P *p,*pr;
136: {
137: V v,v1;
138: P m,m1,q,r;
139: Q tq;
140: int n;
141:
142: if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) {
143: while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) )
144: vl = NEXT(vl);
145: if ( VR(vl) == v ) {
146: *p = p1; *pr = 0;
147: } else {
148: *p = 0; *pr = p1;
149: }
150: } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) {
151: *pr = p1; *p = 0;
152: } else {
153: n++; STOQ(n,tq);
154: pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r);
155: remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr);
156: }
157: }
158:
159: void pmonic(vl, p0, p, pr)
160: VL vl;
161: P p,p0,*pr;
162: {
163: P d,m,pp,c,t;
164:
165: if ( NUM(p) || ( VR(p) == VR(p0) ) )
166: *pr = (P)ONE;
167: else if ( NUM(COEF(DC(p))) )
168: *pr = p;
169: else {
170: ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d);
171: mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr);
172: }
173: }
174:
175: void remsdcp(vl, p, p0, pr)
176: VL vl;
177: P p,p0,*pr;
178: {
179: P q,r,c;
180: DCP dc,dcr,dcr0;
181:
182: if ( !p )
183: *pr = 0;
184: else if ( NUM(p) )
185: *pr = (P)ONE;
186: else if ( VR(p) == VR(p0) ) {
187: divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr);
188: } else {
189: for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
190: divsrp(vl,COEF(dc),p0,&q,&r);
191: if ( r ) {
192: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r;
193: }
194: }
195: if ( !dcr0 )
196: *pr = 0;
197: else {
198: NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr);
199: }
200: }
201: }
202:
203: void pinva(p0, p, pr)
204: P p,p0,*pr;
205: {
206: V v;
207: Q f,q,q1;
208: Q u,t,a,b,s;
209: P c,c1;
210: P tg,th,tr;
211: P z,zz,zzz;
212: UM wg,wh,ws,wt,wm;
213: int n,m,modres,mod,index,i;
214:
215: v = VR(p); n=deg(v,p); m=deg(v,p0);
216: norm(p,&a); norm(p0,&b);
217: STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s);
218: mulq(t,s,&u);
219: factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f);
220: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
221: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
222: wm = W_UMALLOC(m+n);
223: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
224: mod = lprime[index++];
225: if ( !mod )
226: error("pinva : lprime[] exhausted.");
227: if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) )
228: continue;
229: ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
230: if ( !tr )
231: continue;
232: else
233: modres = CONT((MQ)tr);
234: mptoum(tg,wg); mptoum(th,wh);
235: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
236: DEG(wm) = 0; COEF(wm)[0] = modres;
237: mulum(mod,ws,wm,wt);
238: for ( i = DEG(wt); i >= 0; i-- )
239: if ( ( COEF(wt)[i] * 2 ) > mod )
240: COEF(wt)[i] -= mod;
241: chnrem(mod,v,c,q,wt,&c1,&q1);
242: if ( !ucmpp(c,c1) ) {
243: mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
244: if ( NUM(zzz) ) {
245: q = q1; c = c1; break;
246: }
247: }
248: q = q1; c = c1;
249: if ( cmpq(f,q) < 0 )
250: break;
251: }
252: ptozp(c,1,&s,pr);
253: }