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