![[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.14 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.13 2006/01/05 00:21:20 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);
1.14 ! noro 670: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
1.8 noro 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);
1.14 ! noro 711: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
1.8 noro 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 */
1.13 noro 729: nparam = sol->col;
1.8 noro 730: solmat = (Q **)BDY(sol);
731: nd0 = 0;
732: for ( k = 0; k < nparam; k++ ) {
1.13 noro 733: /* construct a new basis element */
1.8 noro 734: m = mb[cinfo[k]];
735: mp0 = 0;
1.3 noro 736: NEXTMP(mp0,mp);
1.8 noro 737: chsgnq(dnsol,&dn1); mp->c = (P)dn1;
738: mp->dl = BDY(m)->dl;
739: /* skip the last parameter */
740: for ( l = rank-2; l >= 0; l-- ) {
741: if ( solmat[l][k] ) {
742: NEXTMP(mp0,mp);
743: mp->c = (P)solmat[l][k];
744: mp->dl = BDY(mb[rinfo[l]])->dl;
745: }
746: }
747: NEXT(mp) = 0; MKDP(n,mp0,u);
748: NEXTNODE(nd0,nd);
749: BDY(nd) = (pointer)u;
750: NEXT(nd) = 0;
1.3 noro 751: }
1.8 noro 752: NEXT(nd) = 0;
753: return nd0;
754: }
1.1 noro 755: }
756:
1.14 ! noro 757: NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv)
! 758: {
! 759: int dim,n,i,j,k,l,nv;
! 760: DP *mb,*ps;
! 761: DP m,d,u,nm;
! 762: N ln,gn,qn;
! 763: DAlg *simp;
! 764: DAlg a0,r;
! 765: Q dn,dnsol,mul,nmc,dn1,iq;
! 766: MAT mobj,sol;
! 767: Q **mat,**solmat;
! 768: MP mp0,mp;
! 769: int *rinfo,*cinfo;
! 770: int rank,nparam;
! 771: NODE nd0,nd,ndt,ind,indt,t,mblist;
! 772: struct oEGT eg0,eg1;
! 773: extern struct oEGT eg_le;
! 774: extern int DP_Print;
! 775: initd(spec);
! 776: dp_ptozp(f,&u); f = u;
! 777:
! 778: n = length(gb);
! 779: ps = (DP *)MALLOC(n*sizeof(DP));
! 780: for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) {
! 781: ps[i] = (DP)BDY(t);
! 782: NEXTNODE(ind,indt);
! 783: STOQ(i,iq); BDY(indt) = iq;
! 784: }
! 785: if ( ind ) NEXT(indt) = 0;
! 786: dp_true_nf(ind,f,ps,1,&nm,&dn);
! 787: if ( !nm ) error("dp_inv_or_split : input is 0");
! 788: f = nm;
! 789:
! 790: dp_mbase(gb,&mblist);
! 791: dim = length(mblist);
! 792: mb = (DP *)MALLOC(dim*sizeof(DP));
! 793: for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) )
! 794: mb[dim-i-1] = (DP)BDY(t);
! 795: nv = mb[0]->nv;
! 796: ln = ONEN;
! 797: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
! 798: for ( i = 0; i < dim; i++ ) {
! 799: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
! 800: m = mb[i];
! 801: for ( j = i-1; j >= 0; j-- )
! 802: if ( dp_redble(m,mb[j]) )
! 803: break;
! 804: if ( j >= 0 ) {
! 805: dp_subd(m,mb[j],&d);
! 806: if ( simp[j] ) {
! 807: muld(CO,d,simp[j]->nm,&u);
! 808: dp_true_nf(ind,u,ps,1,&nm,&dn);
! 809: mulq(simp[j]->dn,dn,&dn1);
! 810: MKDAlg(nm,dn1,simp[i]);
! 811: } else
! 812: simp[i] = 0;
! 813: } else {
! 814: dp_true_nf(ind,f,ps,1,&nm,&dn);
! 815: MKDAlg(nm,dn,simp[i]);
! 816: }
! 817: if ( simp[i] ) {
! 818: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
! 819: muln(NM(simp[i]->dn),qn,&ln);
! 820: }
! 821: }
! 822: NTOQ(ln,1,dn);
! 823: MKMAT(mobj,dim,dim+1);
! 824: mat = (Q **)BDY(mobj);
! 825: mat[0][dim] = dn;
! 826: for ( j = 0; j < dim; j++ ) {
! 827: if ( simp[j] ) {
! 828: divq(dn,simp[j]->dn,&mul);
! 829: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
! 830: if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) {
! 831: mulq(mul,(Q)mp->c,&mat[i][j]);
! 832: mp = NEXT(mp);
! 833: }
! 834: }
! 835: }
! 836: get_eg(&eg0);
! 837: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
! 838: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
! 839: if ( cinfo[0] == dim ) {
! 840: /* the input is invertible */
! 841: solmat = (Q **)BDY(sol);
! 842: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
! 843: if ( solmat[i][0] ) {
! 844: NEXTMP(mp0,mp);
! 845: mp->c = (P)solmat[i][0];
! 846: mp->dl = BDY(mb[i])->dl;
! 847: }
! 848: NEXT(mp) = 0; MKDP(nv,mp0,*inv);
! 849: return 0;
! 850: } else {
! 851: /* the input is not invertible */
! 852: nparam = sol->col;
! 853: solmat = (Q **)BDY(sol);
! 854: nd0 = 0;
! 855: for ( k = 0; k < nparam; k++ ) {
! 856: /* construct a new basis element */
! 857: m = mb[cinfo[k]];
! 858: mp0 = 0;
! 859: NEXTMP(mp0,mp);
! 860: chsgnq(dnsol,&dn1); mp->c = (P)dn1;
! 861: mp->dl = BDY(m)->dl;
! 862: /* skip the last parameter */
! 863: for ( l = rank-2; l >= 0; l-- ) {
! 864: if ( solmat[l][k] ) {
! 865: NEXTMP(mp0,mp);
! 866: mp->c = (P)solmat[l][k];
! 867: mp->dl = BDY(mb[rinfo[l]])->dl;
! 868: }
! 869: }
! 870: NEXT(mp) = 0; MKDP(nv,mp0,u);
! 871: NEXTNODE(nd0,nd);
! 872: BDY(nd) = (pointer)u;
! 873: NEXT(nd) = 0;
! 874: }
! 875: NEXT(nd) = 0;
! 876: return nd0;
! 877: }
! 878: }
! 879:
1.1 noro 880: void chsgndalg(DAlg a,DAlg *c)
881: {
1.3 noro 882: DP nm;
1.4 noro 883: Q t;
1.3 noro 884:
885: if ( !a ) *c = 0;
1.4 noro 886: else if ( NID(a) == N_Q ) {
887: chsgnq((Q)a,&t); *c = (DAlg)t;
888: } else {
1.3 noro 889: chsgnd(a->nm,&nm);
890: MKDAlg(nm,a->dn,*c);
891: }
1.1 noro 892: }
893:
1.3 noro 894: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 895: {
1.3 noro 896: NumberField nf;
897: DAlg t,z,y;
1.4 noro 898: Q q;
1.3 noro 899: N en,qn;
900: int r;
1.8 noro 901: int ret;
1.3 noro 902:
903: if ( !(nf=current_numberfield) )
904: error("pwrdalg : current_numberfield is not set");
1.4 noro 905: if ( !a )
906: *c = !e ? (DAlg)ONE : 0;
907: else if ( NID(a) == N_Q ) {
908: pwrq((Q)a,e,&q); *c = (DAlg)q;
909: } else if ( !e )
1.3 noro 910: *c = nf->one;
911: else if ( UNIQ(e) )
912: *c = a;
913: else {
914: if ( SGN(e) < 0 ) {
1.8 noro 915: ret = invdalg(a,&t);
916: if ( !ret )
917: error("pwrdalg : the denominator is not invertible");
918: a = t;
1.3 noro 919: }
920: en = NM(e);
921: y = nf->one;
922: z = a;
923: while ( 1 ) {
924: r = divin(en,2,&qn); en = qn;
925: if ( r ) {
926: muldalg(z,y,&t); y = t;
927: if ( !en ) {
928: *c = y;
929: return;
930: }
931: }
932: muldalg(z,z,&t); z = t;
933: }
934: }
1.1 noro 935: }
936:
1.3 noro 937: int cmpdalg(DAlg a,DAlg b)
1.1 noro 938: {
1.3 noro 939: DAlg c;
940:
941: subdalg(a,b,&c);
942: if ( !c ) return 0;
943: else
944: return SGN((Q)BDY(c->nm)->c);
1.1 noro 945: }
1.9 noro 946:
947: /* convert da to a univariate poly; return the position of variable */
948:
949: int dalgtoup(DAlg da,P *up,Q *dn)
950: {
951: int nv,i,hi,current_d;
952: DCP dc0,dc;
953: MP h,mp0,mp,t;
954: DL hd,d;
955: DP c;
956: DAlg cc;
957: P v;
958:
959: nv = da->nm->nv;
960: h = BDY(da->nm);
961: *dn = da->dn;
962: hd = h->dl;
963: for ( i = 0; i < nv; i++ )
964: if ( hd->d[i] ) break;
965: hi = i;
966: current_d = hd->d[i];
967: dc0 = 0;
968: mp0 = 0;
969: for ( t = h; t; t = NEXT(t) ) {
970: NEWDL(d,nv);
971: for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
972: for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
973: d->td = t->dl->td - t->dl->d[hi];
974: if ( t->dl->d[hi] != current_d ) {
975: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
976: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
977: current_d = t->dl->d[hi];
978: mp0 = 0;
979: }
980: NEXTMP(mp0,mp);
981: mp->c = t->c; mp->dl = d;
982: }
983: NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
984: NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
985: NEXT(dc) = 0;
986: makevar("x",&v);
987: MKP(VR(v),dc0,*up);
988: return hi;
989: }
990: