Annotation of OpenXM_contrib2/asir2000/engine/dalg.c, Revision 1.18
1.1 noro 1: /*
1.18 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.17 2017/08/31 02:36:21 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: {
1.18 ! noro 18: return current_numberfield;
1.5 noro 19: }
20:
1.1 noro 21: void setfield_dalg(NODE alist)
22: {
1.18 ! noro 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;
! 31: DP one;
! 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((P)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;
! 55: MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
! 56: }
! 57: ptod(ALG,vl,(P)ONE,&one);
! 58: MKDAlg(one,ONE,nf->one);
! 59: nf->ind = b;
! 60: dp_mbase(hlist,&mblist);
! 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++ )
! 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: {
1.18 ! noro 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);
1.10 noro 108: }
109:
1.4 noro 110: void qtodalg(Q q,DAlg *r)
111: {
1.18 ! noro 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,(Obj)q,&nm);
! 125: MKDAlg(nm,ONE,*r);
! 126: } else {
! 127: NTOQ(NM(q),SGN(q),t);
! 128: muldc(CO,nf->one->nm,(Obj)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: {
1.18 ! noro 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: }
1.6 noro 222: }
223:
224: void obj_dalgtoalg(Obj obj,Obj *r)
225: {
1.18 ! noro 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: {
1.18 ! noro 314: P ap,p,p1;
! 315: Q c,c1,d1,dn,nm;
! 316: DP dp;
! 317: DAlg da;
! 318: NumberField nf;
! 319: struct order_spec *current_spec;
! 320: VL vl,tvl,svl;
! 321: V v;
! 322:
! 323: if ( !(nf=current_numberfield) )
! 324: error("algtodalg : current_numberfield is not set");
! 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,(Obj)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,(Obj)c1,&dp);
! 339: MKDAlg(dp,d1,*r);
! 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((Obj)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;
! 370: }
1.1 noro 371: }
372:
1.2 noro 373: void dalgtoalg(DAlg da,Alg *r)
1.1 noro 374: {
1.18 ! noro 375: NumberField nf;
! 376: P p,p1;
! 377: Q inv;
! 378:
! 379: if ( !(nf=current_numberfield) )
! 380: error("dalgtoalg : current_numberfield is not set");
! 381: if ( !da ) *r = 0;
! 382: else {
! 383: dtop(ALG,nf->vl,da->nm,(Obj *)&p);
! 384: invq(da->dn,&inv);
! 385: mulpq(p,(P)inv,&p1);
! 386: MKAlg(p1,*r);
! 387: }
1.1 noro 388: }
389:
390: void simpdalg(DAlg da,DAlg *r)
391: {
1.18 ! noro 392: NumberField nf;
! 393: DP nm;
! 394: DAlg d;
! 395: Q dn,dn1;
! 396: struct order_spec *current_spec;
! 397:
! 398: if ( !(nf=current_numberfield) )
! 399: error("simpdalg : current_numberfield is not set");
! 400: if ( !da ) {
! 401: *r = 0;
! 402: return;
! 403: }
! 404: current_spec = dp_current_spec; initd(nf->spec);
! 405: dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,(P *)&dn);
! 406: if ( !nm ) *r = 0;
! 407: else {
! 408: initd(current_spec);
! 409: mulq(da->dn,dn,&dn1);
! 410: MKDAlg(nm,dn1,d);
! 411: rmcontdalg(d,r);
! 412: }
1.1 noro 413: }
414:
415: void adddalg(DAlg a,DAlg b,DAlg *c)
416: {
1.18 ! noro 417: NumberField nf;
! 418: Q dna,dnb,a1,b1,dn,g;
! 419: N an,bn,gn;
! 420: DAlg t;
! 421: DP ta,tb,nm;
! 422: struct order_spec *current_spec;
! 423:
! 424: if ( !(nf=current_numberfield) )
! 425: error("adddalg : current_numberfield is not set");
! 426: if ( !a )
! 427: *c = b;
! 428: else if ( !b )
! 429: *c = a;
! 430: else {
! 431: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
! 432: dna = a->dn;
! 433: dnb = b->dn;
! 434: gcdn(NM(dna),NM(dnb),&gn);
! 435: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
! 436: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
! 437: /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
! 438: muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
! 439: current_spec = dp_current_spec; initd(nf->spec);
! 440: addd(CO,ta,tb,&nm);
! 441: initd(current_spec);
! 442: if ( !nm )
! 443: *c = 0;
! 444: else {
! 445: mulq(dna,b1,&dn);
! 446: MKDAlg(nm,dn,*c);
! 447: }
! 448: }
1.1 noro 449: }
450:
451: void subdalg(DAlg a,DAlg b,DAlg *c)
452: {
1.18 ! noro 453: NumberField nf;
! 454: Q dna,dnb,a1,b1,dn,g;
! 455: N an,bn,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: gcdn(NM(dna),NM(dnb),&gn);
! 471: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
! 472: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
! 473: /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
! 474: muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
! 475: current_spec = dp_current_spec; initd(nf->spec);
! 476: subd(CO,ta,tb,&nm);
! 477: initd(current_spec);
! 478: if ( !nm )
! 479: *c = 0;
! 480: else {
! 481: mulq(dna,b1,&dn);
! 482: MKDAlg(nm,dn,*c);
! 483: }
! 484: }
1.1 noro 485: }
486:
1.2 noro 487: void muldalg(DAlg a,DAlg b,DAlg *c)
488: {
1.18 ! noro 489: NumberField nf;
! 490: DP nm;
! 491: Q dn;
! 492: DAlg t;
! 493: struct order_spec *current_spec;
! 494:
! 495: if ( !(nf=current_numberfield) )
! 496: error("muldalg : current_numberfield is not set");
! 497: if ( !a || !b )
! 498: *c = 0;
! 499: else {
! 500: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
! 501: current_spec = dp_current_spec; initd(nf->spec);
! 502: muld(CO,a->nm,b->nm,&nm);
! 503: initd(current_spec);
! 504: mulq(a->dn,b->dn,&dn);
! 505: MKDAlg(nm,dn,t);
! 506: simpdalg(t,c);
! 507: }
1.2 noro 508: }
509:
510:
1.1 noro 511: void divdalg(DAlg a,DAlg b,DAlg *c)
512: {
1.18 ! noro 513: DAlg inv,t;
! 514: int ret;
1.3 noro 515:
1.18 ! noro 516: if ( !current_numberfield )
! 517: error("divdalg : current_numberfield is not set");
! 518: if ( !b )
! 519: error("divdalg : division by 0");
! 520: if ( !a )
! 521: c = 0;
! 522: else {
! 523: qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
! 524: ret = invdalg(b,&inv);
! 525: if ( !ret ) {
! 526: error("divdalg : the denominator is not invertible");
! 527: }
! 528: muldalg(a,inv,c);
! 529: }
1.3 noro 530: }
531:
532: void rmcontdalg(DAlg a, DAlg *r)
533: {
1.18 ! noro 534: DP u,u1;
! 535: Q cont,c,d;
! 536: N gn,cn,dn;
! 537:
! 538: if ( !a )
! 539: *r = a;
! 540: else {
! 541: dp_ptozp(a->nm,&u);
! 542: divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
! 543: gcdn(NM(cont),NM(a->dn),&gn);
! 544: divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
! 545: divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
! 546: muldc(CO,u,(Obj)c,&u1);
! 547: MKDAlg(u1,d,*r);
! 548: }
1.1 noro 549: }
550:
1.8 noro 551: int invdalg(DAlg a,DAlg *c)
1.1 noro 552: {
1.18 ! noro 553: NumberField nf;
! 554: int dim,n,i,j,k,l;
! 555: DP *mb;
! 556: DP m,d,u;
! 557: N ln,gn,qn;
! 558: DAlg *simp;
! 559: DAlg t,a0,r;
! 560: Q dn,dnsol,mul,nmc,dn1;
! 561: MAT mobj,sol;
! 562: Q **mat,**solmat;
! 563: MP mp0,mp;
! 564: int *rinfo,*cinfo;
! 565: int rank,nparam;
! 566: NODE nd0,nd,ndt;
! 567: struct order_spec *current_spec;
! 568: struct oEGT eg0,eg1;
! 569: extern struct oEGT eg_le;
! 570:
! 571: if ( !(nf=current_numberfield) )
! 572: error("invdalg : current_numberfield is not set");
! 573: if ( !a )
! 574: error("invdalg : division by 0");
! 575: else if ( NID(a) == N_Q ) {
! 576: invq((Q)a,&dn); *c = (DAlg)dn;
! 577: return 1;
! 578: }
! 579: dim = nf->dim;
! 580: mb = nf->mb;
! 581: n = nf->n;
! 582: ln = ONEN;
! 583: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
! 584: MKDAlg(u,ONE,a0);
! 585: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
! 586: current_spec = dp_current_spec; initd(nf->spec);
! 587: for ( i = 0; i < dim; i++ ) {
! 588: m = mb[i];
! 589: for ( j = i-1; j >= 0; j-- )
! 590: if ( dp_redble(m,mb[j]) )
! 591: break;
! 592: if ( j >= 0 ) {
! 593: dp_subd(m,mb[j],&d);
! 594: muld(CO,d,simp[j]->nm,&u);
! 595: MKDAlg(u,simp[j]->dn,t);
! 596: simpdalg(t,&simp[i]);
! 597: } else {
! 598: MKDAlg(m,ONE,t);
! 599: muldalg(t,a0,&simp[i]);
! 600: }
! 601: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
! 602: muln(NM(simp[i]->dn),qn,&ln);
! 603: }
! 604: initd(current_spec);
! 605: NTOQ(ln,1,dn);
! 606: MKMAT(mobj,dim,dim+1);
! 607: mat = (Q **)BDY(mobj);
! 608: mulq(dn,a->dn,&mat[0][dim]);
! 609: for ( j = 0; j < dim; j++ ) {
! 610: divq(dn,simp[j]->dn,&mul);
! 611: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
! 612: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
! 613: mulq(mul,(Q)mp->c,&mat[i][j]);
! 614: mp = NEXT(mp);
! 615: }
! 616: }
! 617: get_eg(&eg0);
! 618: rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
! 619: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
! 620: if ( cinfo[0] == dim ) {
! 621: /* the input is invertible */
! 622: solmat = (Q **)BDY(sol);
! 623: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
! 624: if ( solmat[i][0] ) {
! 625: NEXTMP(mp0,mp);
! 626: mp->c = (Obj)solmat[i][0];
! 627: mp->dl = BDY(mb[i])->dl;
! 628: }
! 629: NEXT(mp) = 0; MKDP(n,mp0,u);
! 630: mulq(dnsol,nmc,&dn1);
! 631: MKDAlg(u,dn1,r);
! 632: rmcontdalg(r,c);
! 633: return 1;
! 634: } else
! 635: return 0;
1.8 noro 636: }
637:
638: NODE inv_or_split_dalg(DAlg a,DAlg *c)
639: {
1.18 ! noro 640: NumberField nf;
! 641: int dim,n,i,j,k,l;
! 642: DP *mb;
! 643: DP m,d,u;
! 644: N ln,gn,qn;
! 645: DAlg *simp;
! 646: DAlg t,a0,r;
! 647: Q dn,dnsol,mul,nmc,dn1;
! 648: MAT mobj,sol;
! 649: Q **mat,**solmat;
! 650: MP mp0,mp;
! 651: int *rinfo,*cinfo;
! 652: int rank,nparam;
! 653: NODE nd0,nd,ndt;
! 654: struct order_spec *current_spec;
! 655: struct oEGT eg0,eg1;
! 656: extern struct oEGT eg_le;
! 657: extern int DP_Print;
! 658:
! 659: if ( !(nf=current_numberfield) )
! 660: error("invdalg : current_numberfield is not set");
! 661: if ( !a )
! 662: error("invdalg : division by 0");
! 663: else if ( NID(a) == N_Q ) {
! 664: invq((Q)a,&dn); *c = (DAlg)dn;
! 665: return 0;
! 666: }
! 667: dim = nf->dim;
! 668: mb = nf->mb;
! 669: n = nf->n;
! 670: ln = ONEN;
! 671: dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
! 672: MKDAlg(u,ONE,a0);
! 673: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
! 674: current_spec = dp_current_spec; initd(nf->spec);
! 675: for ( i = 0; i < dim; i++ ) {
! 676: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
! 677: m = mb[i];
! 678: for ( j = i-1; j >= 0; j-- )
! 679: if ( dp_redble(m,mb[j]) )
! 680: break;
! 681: if ( j >= 0 ) {
! 682: dp_subd(m,mb[j],&d);
! 683: if ( simp[j] ) {
! 684: muld(CO,d,simp[j]->nm,&u);
! 685: MKDAlg(u,simp[j]->dn,t);
! 686: simpdalg(t,&simp[i]);
! 687: } else
! 688: simp[i] = 0;
! 689: } else {
! 690: MKDAlg(m,ONE,t);
! 691: muldalg(t,a0,&simp[i]);
! 692: }
! 693: if ( simp[i] ) {
! 694: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
! 695: muln(NM(simp[i]->dn),qn,&ln);
! 696: }
! 697: }
! 698: initd(current_spec);
! 699: NTOQ(ln,1,dn);
! 700: MKMAT(mobj,dim,dim+1);
! 701: mat = (Q **)BDY(mobj);
! 702: mulq(dn,a->dn,&mat[0][dim]);
! 703: for ( j = 0; j < dim; j++ ) {
! 704: if ( simp[j] ) {
! 705: divq(dn,simp[j]->dn,&mul);
! 706: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
! 707: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
! 708: mulq(mul,(Q)mp->c,&mat[i][j]);
! 709: mp = NEXT(mp);
! 710: }
! 711: }
! 712: }
! 713: get_eg(&eg0);
! 714: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
! 715: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
! 716: if ( cinfo[0] == dim ) {
! 717: /* the input is invertible */
! 718: solmat = (Q **)BDY(sol);
! 719: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
! 720: if ( solmat[i][0] ) {
! 721: NEXTMP(mp0,mp);
! 722: mp->c = (Obj)solmat[i][0];
! 723: mp->dl = BDY(mb[i])->dl;
! 724: }
! 725: NEXT(mp) = 0; MKDP(n,mp0,u);
! 726: mulq(dnsol,nmc,&dn1);
! 727: MKDAlg(u,dn1,r);
! 728: rmcontdalg(r,c);
! 729: return 0;
! 730: } else {
! 731: /* the input is not invertible */
! 732: nparam = sol->col;
! 733: solmat = (Q **)BDY(sol);
! 734: nd0 = 0;
! 735: for ( k = 0; k < nparam; k++ ) {
! 736: /* construct a new basis element */
! 737: m = mb[cinfo[k]];
! 738: mp0 = 0;
! 739: NEXTMP(mp0,mp);
! 740: chsgnq(dnsol,&dn1); mp->c = (Obj)dn1;
! 741: mp->dl = BDY(m)->dl;
! 742: /* skip the last parameter */
! 743: for ( l = rank-2; l >= 0; l-- ) {
! 744: if ( solmat[l][k] ) {
! 745: NEXTMP(mp0,mp);
! 746: mp->c = (Obj)solmat[l][k];
! 747: mp->dl = BDY(mb[rinfo[l]])->dl;
! 748: }
! 749: }
! 750: NEXT(mp) = 0; MKDP(n,mp0,u);
! 751: NEXTNODE(nd0,nd);
! 752: BDY(nd) = (pointer)u;
! 753: NEXT(nd) = 0;
! 754: }
! 755: NEXT(nd) = 0;
! 756: return nd0;
! 757: }
1.1 noro 758: }
759:
1.14 noro 760: NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv)
761: {
1.18 ! noro 762: int dim,n,i,j,k,l,nv;
! 763: DP *mb,*ps;
! 764: DP m,d,u,nm;
! 765: N ln,gn,qn;
! 766: DAlg *simp;
! 767: DAlg a0,r;
! 768: Q dn,dnsol,mul,nmc,dn1,iq;
! 769: MAT mobj,sol;
! 770: Q **mat,**solmat;
! 771: MP mp0,mp;
! 772: int *rinfo,*cinfo;
! 773: int rank,nparam;
! 774: NODE nd0,nd,ndt,ind,indt,t,mblist;
! 775: struct oEGT eg0,eg1;
! 776: extern struct oEGT eg_le;
! 777: extern int DP_Print;
! 778: initd(spec);
! 779: dp_ptozp(f,&u); f = u;
! 780:
! 781: n = length(gb);
! 782: ps = (DP *)MALLOC(n*sizeof(DP));
! 783: for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) {
! 784: ps[i] = (DP)BDY(t);
! 785: NEXTNODE(ind,indt);
! 786: STOQ(i,iq); BDY(indt) = iq;
! 787: }
! 788: if ( ind ) NEXT(indt) = 0;
! 789: dp_true_nf(ind,f,ps,1,&nm,(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 = ONEN;
! 800: simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
! 801: for ( i = 0; i < dim; i++ ) {
! 802: if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
! 803: m = mb[i];
! 804: for ( j = i-1; j >= 0; j-- )
! 805: if ( dp_redble(m,mb[j]) )
! 806: break;
! 807: if ( j >= 0 ) {
! 808: dp_subd(m,mb[j],&d);
! 809: if ( simp[j] ) {
! 810: muld(CO,d,simp[j]->nm,&u);
! 811: dp_true_nf(ind,u,ps,1,&nm,(P *)&dn);
! 812: mulq(simp[j]->dn,dn,&dn1);
! 813: MKDAlg(nm,dn1,simp[i]);
! 814: } else
! 815: simp[i] = 0;
! 816: } else {
! 817: dp_true_nf(ind,f,ps,1,&nm,(P *)&dn);
! 818: MKDAlg(nm,dn,simp[i]);
! 819: }
! 820: if ( simp[i] ) {
! 821: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
! 822: muln(NM(simp[i]->dn),qn,&ln);
! 823: }
! 824: }
! 825: NTOQ(ln,1,dn);
! 826: MKMAT(mobj,dim,dim+1);
! 827: mat = (Q **)BDY(mobj);
! 828: mat[0][dim] = dn;
! 829: for ( j = 0; j < dim; j++ ) {
! 830: if ( simp[j] ) {
! 831: divq(dn,simp[j]->dn,&mul);
! 832: for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
! 833: if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) {
! 834: mulq(mul,(Q)mp->c,&mat[i][j]);
! 835: mp = NEXT(mp);
! 836: }
! 837: }
! 838: }
! 839: get_eg(&eg0);
! 840: rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
! 841: get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
! 842: if ( cinfo[0] == dim ) {
! 843: /* the input is invertible */
! 844: solmat = (Q **)BDY(sol);
! 845: for ( i = dim-1, mp0 = 0; i >= 0; i-- )
! 846: if ( solmat[i][0] ) {
! 847: NEXTMP(mp0,mp);
! 848: mp->c = (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 = (Q **)BDY(sol);
! 857: nd0 = 0;
! 858: for ( k = 0; k < nparam; k++ ) {
! 859: /* construct a new basis element */
! 860: m = mb[cinfo[k]];
! 861: mp0 = 0;
! 862: NEXTMP(mp0,mp);
! 863: chsgnq(dnsol,&dn1); mp->c = (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: }
1.14 noro 881: }
882:
1.1 noro 883: void chsgndalg(DAlg a,DAlg *c)
884: {
1.18 ! noro 885: DP nm;
! 886: Q t;
1.3 noro 887:
1.18 ! noro 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: }
1.1 noro 895: }
896:
1.3 noro 897: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 898: {
1.18 ! noro 899: NumberField nf;
! 900: DAlg t,z,y;
! 901: Q q;
! 902: N en,qn;
! 903: int r;
! 904: int ret;
! 905:
! 906: if ( !(nf=current_numberfield) )
! 907: error("pwrdalg : current_numberfield is not set");
! 908: if ( !a )
! 909: *c = !e ? (DAlg)ONE : 0;
! 910: else if ( NID(a) == N_Q ) {
! 911: pwrq((Q)a,e,&q); *c = (DAlg)q;
! 912: } else if ( !e )
! 913: *c = nf->one;
! 914: else if ( UNIQ(e) )
! 915: *c = a;
! 916: else {
! 917: if ( SGN(e) < 0 ) {
! 918: ret = invdalg(a,&t);
! 919: if ( !ret )
! 920: error("pwrdalg : the denominator is not invertible");
! 921: a = t;
! 922: }
! 923: en = NM(e);
! 924: y = nf->one;
! 925: z = a;
! 926: while ( 1 ) {
! 927: r = divin(en,2,&qn); en = qn;
! 928: if ( r ) {
! 929: muldalg(z,y,&t); y = t;
! 930: if ( !en ) {
! 931: *c = y;
! 932: return;
! 933: }
! 934: }
! 935: muldalg(z,z,&t); z = t;
! 936: }
! 937: }
1.1 noro 938: }
939:
1.3 noro 940: int cmpdalg(DAlg a,DAlg b)
1.1 noro 941: {
1.18 ! noro 942: DAlg c;
1.3 noro 943:
1.18 ! noro 944: subdalg(a,b,&c);
! 945: if ( !c ) return 0;
! 946: else
! 947: return SGN((Q)BDY(c->nm)->c);
1.1 noro 948: }
1.9 noro 949:
950: /* convert da to a univariate poly; return the position of variable */
951:
952: int dalgtoup(DAlg da,P *up,Q *dn)
953: {
1.18 ! noro 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;
1.9 noro 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 / asir2000 / engine