Annotation of OpenXM_contrib2/asir2018/engine/dalg.c, Revision 1.1
1.1 ! noro 1: /*
! 2: * $OpenXM$
! 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;
! 49: STOQ(2,two);
! 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]);
! 55: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
! 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]);
! 98: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
! 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);
! 437: divz(dna,gn,&a1); divz(dnb,gn,&b1);
! 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);
! 471: divz(dna,gn,&a1); divz(dnb,gn,&b1);
! 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);
! 540: divz((Z)BDY(a->nm)->c,(Z)BDY(u)->c,&cont);
! 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: }
! 599: gcdz(simp[i]->dn,ln,&gn); divz(ln,gn,&qn);
! 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;
! 670: dp_ptozp(a->nm,&u); divz((Z)BDY(a->nm)->c,(Z)BDY(u)->c,&nmc);
! 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] ) {
! 693: gcdz(simp[i]->dn,ln,&gn); divz(ln,gn,&qn);
! 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);
! 786: STOQ(i,iq); BDY(indt) = iq;
! 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] ) {
! 821: gcdz(simp[i]->dn,ln,&gn); divz(ln,gn,&qn);
! 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;
! 925: STOQ(2,two);
! 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);
! 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:
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