![[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.12 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.11 2005/08/24 06:28:39 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.8 noro 11: int invdalg(DAlg a,DAlg *c);
1.3 noro 12: void rmcontdalg(DAlg a, DAlg *c);
1.6 noro 13: void algtodalg(Alg a,DAlg *r);
14: void dalgtoalg(DAlg da,Alg *r);
1.1 noro 15:
1.5 noro 16: NumberField get_numberfield()
17: {
18: return current_numberfield;
19: }
20:
1.1 noro 21: void setfield_dalg(NODE alist)
22: {
23: NumberField nf;
24: VL vl,vl1,vl2;
25: int n,i,dim;
26: Alg *gen;
27: P *defpoly;
28: P p;
29: Q c,iq,two;
30: DP *ps,*mb;
1.3 noro 31: DP one;
1.1 noro 32: NODE t,b,b1,b2,hlist,mblist;
33: struct order_spec *current_spec;
34:
35: nf = (NumberField)MALLOC(sizeof(struct oNumberField));
36: current_numberfield = nf;
37: vl = 0;
38: for ( t = alist; t; t = NEXT(t) ) {
39: clctalg(BDY((Alg)BDY(t)),&vl1);
40: mergev(ALG,vl,vl1,&vl2); vl = vl2;
41: }
42: for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );
43: nf->n = n;
44: nf->vl = vl;
45: nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));
46: nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));
47: current_spec = dp_current_spec;
48: STOQ(2,two);
49: create_order_spec(0,(Obj)two,&nf->spec);
50: initd(nf->spec);
51: for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {
52: ptozp(vl1->v->attr,1,&c,&defpoly[i]);
53: ptod(ALG,vl,defpoly[i],&ps[i]);
54: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
1.3 noro 55: MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
1.1 noro 56: }
1.3 noro 57: ptod(ALG,vl,(P)ONE,&one);
58: MKDAlg(one,ONE,nf->one);
59: nf->ind = b;
60: dp_mbase(hlist,&mblist);
1.1 noro 61: initd(current_spec);
62: nf->dim = dim = length(mblist);
63: nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
64: for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
1.8 noro 65: mb[dim-i-1] = (DP)BDY(t);
1.1 noro 66: }
67:
1.10 noro 68: void setfield_gb(NODE gb,VL vl,struct order_spec *spec)
69: {
70: NumberField nf;
71: VL vl1,vl2;
72: int n,i,dim;
73: Alg *gen;
74: P *defpoly;
75: P p;
76: Q c,iq,two;
77: DP *ps,*mb;
78: DP one;
79: NODE t,b,b1,b2,hlist,mblist;
80: struct order_spec *current_spec;
81:
82: nf = (NumberField)MALLOC(sizeof(struct oNumberField));
83: current_numberfield = nf;
84: for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );
85: nf->n = n;
86: nf->psn = length(gb);
87: nf->vl = vl;
88: nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));
89: nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));
90: current_spec = dp_current_spec;
91: nf->spec = spec;
92: initd(nf->spec);
93: for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {
94: ptozp((P)BDY(t),1,&c,&defpoly[i]);
95: ptod(CO,vl,defpoly[i],&ps[i]);
96: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
97: MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
98: }
99: ptod(ALG,vl,(P)ONE,&one);
100: MKDAlg(one,ONE,nf->one);
101: nf->ind = b;
102: dp_mbase(hlist,&mblist);
103: initd(current_spec);
104: nf->dim = dim = length(mblist);
105: nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
106: for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
107: mb[dim-i-1] = (DP)BDY(t);
108: }
109:
1.4 noro 110: void qtodalg(Q q,DAlg *r)
111: {
112: NumberField nf;
113: Q t;
114: DP nm;
115:
116: if ( !(nf=current_numberfield) )
117: error("qtodalg : current_numberfield is not set");
118: if ( !q )
119: *r = 0;
120: else if ( NID(q) == N_DA )
121: *r = (DAlg)q;
122: else if ( NID(q) == N_Q ) {
123: if ( INT(q) ) {
124: muldc(CO,nf->one->nm,(P)q,&nm);
125: MKDAlg(nm,ONE,*r);
126: } else {
127: NTOQ(NM(q),SGN(q),t);
128: muldc(CO,nf->one->nm,(P)t,&nm);
129: NTOQ(DN(q),1,t);
130: MKDAlg(nm,t,*r);
131: }
132: } else
133: error("qtodalg : invalid argument");
1.6 noro 134: }
135:
136: void obj_algtodalg(Obj obj,Obj *r)
137: {
138: DAlg d;
139: DCP dc,dcr0,dcr;
140: P c,p;
141: Obj t;
142: Obj nm,dn;
143: NODE b,s,s0;
144: R rat;
145: VECT v;
146: MAT mat;
147: LIST list;
148: pointer *a;
149: pointer **m;
150: int len,row,col,i,j,l;
151:
152: if ( !obj ) {
153: *r = 0;
154: return;
155: }
156: switch ( OID(obj) ) {
157: case O_N:
158: algtodalg((Alg)obj,&d); *r = (Obj)d;
159: break;
160: case O_P:
161: for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
162: obj_algtodalg((Obj)COEF(dc),&t);
163: if ( t ) {
164: NEXTDC(dcr0,dcr);
165: COEF(dcr) = (P)t;
166: DEG(dcr) = DEG(dc);
167: }
168: }
169: if ( dcr0 ) {
170: MKP(VR((P)obj),dcr0,p);
171: *r = (Obj)p;
172: } else
173: *r = 0;
174: break;
175: case O_R:
176: obj_algtodalg((Obj)NM((R)obj),&nm);
177: obj_algtodalg((Obj)DN((R)obj),&dn);
178: if ( !dn )
179: error("obj_algtodalg : division by 0");
180: if ( !nm )
181: *r = 0;
182: else {
183: MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
184: }
185: break;
186: case O_LIST:
187: s0 = 0;
188: for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
189: NEXTNODE(s0,s);
190: obj_algtodalg((Obj)BDY(b),&t);
191: BDY(s) = (pointer)t;
192: }
193: NEXT(s) = 0;
194: MKLIST(list,s0);
195: *r = (Obj)list;
196: break;
197: case O_VECT:
198: l = ((VECT)obj)->len;
199: a = BDY((VECT)obj);
200: MKVECT(v,l);
201: for ( i = 0; i < l; i++ ) {
202: obj_algtodalg((Obj)a[i],&t);
203: BDY(v)[i] = (pointer)t;
204: }
205: *r = (Obj)v;
206: break;
207: case O_MAT:
208: row = ((MAT)obj)->row; col = ((MAT)obj)->col;
209: m = BDY((MAT)obj);
210: MKMAT(mat,row,col);
211: for ( i = 0; i < row; i++ )
212: for ( j = 0; j < col; j++ ) {
213: obj_algtodalg((Obj)m[i][j],&t);
214: BDY(mat)[i][j] = (pointer)t;
215: }
216: *r = (Obj)mat;
217: break;
218: default:
219: *r = obj;
220: break;
221: }
222: }
223:
224: void obj_dalgtoalg(Obj obj,Obj *r)
225: {
226: Alg d;
227: DCP dc,dcr0,dcr;
228: P c,p;
229: Obj t;
230: Obj nm,dn;
231: NODE b,s,s0;
232: R rat;
233: VECT v;
234: MAT mat;
235: LIST list;
236: pointer *a;
237: pointer **m;
238: int len,row,col,i,j,l;
239:
240: if ( !obj ) {
241: *r = 0;
242: return;
243: }
244: switch ( OID(obj) ) {
245: case O_N:
246: dalgtoalg((DAlg)obj,&d); *r = (Obj)d;
247: break;
248: case O_P:
249: for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
250: obj_dalgtoalg((Obj)COEF(dc),&t);
251: if ( t ) {
252: NEXTDC(dcr0,dcr);
253: COEF(dcr) = (P)t;
254: DEG(dcr) = DEG(dc);
255: }
256: }
257: if ( dcr0 ) {
258: MKP(VR((P)obj),dcr0,p);
259: *r = (Obj)p;
260: } else
261: *r = 0;
262: break;
263: case O_R:
264: obj_dalgtoalg((Obj)NM((R)obj),&nm);
265: obj_dalgtoalg((Obj)DN((R)obj),&dn);
266: if ( !dn )
267: error("obj_dalgtoalg : division by 0");
268: if ( !nm )
269: *r = 0;
270: else {
271: MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
272: }
273: break;
274: case O_LIST:
275: s0 = 0;
276: for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
277: NEXTNODE(s0,s);
278: obj_dalgtoalg((Obj)BDY(b),&t);
279: BDY(s) = (pointer)t;
280: }
281: NEXT(s) = 0;
282: MKLIST(list,s0);
283: *r = (Obj)list;
284: break;
285: case O_VECT:
286: l = ((VECT)obj)->len;
287: a = BDY((VECT)obj);
288: MKVECT(v,l);
289: for ( i = 0; i < l; i++ ) {
290: obj_dalgtoalg((Obj)a[i],&t);
291: BDY(v)[i] = (pointer)t;
292: }
293: *r = (Obj)v;
294: break;
295: case O_MAT:
296: row = ((MAT)obj)->row; col = ((MAT)obj)->col;
297: m = BDY((MAT)obj);
298: MKMAT(mat,row,col);
299: for ( i = 0; i < row; i++ )
300: for ( j = 0; j < col; j++ ) {
301: obj_dalgtoalg((Obj)m[i][j],&t);
302: BDY(mat)[i][j] = (pointer)t;
303: }
304: *r = (Obj)mat;
305: break;
306: default:
307: *r = obj;
308: break;
309: }
1.4 noro 310: }
311:
1.1 noro 312: void algtodalg(Alg a,DAlg *r)
313: {
314: P ap,p,p1;
1.3 noro 315: Q c,c1,d1,dn,nm;
1.1 noro 316: DP dp;
317: DAlg da;
318: NumberField nf;
319: struct order_spec *current_spec;
1.3 noro 320: VL vl,tvl,svl;
321: V v;
1.1 noro 322:
323: if ( !(nf=current_numberfield) )
324: error("algtodalg : current_numberfield is not set");
1.3 noro 325: if ( !a ) {
326: *r = 0;
327: return;
328: }
329: switch (NID((Num)a) ) {
330: case N_Q:
331: c = (Q)a;
332: if ( INT(c) ) {
333: muldc(CO,nf->one->nm,(P)c,&dp);
334: MKDAlg(dp,ONE,*r);
335: } else {
336: NTOQ(NM(c),SGN(c),c1);
337: NTOQ(DN(c),1,d1);
338: muldc(CO,nf->one->nm,(P)c1,&dp);
1.9 noro 339: MKDAlg(dp,d1,*r);
1.3 noro 340: }
341: break;
342: case N_A:
343: ap = (P)BDY(a);
344: ptozp(ap,1,&c,&p);
345: if ( INT(c) ) {
346: p = ap;
347: dn = ONE;
348: } else {
349: NTOQ(NM(c),SGN(c),nm);
350: NTOQ(DN(c),1,dn);
351: mulpq(p,(P)nm,&p1); p = p1;
352: }
353: current_spec = dp_current_spec; initd(nf->spec);
354: get_vars(p,&vl);
355: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
356: v = tvl->v;
357: for ( svl = nf->vl; svl; svl = NEXT(svl) )
358: if ( v == svl->v )
359: break;
360: if ( !svl )
361: error("algtodalg : incompatible numberfield");
362: }
363: ptod(ALG,nf->vl,p,&dp);
364: MKDAlg(dp,dn,da);
365: simpdalg(da,r);
366: break;
367: default:
368: error("algtodalg : invalid argument");
369: break;
1.1 noro 370: }
371: }
372:
1.2 noro 373: void dalgtoalg(DAlg da,Alg *r)
1.1 noro 374: {
1.2 noro 375: NumberField nf;
376: P p,p1;
377: Q inv;
378:
379: if ( !(nf=current_numberfield) )
1.3 noro 380: error("dalgtoalg : current_numberfield is not set");
1.2 noro 381: dtop(ALG,nf->vl,da->nm,&p);
382: invq(da->dn,&inv);
383: mulpq(p,(P)inv,&p1);
384: MKAlg(p1,*r);
1.1 noro 385: }
386:
387: void simpdalg(DAlg da,DAlg *r)
388: {
1.2 noro 389: NumberField nf;
390: DP nm;
1.3 noro 391: DAlg d;
1.2 noro 392: Q dn,dn1;
1.5 noro 393: struct order_spec *current_spec;
1.2 noro 394:
395: if ( !(nf=current_numberfield) )
1.3 noro 396: error("simpdalg : current_numberfield is not set");
397: if ( !da ) {
398: *r = 0;
399: return;
400: }
1.5 noro 401: current_spec = dp_current_spec; initd(nf->spec);
1.2 noro 402: dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
1.10 noro 403: if ( !nm ) *r = 0;
404: else {
405: initd(current_spec);
406: mulq(da->dn,dn,&dn1);
407: MKDAlg(nm,dn1,d);
408: rmcontdalg(d,r);
409: }
1.1 noro 410: }
411:
412: void adddalg(DAlg a,DAlg b,DAlg *c)
413: {
1.3 noro 414: NumberField nf;
415: Q dna,dnb,a1,b1,dn,g;
416: N an,bn,gn;
1.4 noro 417: DAlg t;
1.3 noro 418: DP ta,tb,nm;
419: struct order_spec *current_spec;
420:
421: if ( !(nf=current_numberfield) )
422: error("adddalg : current_numberfield is not set");
423: if ( !a )
424: *c = b;
425: else if ( !b )
426: *c = a;
427: else {
1.4 noro 428: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 429: dna = a->dn;
430: dnb = b->dn;
431: gcdn(NM(dna),NM(dnb),&gn);
432: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
433: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
434: /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
435: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
436: current_spec = dp_current_spec; initd(nf->spec);
437: addd(CO,ta,tb,&nm);
438: initd(current_spec);
439: if ( !nm )
440: *c = 0;
441: else {
442: mulq(dna,b1,&dn);
443: MKDAlg(nm,dn,*c);
444: }
445: }
1.1 noro 446: }
447:
448: void subdalg(DAlg a,DAlg b,DAlg *c)
449: {
1.3 noro 450: NumberField nf;
451: Q dna,dnb,a1,b1,dn,g;
452: N an,bn,gn;
453: DP ta,tb,nm;
1.4 noro 454: DAlg t;
1.3 noro 455: struct order_spec *current_spec;
456:
457: if ( !(nf=current_numberfield) )
458: error("subdalg : current_numberfield is not set");
459: if ( !a )
460: *c = b;
461: else if ( !b )
462: *c = a;
463: else {
1.4 noro 464: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 465: dna = a->dn;
466: dnb = b->dn;
467: gcdn(NM(dna),NM(dnb),&gn);
468: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
469: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
470: /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
471: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
472: current_spec = dp_current_spec; initd(nf->spec);
473: subd(CO,ta,tb,&nm);
474: initd(current_spec);
475: if ( !nm )
476: *c = 0;
477: else {
478: mulq(dna,b1,&dn);
479: MKDAlg(nm,dn,*c);
480: }
481: }
1.1 noro 482: }
483:
1.2 noro 484: void muldalg(DAlg a,DAlg b,DAlg *c)
485: {
1.3 noro 486: NumberField nf;
487: DP nm;
488: Q dn;
489: DAlg t;
490: struct order_spec *current_spec;
491:
492: if ( !(nf=current_numberfield) )
493: error("muldalg : current_numberfield is not set");
494: if ( !a || !b )
495: *c = 0;
496: else {
1.4 noro 497: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 498: current_spec = dp_current_spec; initd(nf->spec);
499: muld(CO,a->nm,b->nm,&nm);
500: initd(current_spec);
501: mulq(a->dn,b->dn,&dn);
502: MKDAlg(nm,dn,t);
503: simpdalg(t,c);
504: }
1.2 noro 505: }
506:
507:
1.1 noro 508: void divdalg(DAlg a,DAlg b,DAlg *c)
509: {
1.4 noro 510: DAlg inv,t;
1.8 noro 511: int ret;
1.3 noro 512:
1.1 noro 513: if ( !current_numberfield )
1.3 noro 514: error("divdalg : current_numberfield is not set");
515: if ( !b )
516: error("divdalg : division by 0");
517: if ( !a )
518: c = 0;
519: else {
1.4 noro 520: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.8 noro 521: ret = invdalg(b,&inv);
522: if ( !ret ) {
523: error("divdalg : the denominator is not invertible");
524: }
1.3 noro 525: muldalg(a,inv,c);
526: }
527: }
528:
529: void rmcontdalg(DAlg a, DAlg *r)
530: {
531: DP u,u1;
532: Q cont,c,d;
533: N gn,cn,dn;
534:
535: if ( !a )
536: *r = a;
537: else {
538: dp_ptozp(a->nm,&u);
539: divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
540: gcdn(NM(cont),NM(a->dn),&gn);
541: divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
542: divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
543: muldc(CO,u,(P)c,&u1);
544: MKDAlg(u1,d,*r);
545: }
1.1 noro 546: }
547:
1.8 noro 548: int invdalg(DAlg a,DAlg *c)
1.1 noro 549: {
1.3 noro 550: NumberField nf;
1.8 noro 551: int dim,n,i,j,k,l;
1.3 noro 552: DP *mb;
553: DP m,d,u;
554: N ln,gn,qn;
555: DAlg *simp;
556: DAlg t,a0,r;
1.7 noro 557: Q dn,dnsol,mul,nmc,dn1;
1.3 noro 558: MAT mobj,sol;
559: Q **mat,**solmat;
560: MP mp0,mp;
561: int *rinfo,*cinfo;
1.8 noro 562: int rank,nparam;
563: NODE nd0,nd,ndt;
1.3 noro 564: struct order_spec *current_spec;
1.5 noro 565: struct oEGT eg0,eg1;
566: extern struct oEGT eg_le;
1.3 noro 567:
568: if ( !(nf=current_numberfield) )
569: error("invdalg : current_numberfield is not set");
570: if ( !a )
571: error("invdalg : division by 0");
1.4 noro 572: else if ( NID(a) == N_Q ) {
573: invq((Q)a,&dn); *c = (DAlg)dn;
574: return;
575: }
1.3 noro 576: dim = nf->dim;
577: mb = nf->mb;
578: n = nf->n;
579: ln = ONEN;
1.7 noro 580: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
581: MKDAlg(u,ONE,a0);
1.3 noro 582: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
583: current_spec = dp_current_spec; initd(nf->spec);
1.8 noro 584: for ( i = 0; i < dim; i++ ) {
1.3 noro 585: m = mb[i];
1.8 noro 586: for ( j = i-1; j >= 0; j-- )
1.3 noro 587: if ( dp_redble(m,mb[j]) )
588: break;
1.8 noro 589: if ( j >= 0 ) {
1.3 noro 590: dp_subd(m,mb[j],&d);
591: muld(CO,d,simp[j]->nm,&u);
592: MKDAlg(u,simp[j]->dn,t);
593: simpdalg(t,&simp[i]);
594: } else {
595: MKDAlg(m,ONE,t);
596: muldalg(t,a0,&simp[i]);
597: }
598: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
599: muln(NM(simp[i]->dn),qn,&ln);
600: }
601: initd(current_spec);
602: NTOQ(ln,1,dn);
603: MKMAT(mobj,dim,dim+1);
604: mat = (Q **)BDY(mobj);
1.8 noro 605: mulq(dn,a->dn,&mat[0][dim]);
1.3 noro 606: for ( j = 0; j < dim; j++ ) {
607: divq(dn,simp[j]->dn,&mul);
1.8 noro 608: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
1.3 noro 609: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
610: mulq(mul,(Q)mp->c,&mat[i][j]);
611: mp = NEXT(mp);
612: }
613: }
1.5 noro 614: get_eg(&eg0);
1.8 noro 615: rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
1.5 noro 616: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
1.8 noro 617: if ( cinfo[0] == dim ) {
618: /* the input is invertible */
619: solmat = (Q **)BDY(sol);
620: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
621: if ( solmat[i][0] ) {
622: NEXTMP(mp0,mp);
623: mp->c = (P)solmat[i][0];
624: mp->dl = BDY(mb[i])->dl;
625: }
626: NEXT(mp) = 0; MKDP(n,mp0,u);
627: mulq(dnsol,nmc,&dn1);
628: MKDAlg(u,dn1,r);
629: rmcontdalg(r,c);
630: return 1;
631: } else
632: return 0;
633: }
634:
635: NODE inv_or_split_dalg(DAlg a,DAlg *c)
636: {
637: NumberField nf;
638: int dim,n,i,j,k,l;
639: DP *mb;
640: DP m,d,u;
641: N ln,gn,qn;
642: DAlg *simp;
643: DAlg t,a0,r;
644: Q dn,dnsol,mul,nmc,dn1;
645: MAT mobj,sol;
646: Q **mat,**solmat;
647: MP mp0,mp;
648: int *rinfo,*cinfo;
649: int rank,nparam;
650: NODE nd0,nd,ndt;
651: struct order_spec *current_spec;
652: struct oEGT eg0,eg1;
653: extern struct oEGT eg_le;
1.12 ! noro 654: extern int DP_Print;
1.8 noro 655:
656: if ( !(nf=current_numberfield) )
657: error("invdalg : current_numberfield is not set");
658: if ( !a )
659: error("invdalg : division by 0");
660: else if ( NID(a) == N_Q ) {
661: invq((Q)a,&dn); *c = (DAlg)dn;
662: return;
663: }
664: dim = nf->dim;
665: mb = nf->mb;
666: n = nf->n;
667: ln = ONEN;
668: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
669: MKDAlg(u,ONE,a0);
670: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
671: current_spec = dp_current_spec; initd(nf->spec);
672: for ( i = 0; i < dim; i++ ) {
1.12 ! noro 673: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
1.8 noro 674: m = mb[i];
675: for ( j = i-1; j >= 0; j-- )
676: if ( dp_redble(m,mb[j]) )
677: break;
678: if ( j >= 0 ) {
679: dp_subd(m,mb[j],&d);
1.11 noro 680: if ( simp[j] ) {
681: muld(CO,d,simp[j]->nm,&u);
682: MKDAlg(u,simp[j]->dn,t);
683: simpdalg(t,&simp[i]);
684: } else
685: simp[i] = 0;
1.8 noro 686: } else {
687: MKDAlg(m,ONE,t);
688: muldalg(t,a0,&simp[i]);
689: }
1.11 noro 690: if ( simp[i] ) {
691: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
692: muln(NM(simp[i]->dn),qn,&ln);
693: }
1.8 noro 694: }
695: initd(current_spec);
696: NTOQ(ln,1,dn);
697: MKMAT(mobj,dim,dim+1);
698: mat = (Q **)BDY(mobj);
699: mulq(dn,a->dn,&mat[0][dim]);
700: for ( j = 0; j < dim; j++ ) {
1.11 noro 701: if ( simp[j] ) {
702: divq(dn,simp[j]->dn,&mul);
703: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
704: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
705: mulq(mul,(Q)mp->c,&mat[i][j]);
706: mp = NEXT(mp);
707: }
708: }
1.8 noro 709: }
710: get_eg(&eg0);
711: rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
712: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
713: if ( cinfo[0] == dim ) {
714: /* the input is invertible */
715: solmat = (Q **)BDY(sol);
716: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
717: if ( solmat[i][0] ) {
718: NEXTMP(mp0,mp);
719: mp->c = (P)solmat[i][0];
720: mp->dl = BDY(mb[i])->dl;
721: }
722: NEXT(mp) = 0; MKDP(n,mp0,u);
723: mulq(dnsol,nmc,&dn1);
724: MKDAlg(u,dn1,r);
725: rmcontdalg(r,c);
726: return 0;
727: } else {
728: /* the input is not invertible */
729: nparam = (dim+1)-rank;
730: /* the index 'dim' should not be in cinfo[] */
731: solmat = (Q **)BDY(sol);
732: for ( k = 0; k < nparam; k++ )
733: if ( cinfo[k] == dim )
734: error("invdalg : cannot happen");
735: nd0 = 0;
736: for ( k = 0; k < nparam; k++ ) {
737: m = mb[cinfo[k]];
738: for ( ndt = nd0; ndt; ndt = NEXT(ndt) ) {
739: if ( dp_redble(m,(DP)BDY(ndt)) ) break;
740: }
741: /* skip a redundunt basis element */
742: if ( ndt ) continue;
743: /* construct a new basis element */
744: mp0 = 0;
1.3 noro 745: NEXTMP(mp0,mp);
1.8 noro 746: chsgnq(dnsol,&dn1); mp->c = (P)dn1;
747: mp->dl = BDY(m)->dl;
748: /* skip the last parameter */
749: for ( l = rank-2; l >= 0; l-- ) {
750: if ( solmat[l][k] ) {
751: NEXTMP(mp0,mp);
752: mp->c = (P)solmat[l][k];
753: mp->dl = BDY(mb[rinfo[l]])->dl;
754: }
755: }
756: NEXT(mp) = 0; MKDP(n,mp0,u);
757: NEXTNODE(nd0,nd);
758: BDY(nd) = (pointer)u;
759: NEXT(nd) = 0;
1.3 noro 760: }
1.8 noro 761: NEXT(nd) = 0;
762: return nd0;
763: }
1.1 noro 764: }
765:
766: void chsgndalg(DAlg a,DAlg *c)
767: {
1.3 noro 768: DP nm;
1.4 noro 769: Q t;
1.3 noro 770:
771: if ( !a ) *c = 0;
1.4 noro 772: else if ( NID(a) == N_Q ) {
773: chsgnq((Q)a,&t); *c = (DAlg)t;
774: } else {
1.3 noro 775: chsgnd(a->nm,&nm);
776: MKDAlg(nm,a->dn,*c);
777: }
1.1 noro 778: }
779:
1.3 noro 780: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 781: {
1.3 noro 782: NumberField nf;
783: DAlg t,z,y;
1.4 noro 784: Q q;
1.3 noro 785: N en,qn;
786: int r;
1.8 noro 787: int ret;
1.3 noro 788:
789: if ( !(nf=current_numberfield) )
790: error("pwrdalg : current_numberfield is not set");
1.4 noro 791: if ( !a )
792: *c = !e ? (DAlg)ONE : 0;
793: else if ( NID(a) == N_Q ) {
794: pwrq((Q)a,e,&q); *c = (DAlg)q;
795: } else if ( !e )
1.3 noro 796: *c = nf->one;
797: else if ( UNIQ(e) )
798: *c = a;
799: else {
800: if ( SGN(e) < 0 ) {
1.8 noro 801: ret = invdalg(a,&t);
802: if ( !ret )
803: error("pwrdalg : the denominator is not invertible");
804: a = t;
1.3 noro 805: }
806: en = NM(e);
807: y = nf->one;
808: z = a;
809: while ( 1 ) {
810: r = divin(en,2,&qn); en = qn;
811: if ( r ) {
812: muldalg(z,y,&t); y = t;
813: if ( !en ) {
814: *c = y;
815: return;
816: }
817: }
818: muldalg(z,z,&t); z = t;
819: }
820: }
1.1 noro 821: }
822:
1.3 noro 823: int cmpdalg(DAlg a,DAlg b)
1.1 noro 824: {
1.3 noro 825: DAlg c;
826:
827: subdalg(a,b,&c);
828: if ( !c ) return 0;
829: else
830: return SGN((Q)BDY(c->nm)->c);
1.1 noro 831: }
1.9 noro 832:
833: /* convert da to a univariate poly; return the position of variable */
834:
835: int dalgtoup(DAlg da,P *up,Q *dn)
836: {
837: int nv,i,hi,current_d;
838: DCP dc0,dc;
839: MP h,mp0,mp,t;
840: DL hd,d;
841: DP c;
842: DAlg cc;
843: P v;
844:
845: nv = da->nm->nv;
846: h = BDY(da->nm);
847: *dn = da->dn;
848: hd = h->dl;
849: for ( i = 0; i < nv; i++ )
850: if ( hd->d[i] ) break;
851: hi = i;
852: current_d = hd->d[i];
853: dc0 = 0;
854: mp0 = 0;
855: for ( t = h; t; t = NEXT(t) ) {
856: NEWDL(d,nv);
857: for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
858: for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
859: d->td = t->dl->td - t->dl->d[hi];
860: if ( t->dl->d[hi] != current_d ) {
861: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
862: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
863: current_d = t->dl->d[hi];
864: mp0 = 0;
865: }
866: NEXTMP(mp0,mp);
867: mp->c = t->c; mp->dl = d;
868: }
869: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
870: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
871: NEXT(dc) = 0;
872: makevar("x",&v);
873: MKP(VR(v),dc0,*up);
874: return hi;
875: }
876: