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