![[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.3 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.2 2004/12/02 08:39:54 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.3 ! noro 11: void invdalg(DAlg a,DAlg *c);
! 12: void rmcontdalg(DAlg a, DAlg *c);
1.1 noro 13:
14: void setfield_dalg(NODE alist)
15: {
16: NumberField nf;
17: VL vl,vl1,vl2;
18: int n,i,dim;
19: Alg *gen;
20: P *defpoly;
21: P p;
22: Q c,iq,two;
23: DP *ps,*mb;
1.3 ! noro 24: DP one;
1.1 noro 25: NODE t,b,b1,b2,hlist,mblist;
26: struct order_spec *current_spec;
27:
28: nf = (NumberField)MALLOC(sizeof(struct oNumberField));
29: current_numberfield = nf;
30: vl = 0;
31: for ( t = alist; t; t = NEXT(t) ) {
32: clctalg(BDY((Alg)BDY(t)),&vl1);
33: mergev(ALG,vl,vl1,&vl2); vl = vl2;
34: }
35: for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );
36: nf->n = n;
37: nf->vl = vl;
38: nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));
39: nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));
40: current_spec = dp_current_spec;
41: STOQ(2,two);
42: create_order_spec(0,(Obj)two,&nf->spec);
43: initd(nf->spec);
44: for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {
45: ptozp(vl1->v->attr,1,&c,&defpoly[i]);
46: ptod(ALG,vl,defpoly[i],&ps[i]);
47: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
1.3 ! noro 48: MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
1.1 noro 49: }
1.3 ! noro 50: ptod(ALG,vl,(P)ONE,&one);
! 51: MKDAlg(one,ONE,nf->one);
! 52: nf->ind = b;
! 53: dp_mbase(hlist,&mblist);
1.1 noro 54: initd(current_spec);
55: nf->dim = dim = length(mblist);
56: nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
57: for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
1.3 ! noro 58: mb[i] = (DP)BDY(t);
1.1 noro 59: }
60:
61: void algtodalg(Alg a,DAlg *r)
62: {
63: P ap,p,p1;
1.3 ! noro 64: Q c,c1,d1,dn,nm;
1.1 noro 65: DP dp;
66: DAlg da;
67: NumberField nf;
68: struct order_spec *current_spec;
1.3 ! noro 69: VL vl,tvl,svl;
! 70: V v;
1.1 noro 71:
72: if ( !(nf=current_numberfield) )
73: error("algtodalg : current_numberfield is not set");
1.3 ! noro 74: if ( !a ) {
! 75: *r = 0;
! 76: return;
! 77: }
! 78: switch (NID((Num)a) ) {
! 79: case N_Q:
! 80: c = (Q)a;
! 81: if ( INT(c) ) {
! 82: muldc(CO,nf->one->nm,(P)c,&dp);
! 83: MKDAlg(dp,ONE,*r);
! 84: } else {
! 85: NTOQ(NM(c),SGN(c),c1);
! 86: NTOQ(DN(c),1,d1);
! 87: muldc(CO,nf->one->nm,(P)c1,&dp);
! 88: MKDAlg(dp,c1,*r);
! 89: }
! 90: break;
! 91: case N_A:
! 92: ap = (P)BDY(a);
! 93: ptozp(ap,1,&c,&p);
! 94: if ( INT(c) ) {
! 95: p = ap;
! 96: dn = ONE;
! 97: } else {
! 98: NTOQ(NM(c),SGN(c),nm);
! 99: NTOQ(DN(c),1,dn);
! 100: mulpq(p,(P)nm,&p1); p = p1;
! 101: }
! 102: current_spec = dp_current_spec; initd(nf->spec);
! 103: get_vars(p,&vl);
! 104: for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
! 105: v = tvl->v;
! 106: for ( svl = nf->vl; svl; svl = NEXT(svl) )
! 107: if ( v == svl->v )
! 108: break;
! 109: if ( !svl )
! 110: error("algtodalg : incompatible numberfield");
! 111: }
! 112: ptod(ALG,nf->vl,p,&dp);
! 113: MKDAlg(dp,dn,da);
! 114: simpdalg(da,r);
! 115: break;
! 116: default:
! 117: error("algtodalg : invalid argument");
! 118: break;
1.1 noro 119: }
120: }
121:
1.2 noro 122: void dalgtoalg(DAlg da,Alg *r)
1.1 noro 123: {
1.2 noro 124: NumberField nf;
125: P p,p1;
126: Q inv;
127:
128: if ( !(nf=current_numberfield) )
1.3 ! noro 129: error("dalgtoalg : current_numberfield is not set");
1.2 noro 130: dtop(ALG,nf->vl,da->nm,&p);
131: invq(da->dn,&inv);
132: mulpq(p,(P)inv,&p1);
133: MKAlg(p1,*r);
1.1 noro 134: }
135:
136: void simpdalg(DAlg da,DAlg *r)
137: {
1.2 noro 138: NumberField nf;
139: DP nm;
1.3 ! noro 140: DAlg d;
1.2 noro 141: Q dn,dn1;
142:
143: if ( !(nf=current_numberfield) )
1.3 ! noro 144: error("simpdalg : current_numberfield is not set");
! 145: if ( !da ) {
! 146: *r = 0;
! 147: return;
! 148: }
1.2 noro 149: dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
150: mulq(da->dn,dn,&dn1);
1.3 ! noro 151: MKDAlg(nm,dn1,d);
! 152: rmcontdalg(d,r);
1.1 noro 153: }
154:
155: void adddalg(DAlg a,DAlg b,DAlg *c)
156: {
1.3 ! noro 157: NumberField nf;
! 158: Q dna,dnb,a1,b1,dn,g;
! 159: N an,bn,gn;
! 160: DP ta,tb,nm;
! 161: struct order_spec *current_spec;
! 162:
! 163: if ( !(nf=current_numberfield) )
! 164: error("adddalg : current_numberfield is not set");
! 165: if ( !a )
! 166: *c = b;
! 167: else if ( !b )
! 168: *c = a;
! 169: else {
! 170: dna = a->dn;
! 171: dnb = b->dn;
! 172: gcdn(NM(dna),NM(dnb),&gn);
! 173: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
! 174: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
! 175: /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
! 176: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
! 177: current_spec = dp_current_spec; initd(nf->spec);
! 178: addd(CO,ta,tb,&nm);
! 179: initd(current_spec);
! 180: if ( !nm )
! 181: *c = 0;
! 182: else {
! 183: mulq(dna,b1,&dn);
! 184: MKDAlg(nm,dn,*c);
! 185: }
! 186: }
1.1 noro 187: }
188:
189: void subdalg(DAlg a,DAlg b,DAlg *c)
190: {
1.3 ! noro 191: NumberField nf;
! 192: Q dna,dnb,a1,b1,dn,g;
! 193: N an,bn,gn;
! 194: DP ta,tb,nm;
! 195: struct order_spec *current_spec;
! 196:
! 197: if ( !(nf=current_numberfield) )
! 198: error("subdalg : current_numberfield is not set");
! 199: if ( !a )
! 200: *c = b;
! 201: else if ( !b )
! 202: *c = a;
! 203: else {
! 204: dna = a->dn;
! 205: dnb = b->dn;
! 206: gcdn(NM(dna),NM(dnb),&gn);
! 207: divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
! 208: NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
! 209: /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
! 210: muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
! 211: current_spec = dp_current_spec; initd(nf->spec);
! 212: subd(CO,ta,tb,&nm);
! 213: initd(current_spec);
! 214: if ( !nm )
! 215: *c = 0;
! 216: else {
! 217: mulq(dna,b1,&dn);
! 218: MKDAlg(nm,dn,*c);
! 219: }
! 220: }
1.1 noro 221: }
222:
1.2 noro 223: void muldalg(DAlg a,DAlg b,DAlg *c)
224: {
1.3 ! noro 225: NumberField nf;
! 226: DP nm;
! 227: Q dn;
! 228: DAlg t;
! 229: struct order_spec *current_spec;
! 230:
! 231: if ( !(nf=current_numberfield) )
! 232: error("muldalg : current_numberfield is not set");
! 233: if ( !a || !b )
! 234: *c = 0;
! 235: else {
! 236: current_spec = dp_current_spec; initd(nf->spec);
! 237: muld(CO,a->nm,b->nm,&nm);
! 238: initd(current_spec);
! 239: mulq(a->dn,b->dn,&dn);
! 240: MKDAlg(nm,dn,t);
! 241: simpdalg(t,c);
! 242: }
1.2 noro 243: }
244:
245:
1.1 noro 246: void divdalg(DAlg a,DAlg b,DAlg *c)
247: {
1.3 ! noro 248: DAlg inv;
! 249:
1.1 noro 250: if ( !current_numberfield )
1.3 ! noro 251: error("divdalg : current_numberfield is not set");
! 252: if ( !b )
! 253: error("divdalg : division by 0");
! 254: if ( !a )
! 255: c = 0;
! 256: else {
! 257: invdalg(b,&inv);
! 258: muldalg(a,inv,c);
! 259: }
! 260: }
! 261:
! 262: void rmcontdalg(DAlg a, DAlg *r)
! 263: {
! 264: DP u,u1;
! 265: Q cont,c,d;
! 266: N gn,cn,dn;
! 267:
! 268: if ( !a )
! 269: *r = a;
! 270: else {
! 271: dp_ptozp(a->nm,&u);
! 272: divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
! 273: gcdn(NM(cont),NM(a->dn),&gn);
! 274: divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
! 275: divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
! 276: muldc(CO,u,(P)c,&u1);
! 277: MKDAlg(u1,d,*r);
! 278: }
1.1 noro 279: }
280:
281: void invdalg(DAlg a,DAlg *c)
282: {
1.3 ! noro 283: NumberField nf;
! 284: int dim,n,i,j;
! 285: DP *mb;
! 286: DP m,d,u;
! 287: N ln,gn,qn;
! 288: DAlg *simp;
! 289: DAlg t,a0,r;
! 290: Q dn,dnsol,mul;
! 291: MAT mobj,sol;
! 292: Q **mat,**solmat;
! 293: MP mp0,mp;
! 294: int *rinfo,*cinfo;
! 295: struct order_spec *current_spec;
! 296:
! 297: if ( !(nf=current_numberfield) )
! 298: error("invdalg : current_numberfield is not set");
! 299: if ( !a )
! 300: error("invdalg : division by 0");
! 301: dim = nf->dim;
! 302: mb = nf->mb;
! 303: n = nf->n;
! 304: ln = ONEN;
! 305: MKDAlg(a->nm,ONE,a0);
! 306: simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
! 307: current_spec = dp_current_spec; initd(nf->spec);
! 308: for ( i = dim-1; i >= 0; i-- ) {
! 309: m = mb[i];
! 310: for ( j = i+1; j < dim; j++ )
! 311: if ( dp_redble(m,mb[j]) )
! 312: break;
! 313: if ( j < dim ) {
! 314: dp_subd(m,mb[j],&d);
! 315: muld(CO,d,simp[j]->nm,&u);
! 316: MKDAlg(u,simp[j]->dn,t);
! 317: simpdalg(t,&simp[i]);
! 318: } else {
! 319: MKDAlg(m,ONE,t);
! 320: muldalg(t,a0,&simp[i]);
! 321: }
! 322: gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
! 323: muln(NM(simp[i]->dn),qn,&ln);
! 324: }
! 325: initd(current_spec);
! 326: NTOQ(ln,1,dn);
! 327: MKMAT(mobj,dim,dim+1);
! 328: mat = (Q **)BDY(mobj);
! 329: mulq(dn,a->dn,&mat[dim-1][dim]);
! 330: for ( j = 0; j < dim; j++ ) {
! 331: divq(dn,simp[j]->dn,&mul);
! 332: for ( i = 0, mp = BDY(simp[j]->nm); mp && i < dim; i++ )
! 333: if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
! 334: mulq(mul,(Q)mp->c,&mat[i][j]);
! 335: mp = NEXT(mp);
! 336: }
! 337: }
! 338: generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
! 339: solmat = (Q **)BDY(sol);
! 340: for ( i = 0, mp0 = 0; i < dim; i++ )
! 341: if ( solmat[i][0] ) {
! 342: NEXTMP(mp0,mp);
! 343: mp->c = (P)solmat[i][0];
! 344: mp->dl = BDY(mb[i])->dl;
! 345: }
! 346: NEXT(mp) = 0; MKDP(n,mp0,u);
! 347: MKDAlg(u,dnsol,r);
! 348: rmcontdalg(r,c);
1.1 noro 349: }
350:
351: void chsgndalg(DAlg a,DAlg *c)
352: {
1.3 ! noro 353: DP nm;
! 354:
! 355: if ( !a ) *c = 0;
! 356: else {
! 357: chsgnd(a->nm,&nm);
! 358: MKDAlg(nm,a->dn,*c);
! 359: }
1.1 noro 360: }
361:
1.3 ! noro 362: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1 noro 363: {
1.3 ! noro 364: NumberField nf;
! 365: DAlg t,z,y;
! 366: N en,qn;
! 367: int r;
! 368:
! 369: if ( !(nf=current_numberfield) )
! 370: error("pwrdalg : current_numberfield is not set");
! 371: if ( !e )
! 372: *c = nf->one;
! 373: else if ( UNIQ(e) )
! 374: *c = a;
! 375: else {
! 376: if ( SGN(e) < 0 ) {
! 377: invdalg(a,&t); a = t;
! 378: }
! 379: en = NM(e);
! 380: y = nf->one;
! 381: z = a;
! 382: while ( 1 ) {
! 383: r = divin(en,2,&qn); en = qn;
! 384: if ( r ) {
! 385: muldalg(z,y,&t); y = t;
! 386: if ( !en ) {
! 387: *c = y;
! 388: return;
! 389: }
! 390: }
! 391: muldalg(z,z,&t); z = t;
! 392: }
! 393: }
1.1 noro 394: }
395:
1.3 ! noro 396: int cmpdalg(DAlg a,DAlg b)
1.1 noro 397: {
1.3 ! noro 398: DAlg c;
! 399:
! 400: subdalg(a,b,&c);
! 401: if ( !c ) return 0;
! 402: else
! 403: return SGN((Q)BDY(c->nm)->c);
1.1 noro 404: }