![[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.16 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.15 2007/02/13 07:12:54 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.15 noro 381: if ( !da ) *r = 0;
382: else {
383: dtop(ALG,nf->vl,da->nm,&p);
384: invq(da->dn,&inv);
385: mulpq(p,(P)inv,&p1);
386: MKAlg(p1,*r);
387: }
1.1 noro 388: }
389:
390: void simpdalg(DAlg da,DAlg *r)
391: {
1.2 noro 392: NumberField nf;
393: DP nm;
1.3 noro 394: DAlg d;
1.2 noro 395: Q dn,dn1;
1.5 noro 396: struct order_spec *current_spec;
1.2 noro 397:
398: if ( !(nf=current_numberfield) )
1.3 noro 399: error("simpdalg : current_numberfield is not set");
400: if ( !da ) {
401: *r = 0;
402: return;
403: }
1.5 noro 404: current_spec = dp_current_spec; initd(nf->spec);
1.2 noro 405: dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
1.10 noro 406: if ( !nm ) *r = 0;
407: else {
408: initd(current_spec);
409: mulq(da->dn,dn,&dn1);
410: MKDAlg(nm,dn1,d);
411: rmcontdalg(d,r);
412: }
1.1 noro 413: }
414:
415: void adddalg(DAlg a,DAlg b,DAlg *c)
416: {
1.3 noro 417: NumberField nf;
418: Q dna,dnb,a1,b1,dn,g;
419: N an,bn,gn;
1.4 noro 420: DAlg t;
1.3 noro 421: DP ta,tb,nm;
422: struct order_spec *current_spec;
423:
424: if ( !(nf=current_numberfield) )
425: error("adddalg : current_numberfield is not set");
426: if ( !a )
427: *c = b;
428: else if ( !b )
429: *c = a;
430: else {
1.4 noro 431: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 432: dna = a->dn;
433: dnb = b->dn;
434: gcdn(NM(dna),NM(dnb),&gn);
435: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
436: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
437: /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
438: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
439: current_spec = dp_current_spec; initd(nf->spec);
440: addd(CO,ta,tb,&nm);
441: initd(current_spec);
442: if ( !nm )
443: *c = 0;
444: else {
445: mulq(dna,b1,&dn);
446: MKDAlg(nm,dn,*c);
447: }
448: }
1.1 noro 449: }
450:
451: void subdalg(DAlg a,DAlg b,DAlg *c)
452: {
1.3 noro 453: NumberField nf;
454: Q dna,dnb,a1,b1,dn,g;
455: N an,bn,gn;
456: DP ta,tb,nm;
1.4 noro 457: DAlg t;
1.3 noro 458: struct order_spec *current_spec;
459:
460: if ( !(nf=current_numberfield) )
461: error("subdalg : current_numberfield is not set");
462: if ( !a )
463: *c = b;
464: else if ( !b )
465: *c = a;
466: else {
1.4 noro 467: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 468: dna = a->dn;
469: dnb = b->dn;
470: gcdn(NM(dna),NM(dnb),&gn);
471: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
472: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
473: /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
474: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
475: current_spec = dp_current_spec; initd(nf->spec);
476: subd(CO,ta,tb,&nm);
477: initd(current_spec);
478: if ( !nm )
479: *c = 0;
480: else {
481: mulq(dna,b1,&dn);
482: MKDAlg(nm,dn,*c);
483: }
484: }
1.1 noro 485: }
486:
1.2 noro 487: void muldalg(DAlg a,DAlg b,DAlg *c)
488: {
1.3 noro 489: NumberField nf;
490: DP nm;
491: Q dn;
492: DAlg t;
493: struct order_spec *current_spec;
494:
495: if ( !(nf=current_numberfield) )
496: error("muldalg : current_numberfield is not set");
497: if ( !a || !b )
498: *c = 0;
499: else {
1.4 noro 500: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3 noro 501: current_spec = dp_current_spec; initd(nf->spec);
502: muld(CO,a->nm,b->nm,&nm);
503: initd(current_spec);
504: mulq(a->dn,b->dn,&dn);
505: MKDAlg(nm,dn,t);
506: simpdalg(t,c);
507: }
1.2 noro 508: }
509:
510:
1.1 noro 511: void divdalg(DAlg a,DAlg b,DAlg *c)
512: {
1.4 noro 513: DAlg inv,t;
1.8 noro 514: int ret;
1.3 noro 515:
1.1 noro 516: if ( !current_numberfield )
1.3 noro 517: error("divdalg : current_numberfield is not set");
518: if ( !b )
519: error("divdalg : division by 0");
520: if ( !a )
521: c = 0;
522: else {
1.4 noro 523: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.8 noro 524: ret = invdalg(b,&inv);
525: if ( !ret ) {
526: error("divdalg : the denominator is not invertible");
527: }
1.3 noro 528: muldalg(a,inv,c);
529: }
530: }
531:
532: void rmcontdalg(DAlg a, DAlg *r)
533: {
534: DP u,u1;
535: Q cont,c,d;
536: N gn,cn,dn;
537:
538: if ( !a )
539: *r = a;
540: else {
541: dp_ptozp(a->nm,&u);
542: divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
543: gcdn(NM(cont),NM(a->dn),&gn);
544: divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
545: divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
546: muldc(CO,u,(P)c,&u1);
547: MKDAlg(u1,d,*r);
548: }
1.1 noro 549: }
550:
1.8 noro 551: int invdalg(DAlg a,DAlg *c)
1.1 noro 552: {
1.3 noro 553: NumberField nf;
1.8 noro 554: int dim,n,i,j,k,l;
1.3 noro 555: DP *mb;
556: DP m,d,u;
557: N ln,gn,qn;
558: DAlg *simp;
559: DAlg t,a0,r;
1.7 noro 560: Q dn,dnsol,mul,nmc,dn1;
1.3 noro 561: MAT mobj,sol;
562: Q **mat,**solmat;
563: MP mp0,mp;
564: int *rinfo,*cinfo;
1.8 noro 565: int rank,nparam;
566: NODE nd0,nd,ndt;
1.3 noro 567: struct order_spec *current_spec;
1.5 noro 568: struct oEGT eg0,eg1;
569: extern struct oEGT eg_le;
1.3 noro 570:
571: if ( !(nf=current_numberfield) )
572: error("invdalg : current_numberfield is not set");
573: if ( !a )
574: error("invdalg : division by 0");
1.4 noro 575: else if ( NID(a) == N_Q ) {
576: invq((Q)a,&dn); *c = (DAlg)dn;
1.16 ! noro 577: return 1;
1.4 noro 578: }
1.3 noro 579: dim = nf->dim;
580: mb = nf->mb;
581: n = nf->n;
582: ln = ONEN;
1.7 noro 583: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
584: MKDAlg(u,ONE,a0);
1.3 noro 585: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
586: current_spec = dp_current_spec; initd(nf->spec);
1.8 noro 587: for ( i = 0; i < dim; i++ ) {
1.3 noro 588: m = mb[i];
1.8 noro 589: for ( j = i-1; j >= 0; j-- )
1.3 noro 590: if ( dp_redble(m,mb[j]) )
591: break;
1.8 noro 592: if ( j >= 0 ) {
1.3 noro 593: dp_subd(m,mb[j],&d);
594: muld(CO,d,simp[j]->nm,&u);
595: MKDAlg(u,simp[j]->dn,t);
596: simpdalg(t,&simp[i]);
597: } else {
598: MKDAlg(m,ONE,t);
599: muldalg(t,a0,&simp[i]);
600: }
601: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
602: muln(NM(simp[i]->dn),qn,&ln);
603: }
604: initd(current_spec);
605: NTOQ(ln,1,dn);
606: MKMAT(mobj,dim,dim+1);
607: mat = (Q **)BDY(mobj);
1.8 noro 608: mulq(dn,a->dn,&mat[0][dim]);
1.3 noro 609: for ( j = 0; j < dim; j++ ) {
610: divq(dn,simp[j]->dn,&mul);
1.8 noro 611: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
1.3 noro 612: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
613: mulq(mul,(Q)mp->c,&mat[i][j]);
614: mp = NEXT(mp);
615: }
616: }
1.5 noro 617: get_eg(&eg0);
1.8 noro 618: rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
1.5 noro 619: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
1.8 noro 620: if ( cinfo[0] == dim ) {
621: /* the input is invertible */
622: solmat = (Q **)BDY(sol);
623: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
624: if ( solmat[i][0] ) {
625: NEXTMP(mp0,mp);
626: mp->c = (P)solmat[i][0];
627: mp->dl = BDY(mb[i])->dl;
628: }
629: NEXT(mp) = 0; MKDP(n,mp0,u);
630: mulq(dnsol,nmc,&dn1);
631: MKDAlg(u,dn1,r);
632: rmcontdalg(r,c);
633: return 1;
634: } else
635: return 0;
636: }
637:
638: NODE inv_or_split_dalg(DAlg a,DAlg *c)
639: {
640: NumberField nf;
641: int dim,n,i,j,k,l;
642: DP *mb;
643: DP m,d,u;
644: N ln,gn,qn;
645: DAlg *simp;
646: DAlg t,a0,r;
647: Q dn,dnsol,mul,nmc,dn1;
648: MAT mobj,sol;
649: Q **mat,**solmat;
650: MP mp0,mp;
651: int *rinfo,*cinfo;
652: int rank,nparam;
653: NODE nd0,nd,ndt;
654: struct order_spec *current_spec;
655: struct oEGT eg0,eg1;
656: extern struct oEGT eg_le;
1.12 noro 657: extern int DP_Print;
1.8 noro 658:
659: if ( !(nf=current_numberfield) )
660: error("invdalg : current_numberfield is not set");
661: if ( !a )
662: error("invdalg : division by 0");
663: else if ( NID(a) == N_Q ) {
664: invq((Q)a,&dn); *c = (DAlg)dn;
1.16 ! noro 665: return 1;
1.8 noro 666: }
667: dim = nf->dim;
668: mb = nf->mb;
669: n = nf->n;
670: ln = ONEN;
671: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
672: MKDAlg(u,ONE,a0);
1.14 noro 673: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
1.8 noro 674: current_spec = dp_current_spec; initd(nf->spec);
675: for ( i = 0; i < dim; i++ ) {
1.12 noro 676: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
1.8 noro 677: m = mb[i];
678: for ( j = i-1; j >= 0; j-- )
679: if ( dp_redble(m,mb[j]) )
680: break;
681: if ( j >= 0 ) {
682: dp_subd(m,mb[j],&d);
1.11 noro 683: if ( simp[j] ) {
684: muld(CO,d,simp[j]->nm,&u);
685: MKDAlg(u,simp[j]->dn,t);
686: simpdalg(t,&simp[i]);
687: } else
688: simp[i] = 0;
1.8 noro 689: } else {
690: MKDAlg(m,ONE,t);
691: muldalg(t,a0,&simp[i]);
692: }
1.11 noro 693: if ( simp[i] ) {
694: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
695: muln(NM(simp[i]->dn),qn,&ln);
696: }
1.8 noro 697: }
698: initd(current_spec);
699: NTOQ(ln,1,dn);
700: MKMAT(mobj,dim,dim+1);
701: mat = (Q **)BDY(mobj);
702: mulq(dn,a->dn,&mat[0][dim]);
703: for ( j = 0; j < dim; j++ ) {
1.11 noro 704: if ( simp[j] ) {
705: divq(dn,simp[j]->dn,&mul);
706: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
707: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
708: mulq(mul,(Q)mp->c,&mat[i][j]);
709: mp = NEXT(mp);
710: }
711: }
1.8 noro 712: }
713: get_eg(&eg0);
1.14 noro 714: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
1.8 noro 715: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
716: if ( cinfo[0] == dim ) {
717: /* the input is invertible */
718: solmat = (Q **)BDY(sol);
719: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
720: if ( solmat[i][0] ) {
721: NEXTMP(mp0,mp);
722: mp->c = (P)solmat[i][0];
723: mp->dl = BDY(mb[i])->dl;
724: }
725: NEXT(mp) = 0; MKDP(n,mp0,u);
726: mulq(dnsol,nmc,&dn1);
727: MKDAlg(u,dn1,r);
728: rmcontdalg(r,c);
729: return 0;
730: } else {
731: /* the input is not invertible */
1.13 noro 732: nparam = sol->col;
1.8 noro 733: solmat = (Q **)BDY(sol);
734: nd0 = 0;
735: for ( k = 0; k < nparam; k++ ) {
1.13 noro 736: /* construct a new basis element */
1.8 noro 737: m = mb[cinfo[k]];
738: mp0 = 0;
1.3 noro 739: NEXTMP(mp0,mp);
1.8 noro 740: chsgnq(dnsol,&dn1); mp->c = (P)dn1;
741: mp->dl = BDY(m)->dl;
742: /* skip the last parameter */
743: for ( l = rank-2; l >= 0; l-- ) {
744: if ( solmat[l][k] ) {
745: NEXTMP(mp0,mp);
746: mp->c = (P)solmat[l][k];
747: mp->dl = BDY(mb[rinfo[l]])->dl;
748: }
749: }
750: NEXT(mp) = 0; MKDP(n,mp0,u);
751: NEXTNODE(nd0,nd);
752: BDY(nd) = (pointer)u;
753: NEXT(nd) = 0;
1.3 noro 754: }
1.8 noro 755: NEXT(nd) = 0;
756: return nd0;
757: }
1.1 noro 758: }
759:
1.14 noro 760: NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv)
761: {
762: int dim,n,i,j,k,l,nv;
763: DP *mb,*ps;
764: DP m,d,u,nm;
765: N ln,gn,qn;
766: DAlg *simp;
767: DAlg a0,r;
768: Q dn,dnsol,mul,nmc,dn1,iq;
769: MAT mobj,sol;
770: Q **mat,**solmat;
771: MP mp0,mp;
772: int *rinfo,*cinfo;
773: int rank,nparam;
774: NODE nd0,nd,ndt,ind,indt,t,mblist;
775: struct oEGT eg0,eg1;
776: extern struct oEGT eg_le;
777: extern int DP_Print;
778: initd(spec);
779: dp_ptozp(f,&u); f = u;
780:
781: n = length(gb);
782: ps = (DP *)MALLOC(n*sizeof(DP));
783: for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) {
784: ps[i] = (DP)BDY(t);
785: NEXTNODE(ind,indt);
786: STOQ(i,iq); BDY(indt) = iq;
787: }
788: if ( ind ) NEXT(indt) = 0;
789: dp_true_nf(ind,f,ps,1,&nm,&dn);
790: if ( !nm ) error("dp_inv_or_split : input is 0");
791: f = nm;
792:
793: dp_mbase(gb,&mblist);
794: dim = length(mblist);
795: mb = (DP *)MALLOC(dim*sizeof(DP));
796: for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) )
797: mb[dim-i-1] = (DP)BDY(t);
798: nv = mb[0]->nv;
799: ln = ONEN;
800: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
801: for ( i = 0; i < dim; i++ ) {
802: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
803: m = mb[i];
804: for ( j = i-1; j >= 0; j-- )
805: if ( dp_redble(m,mb[j]) )
806: break;
807: if ( j >= 0 ) {
808: dp_subd(m,mb[j],&d);
809: if ( simp[j] ) {
810: muld(CO,d,simp[j]->nm,&u);
811: dp_true_nf(ind,u,ps,1,&nm,&dn);
812: mulq(simp[j]->dn,dn,&dn1);
813: MKDAlg(nm,dn1,simp[i]);
814: } else
815: simp[i] = 0;
816: } else {
817: dp_true_nf(ind,f,ps,1,&nm,&dn);
818: MKDAlg(nm,dn,simp[i]);
819: }
820: if ( simp[i] ) {
821: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
822: muln(NM(simp[i]->dn),qn,&ln);
823: }
824: }
825: NTOQ(ln,1,dn);
826: MKMAT(mobj,dim,dim+1);
827: mat = (Q **)BDY(mobj);
828: mat[0][dim] = dn;
829: for ( j = 0; j < dim; j++ ) {
830: if ( simp[j] ) {
831: divq(dn,simp[j]->dn,&mul);
832: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
833: if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) {
834: mulq(mul,(Q)mp->c,&mat[i][j]);
835: mp = NEXT(mp);
836: }
837: }
838: }
839: get_eg(&eg0);
840: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
841: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
842: if ( cinfo[0] == dim ) {
843: /* the input is invertible */
844: solmat = (Q **)BDY(sol);
845: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
846: if ( solmat[i][0] ) {
847: NEXTMP(mp0,mp);
848: mp->c = (P)solmat[i][0];
849: mp->dl = BDY(mb[i])->dl;
850: }
851: NEXT(mp) = 0; MKDP(nv,mp0,*inv);
852: return 0;
853: } else {
854: /* the input is not invertible */
855: nparam = sol->col;
856: solmat = (Q **)BDY(sol);
857: nd0 = 0;
858: for ( k = 0; k < nparam; k++ ) {
859: /* construct a new basis element */
860: m = mb[cinfo[k]];
861: mp0 = 0;
862: NEXTMP(mp0,mp);
863: chsgnq(dnsol,&dn1); mp->c = (P)dn1;
864: mp->dl = BDY(m)->dl;
865: /* skip the last parameter */
866: for ( l = rank-2; l >= 0; l-- ) {
867: if ( solmat[l][k] ) {
868: NEXTMP(mp0,mp);
869: mp->c = (P)solmat[l][k];
870: mp->dl = BDY(mb[rinfo[l]])->dl;
871: }
872: }
873: NEXT(mp) = 0; MKDP(nv,mp0,u);
874: NEXTNODE(nd0,nd);
875: BDY(nd) = (pointer)u;
876: NEXT(nd) = 0;
877: }
878: NEXT(nd) = 0;
879: return nd0;
880: }
881: }
882:
1.1 noro 883: void chsgndalg(DAlg a,DAlg *c)
884: {
1.3 noro 885: DP nm;
1.4 noro 886: Q t;
1.3 noro 887:
888: if ( !a ) *c = 0;
1.4 noro 889: else if ( NID(a) == N_Q ) {
890: chsgnq((Q)a,&t); *c = (DAlg)t;
891: } else {
1.3 noro 892: chsgnd(a->nm,&nm);
893: MKDAlg(nm,a->dn,*c);
894: }
1.1 noro 895: }
896:
1.3 noro 897: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 898: {
1.3 noro 899: NumberField nf;
900: DAlg t,z,y;
1.4 noro 901: Q q;
1.3 noro 902: N en,qn;
903: int r;
1.8 noro 904: int ret;
1.3 noro 905:
906: if ( !(nf=current_numberfield) )
907: error("pwrdalg : current_numberfield is not set");
1.4 noro 908: if ( !a )
909: *c = !e ? (DAlg)ONE : 0;
910: else if ( NID(a) == N_Q ) {
911: pwrq((Q)a,e,&q); *c = (DAlg)q;
912: } else if ( !e )
1.3 noro 913: *c = nf->one;
914: else if ( UNIQ(e) )
915: *c = a;
916: else {
917: if ( SGN(e) < 0 ) {
1.8 noro 918: ret = invdalg(a,&t);
919: if ( !ret )
920: error("pwrdalg : the denominator is not invertible");
921: a = t;
1.3 noro 922: }
923: en = NM(e);
924: y = nf->one;
925: z = a;
926: while ( 1 ) {
927: r = divin(en,2,&qn); en = qn;
928: if ( r ) {
929: muldalg(z,y,&t); y = t;
930: if ( !en ) {
931: *c = y;
932: return;
933: }
934: }
935: muldalg(z,z,&t); z = t;
936: }
937: }
1.1 noro 938: }
939:
1.3 noro 940: int cmpdalg(DAlg a,DAlg b)
1.1 noro 941: {
1.3 noro 942: DAlg c;
943:
944: subdalg(a,b,&c);
945: if ( !c ) return 0;
946: else
947: return SGN((Q)BDY(c->nm)->c);
1.1 noro 948: }
1.9 noro 949:
950: /* convert da to a univariate poly; return the position of variable */
951:
952: int dalgtoup(DAlg da,P *up,Q *dn)
953: {
954: int nv,i,hi,current_d;
955: DCP dc0,dc;
956: MP h,mp0,mp,t;
957: DL hd,d;
958: DP c;
959: DAlg cc;
960: P v;
961:
962: nv = da->nm->nv;
963: h = BDY(da->nm);
964: *dn = da->dn;
965: hd = h->dl;
966: for ( i = 0; i < nv; i++ )
967: if ( hd->d[i] ) break;
968: hi = i;
969: current_d = hd->d[i];
970: dc0 = 0;
971: mp0 = 0;
972: for ( t = h; t; t = NEXT(t) ) {
973: NEWDL(d,nv);
974: for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
975: for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
976: d->td = t->dl->td - t->dl->d[hi];
977: if ( t->dl->d[hi] != current_d ) {
978: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
979: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
980: current_d = t->dl->d[hi];
981: mp0 = 0;
982: }
983: NEXTMP(mp0,mp);
984: mp->c = t->c; mp->dl = d;
985: }
986: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
987: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
988: NEXT(dc) = 0;
989: makevar("x",&v);
990: MKP(VR(v),dc0,*up);
991: return hi;
992: }
993: