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