![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dalg.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 dalg.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /*
1.4 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.3 2004/12/02 13:48:43 noro Exp $
1.1 noro 3: */
4:
5: #include "ca.h"
6: #include "base.h"
7:
8: static NumberField current_numberfield;
9: extern struct order_spec *dp_current_spec;
1.2 noro 10: void simpdalg(DAlg da,DAlg *r);
1.3 noro 11: void invdalg(DAlg a,DAlg *c);
12: void rmcontdalg(DAlg a, DAlg *c);
1.1 noro 13:
14: void setfield_dalg(NODE alist)
15: {
16: NumberField nf;
17: VL vl,vl1,vl2;
18: int n,i,dim;
19: Alg *gen;
20: P *defpoly;
21: P p;
22: Q c,iq,two;
23: DP *ps,*mb;
1.3 noro 24: DP one;
1.1 noro 25: NODE t,b,b1,b2,hlist,mblist;
26: struct order_spec *current_spec;
27:
28: nf = (NumberField)MALLOC(sizeof(struct oNumberField));
29: current_numberfield = nf;
30: vl = 0;
31: for ( t = alist; t; t = NEXT(t) ) {
32: clctalg(BDY((Alg)BDY(t)),&vl1);
33: mergev(ALG,vl,vl1,&vl2); vl = vl2;
34: }
35: for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );
36: nf->n = n;
37: nf->vl = vl;
38: nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));
39: nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));
40: current_spec = dp_current_spec;
41: STOQ(2,two);
42: create_order_spec(0,(Obj)two,&nf->spec);
43: initd(nf->spec);
44: for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {
45: ptozp(vl1->v->attr,1,&c,&defpoly[i]);
46: ptod(ALG,vl,defpoly[i],&ps[i]);
47: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
1.3 noro 48: MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
1.1 noro 49: }
1.3 noro 50: ptod(ALG,vl,(P)ONE,&one);
51: MKDAlg(one,ONE,nf->one);
52: nf->ind = b;
53: dp_mbase(hlist,&mblist);
1.1 noro 54: initd(current_spec);
55: nf->dim = dim = length(mblist);
56: nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
57: for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
1.3 noro 58: mb[i] = (DP)BDY(t);
1.1 noro 59: }
60:
1.4 ! noro 61: void qtodalg(Q q,DAlg *r)
! 62: {
! 63: NumberField nf;
! 64: Q t;
! 65: DP nm;
! 66:
! 67: if ( !(nf=current_numberfield) )
! 68: error("qtodalg : current_numberfield is not set");
! 69: if ( !q )
! 70: *r = 0;
! 71: else if ( NID(q) == N_DA )
! 72: *r = (DAlg)q;
! 73: else if ( NID(q) == N_Q ) {
! 74: if ( INT(q) ) {
! 75: muldc(CO,nf->one->nm,(P)q,&nm);
! 76: MKDAlg(nm,ONE,*r);
! 77: } else {
! 78: NTOQ(NM(q),SGN(q),t);
! 79: muldc(CO,nf->one->nm,(P)t,&nm);
! 80: NTOQ(DN(q),1,t);
! 81: MKDAlg(nm,t,*r);
! 82: }
! 83: } else
! 84: error("qtodalg : invalid argument");
! 85: }
! 86:
1.1 noro 87: void algtodalg(Alg a,DAlg *r)
88: {
89: P ap,p,p1;
1.3 noro 90: Q c,c1,d1,dn,nm;
1.1 noro 91: DP dp;
92: DAlg da;
93: NumberField nf;
94: struct order_spec *current_spec;
1.3 noro 95: VL vl,tvl,svl;
96: V v;
1.1 noro 97:
98: if ( !(nf=current_numberfield) )
99: error("algtodalg : current_numberfield is not set");
1.3 noro 100: if ( !a ) {
101: *r = 0;
102: return;
103: }
104: switch (NID((Num)a) ) {
105: case N_Q:
106: c = (Q)a;
107: if ( INT(c) ) {
108: muldc(CO,nf->one->nm,(P)c,&dp);
109: MKDAlg(dp,ONE,*r);
110: } else {
111: NTOQ(NM(c),SGN(c),c1);
112: NTOQ(DN(c),1,d1);
113: muldc(CO,nf->one->nm,(P)c1,&dp);
114: MKDAlg(dp,c1,*r);
115: }
116: break;
117: case N_A:
118: ap = (P)BDY(a);
119: ptozp(ap,1,&c,&p);
120: if ( INT(c) ) {
121: p = ap;
122: dn = ONE;
123: } else {
124: NTOQ(NM(c),SGN(c),nm);
125: NTOQ(DN(c),1,dn);
126: mulpq(p,(P)nm,&p1); p = p1;
127: }
128: current_spec = dp_current_spec; initd(nf->spec);
129: get_vars(p,&vl);
130: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
131: v = tvl->v;
132: for ( svl = nf->vl; svl; svl = NEXT(svl) )
133: if ( v == svl->v )
134: break;
135: if ( !svl )
136: error("algtodalg : incompatible numberfield");
137: }
138: ptod(ALG,nf->vl,p,&dp);
139: MKDAlg(dp,dn,da);
140: simpdalg(da,r);
141: break;
142: default:
143: error("algtodalg : invalid argument");
144: break;
1.1 noro 145: }
146: }
147:
1.2 noro 148: void dalgtoalg(DAlg da,Alg *r)
1.1 noro 149: {
1.2 noro 150: NumberField nf;
151: P p,p1;
152: Q inv;
153:
154: if ( !(nf=current_numberfield) )
1.3 noro 155: error("dalgtoalg : current_numberfield is not set");
1.2 noro 156: dtop(ALG,nf->vl,da->nm,&p);
157: invq(da->dn,&inv);
158: mulpq(p,(P)inv,&p1);
159: MKAlg(p1,*r);
1.1 noro 160: }
161:
162: void simpdalg(DAlg da,DAlg *r)
163: {
1.2 noro 164: NumberField nf;
165: DP nm;
1.3 noro 166: DAlg d;
1.2 noro 167: Q dn,dn1;
168:
169: if ( !(nf=current_numberfield) )
1.3 noro 170: error("simpdalg : current_numberfield is not set");
171: if ( !da ) {
172: *r = 0;
173: return;
174: }
1.2 noro 175: dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
176: mulq(da->dn,dn,&dn1);
1.3 noro 177: MKDAlg(nm,dn1,d);
178: rmcontdalg(d,r);
1.1 noro 179: }
180:
181: void adddalg(DAlg a,DAlg b,DAlg *c)
182: {
1.3 noro 183: NumberField nf;
184: Q dna,dnb,a1,b1,dn,g;
185: N an,bn,gn;
1.4 ! noro 186: DAlg t;
1.3 noro 187: DP ta,tb,nm;
188: struct order_spec *current_spec;
189:
190: if ( !(nf=current_numberfield) )
191: error("adddalg : current_numberfield is not set");
192: if ( !a )
193: *c = b;
194: else if ( !b )
195: *c = a;
196: else {
1.4 ! noro 197: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 198: dna = a->dn;
199: dnb = b->dn;
200: gcdn(NM(dna),NM(dnb),&gn);
201: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
202: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
203: /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
204: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
205: current_spec = dp_current_spec; initd(nf->spec);
206: addd(CO,ta,tb,&nm);
207: initd(current_spec);
208: if ( !nm )
209: *c = 0;
210: else {
211: mulq(dna,b1,&dn);
212: MKDAlg(nm,dn,*c);
213: }
214: }
1.1 noro 215: }
216:
217: void subdalg(DAlg a,DAlg b,DAlg *c)
218: {
1.3 noro 219: NumberField nf;
220: Q dna,dnb,a1,b1,dn,g;
221: N an,bn,gn;
222: DP ta,tb,nm;
1.4 ! noro 223: DAlg t;
1.3 noro 224: struct order_spec *current_spec;
225:
226: if ( !(nf=current_numberfield) )
227: error("subdalg : current_numberfield is not set");
228: if ( !a )
229: *c = b;
230: else if ( !b )
231: *c = a;
232: else {
1.4 ! noro 233: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 234: dna = a->dn;
235: dnb = b->dn;
236: gcdn(NM(dna),NM(dnb),&gn);
237: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
238: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
239: /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
240: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
241: current_spec = dp_current_spec; initd(nf->spec);
242: subd(CO,ta,tb,&nm);
243: initd(current_spec);
244: if ( !nm )
245: *c = 0;
246: else {
247: mulq(dna,b1,&dn);
248: MKDAlg(nm,dn,*c);
249: }
250: }
1.1 noro 251: }
252:
1.2 noro 253: void muldalg(DAlg a,DAlg b,DAlg *c)
254: {
1.3 noro 255: NumberField nf;
256: DP nm;
257: Q dn;
258: DAlg t;
259: struct order_spec *current_spec;
260:
261: if ( !(nf=current_numberfield) )
262: error("muldalg : current_numberfield is not set");
263: if ( !a || !b )
264: *c = 0;
265: else {
1.4 ! noro 266: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 267: current_spec = dp_current_spec; initd(nf->spec);
268: muld(CO,a->nm,b->nm,&nm);
269: initd(current_spec);
270: mulq(a->dn,b->dn,&dn);
271: MKDAlg(nm,dn,t);
272: simpdalg(t,c);
273: }
1.2 noro 274: }
275:
276:
1.1 noro 277: void divdalg(DAlg a,DAlg b,DAlg *c)
278: {
1.4 ! noro 279: DAlg inv,t;
1.3 noro 280:
1.1 noro 281: if ( !current_numberfield )
1.3 noro 282: error("divdalg : current_numberfield is not set");
283: if ( !b )
284: error("divdalg : division by 0");
285: if ( !a )
286: c = 0;
287: else {
1.4 ! noro 288: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 289: invdalg(b,&inv);
290: muldalg(a,inv,c);
291: }
292: }
293:
294: void rmcontdalg(DAlg a, DAlg *r)
295: {
296: DP u,u1;
297: Q cont,c,d;
298: N gn,cn,dn;
299:
300: if ( !a )
301: *r = a;
302: else {
303: dp_ptozp(a->nm,&u);
304: divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
305: gcdn(NM(cont),NM(a->dn),&gn);
306: divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
307: divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
308: muldc(CO,u,(P)c,&u1);
309: MKDAlg(u1,d,*r);
310: }
1.1 noro 311: }
312:
313: void invdalg(DAlg a,DAlg *c)
314: {
1.3 noro 315: NumberField nf;
316: int dim,n,i,j;
317: DP *mb;
318: DP m,d,u;
319: N ln,gn,qn;
320: DAlg *simp;
321: DAlg t,a0,r;
322: Q dn,dnsol,mul;
323: MAT mobj,sol;
324: Q **mat,**solmat;
325: MP mp0,mp;
326: int *rinfo,*cinfo;
327: struct order_spec *current_spec;
328:
329: if ( !(nf=current_numberfield) )
330: error("invdalg : current_numberfield is not set");
331: if ( !a )
332: error("invdalg : division by 0");
1.4 ! noro 333: else if ( NID(a) == N_Q ) {
! 334: invq((Q)a,&dn); *c = (DAlg)dn;
! 335: return;
! 336: }
1.3 noro 337: dim = nf->dim;
338: mb = nf->mb;
339: n = nf->n;
340: ln = ONEN;
341: MKDAlg(a->nm,ONE,a0);
342: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
343: current_spec = dp_current_spec; initd(nf->spec);
344: for ( i = dim-1; i >= 0; i-- ) {
345: m = mb[i];
346: for ( j = i+1; j < dim; j++ )
347: if ( dp_redble(m,mb[j]) )
348: break;
349: if ( j < dim ) {
350: dp_subd(m,mb[j],&d);
351: muld(CO,d,simp[j]->nm,&u);
352: MKDAlg(u,simp[j]->dn,t);
353: simpdalg(t,&simp[i]);
354: } else {
355: MKDAlg(m,ONE,t);
356: muldalg(t,a0,&simp[i]);
357: }
358: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
359: muln(NM(simp[i]->dn),qn,&ln);
360: }
361: initd(current_spec);
362: NTOQ(ln,1,dn);
363: MKMAT(mobj,dim,dim+1);
364: mat = (Q **)BDY(mobj);
365: mulq(dn,a->dn,&mat[dim-1][dim]);
366: for ( j = 0; j < dim; j++ ) {
367: divq(dn,simp[j]->dn,&mul);
368: for ( i = 0, mp = BDY(simp[j]->nm); mp && i < dim; i++ )
369: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
370: mulq(mul,(Q)mp->c,&mat[i][j]);
371: mp = NEXT(mp);
372: }
373: }
374: generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
375: solmat = (Q **)BDY(sol);
376: for ( i = 0, mp0 = 0; i < dim; i++ )
377: if ( solmat[i][0] ) {
378: NEXTMP(mp0,mp);
379: mp->c = (P)solmat[i][0];
380: mp->dl = BDY(mb[i])->dl;
381: }
382: NEXT(mp) = 0; MKDP(n,mp0,u);
383: MKDAlg(u,dnsol,r);
384: rmcontdalg(r,c);
1.1 noro 385: }
386:
387: void chsgndalg(DAlg a,DAlg *c)
388: {
1.3 noro 389: DP nm;
1.4 ! noro 390: Q t;
1.3 noro 391:
392: if ( !a ) *c = 0;
1.4 ! noro 393: else if ( NID(a) == N_Q ) {
! 394: chsgnq((Q)a,&t); *c = (DAlg)t;
! 395: } else {
1.3 noro 396: chsgnd(a->nm,&nm);
397: MKDAlg(nm,a->dn,*c);
398: }
1.1 noro 399: }
400:
1.3 noro 401: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 402: {
1.3 noro 403: NumberField nf;
404: DAlg t,z,y;
1.4 ! noro 405: Q q;
1.3 noro 406: N en,qn;
407: int r;
408:
409: if ( !(nf=current_numberfield) )
410: error("pwrdalg : current_numberfield is not set");
1.4 ! noro 411: if ( !a )
! 412: *c = !e ? (DAlg)ONE : 0;
! 413: else if ( NID(a) == N_Q ) {
! 414: pwrq((Q)a,e,&q); *c = (DAlg)q;
! 415: } else if ( !e )
1.3 noro 416: *c = nf->one;
417: else if ( UNIQ(e) )
418: *c = a;
419: else {
420: if ( SGN(e) < 0 ) {
421: invdalg(a,&t); a = t;
422: }
423: en = NM(e);
424: y = nf->one;
425: z = a;
426: while ( 1 ) {
427: r = divin(en,2,&qn); en = qn;
428: if ( r ) {
429: muldalg(z,y,&t); y = t;
430: if ( !en ) {
431: *c = y;
432: return;
433: }
434: }
435: muldalg(z,z,&t); z = t;
436: }
437: }
1.1 noro 438: }
439:
1.3 noro 440: int cmpdalg(DAlg a,DAlg b)
1.1 noro 441: {
1.3 noro 442: DAlg c;
443:
444: subdalg(a,b,&c);
445: if ( !c ) return 0;
446: else
447: return SGN((Q)BDY(c->nm)->c);
1.1 noro 448: }