Annotation of OpenXM_contrib2/asir2018/builtin/algnum.c, Revision 1.1
1.1 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM$
! 49: */
! 50: #include "ca.h"
! 51: #include "parse.h"
! 52:
! 53: void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();
! 54: void Palg(), Palgv(), Pgetalgtree();
! 55: void Pinvalg_le();
! 56: void Pset_field(),Palgtodalg(),Pdalgtoalg();
! 57: void Pinv_or_split_dalg();
! 58: void Pdalgtoup();
! 59: void Pget_field_defpoly();
! 60: void Pget_field_generator();
! 61:
! 62: void mkalg(P,Alg *);
! 63: int cmpalgp(P,P);
! 64: void algptop(P,P *);
! 65: void algtorat(Num,Obj *);
! 66: void rattoalg(Obj,Alg *);
! 67: void ptoalgp(P,P *);
! 68: void clctalg(P,VL *);
! 69: void get_algtree(Obj f,VL *r);
! 70: void Pinvalg_chrem();
! 71: void Pdalgtodp();
! 72: void Pdptodalg();
! 73:
! 74: struct ftab alg_tab[] = {
! 75: {"set_field",Pset_field,-3},
! 76: {"get_field_defpoly",Pget_field_defpoly,1},
! 77: {"get_field_generator",Pget_field_generator,1},
! 78: {"algtodalg",Palgtodalg,1},
! 79: {"dalgtoalg",Pdalgtoalg,1},
! 80: {"dalgtodp",Pdalgtodp,1},
! 81: {"dalgtoup",Pdalgtoup,1},
! 82: {"dptodalg",Pdptodalg,1},
! 83: {"inv_or_split_dalg",Pinv_or_split_dalg,1},
! 84: {"invalg_chrem",Pinvalg_chrem,2},
! 85: {"invalg_le",Pinvalg_le,1},
! 86: {"defpoly",Pdefpoly,1},
! 87: {"newalg",Pnewalg,1},
! 88: {"mainalg",Pmainalg,1},
! 89: {"algtorat",Palgtorat,1},
! 90: {"rattoalg",Prattoalg,1},
! 91: {"getalg",Pgetalg,1},
! 92: {"getalgtree",Pgetalgtree,1},
! 93: {"alg",Palg,1},
! 94: {"algv",Palgv,1},
! 95: {0,0,0},
! 96: };
! 97:
! 98: static int UCN,ACNT;
! 99:
! 100: void Pset_field(NODE arg,Q *rp)
! 101: {
! 102: int ac;
! 103: NODE a0,a1;
! 104: VL vl0,vl;
! 105: struct order_spec *spec;
! 106:
! 107: if ( (ac = argc(arg)) == 1 )
! 108: setfield_dalg(BDY((LIST)ARG0(arg)));
! 109: else if ( ac == 3 ) {
! 110: a0 = BDY((LIST)ARG0(arg));
! 111: a1 = BDY((LIST)ARG1(arg));
! 112: for ( vl0 = 0; a1; a1 = NEXT(a1) ) {
! 113: NEXTVL(vl0,vl);
! 114: vl->v = VR((P)BDY(a1));
! 115: }
! 116: if ( vl0 ) NEXT(vl) = 0;
! 117: create_order_spec(0,ARG2(arg),&spec);
! 118: setfield_gb(a0,vl0,spec);
! 119: }
! 120: *rp = 0;
! 121: }
! 122:
! 123: void Palgtodalg(NODE arg,DAlg *rp)
! 124: {
! 125: algtodalg((Alg)ARG0(arg),rp);
! 126: }
! 127:
! 128: void Pdalgtoalg(NODE arg,Alg *rp)
! 129: {
! 130: dalgtoalg((DAlg)ARG0(arg),rp);
! 131: }
! 132:
! 133: void Pdalgtodp(NODE arg,LIST *r)
! 134: {
! 135: NODE b;
! 136: DP nm;
! 137: Z dn;
! 138: DAlg da;
! 139:
! 140: da = (DAlg)ARG0(arg);
! 141: nm = da->nm;
! 142: dn = da->dn;
! 143: b = mknode(2,nm,dn);
! 144: MKLIST(*r,b);
! 145: }
! 146:
! 147: void Pdptodalg(NODE arg,DAlg *r)
! 148: {
! 149: DP d,nm,nm1;
! 150: MP m;
! 151: Q c;
! 152: Z a;
! 153: DAlg t;
! 154:
! 155: d = (DP)ARG0(arg);
! 156: if ( !d ) *r = 0;
! 157: else {
! 158: for ( m = BDY(d); m; m = NEXT(m) )
! 159: if ( !INT((Q)m->c) ) break;
! 160: if ( !m ) {
! 161: MKDAlg(d,ONE,t);
! 162: } else {
! 163: dp_ptozp(d,&nm);
! 164: divq((Q)BDY(d)->c,(Q)BDY(nm)->c,&c);
! 165: nmq(c,&a);
! 166: muldc(CO,nm,(Obj)a,&nm1);
! 167: dnq(c,&a);
! 168: MKDAlg(nm1,a,t);
! 169: }
! 170: simpdalg(t,r);
! 171: }
! 172: }
! 173:
! 174: void Pdalgtoup(NODE arg,LIST *r)
! 175: {
! 176: NODE b;
! 177: int pos;
! 178: P up;
! 179: DP nm;
! 180: Z dn;
! 181: Z q;
! 182:
! 183: pos = dalgtoup((DAlg)ARG0(arg),&up,&dn);
! 184: STOQ(pos,q);
! 185: b = mknode(3,up,dn,q);
! 186: MKLIST(*r,b);
! 187: }
! 188:
! 189: NODE inv_or_split_dalg(DAlg,DAlg *);
! 190: NumberField get_numberfield();
! 191:
! 192: void Pget_field_defpoly(NODE arg,DAlg *r)
! 193: {
! 194: NumberField nf;
! 195: DP d;
! 196:
! 197: nf = get_numberfield();
! 198: d = nf->ps[QTOS((Q)ARG0(arg))];
! 199: MKDAlg(d,ONE,*r);
! 200: }
! 201:
! 202: void Pget_field_generator(NODE arg,DAlg *r)
! 203: {
! 204: int index,n,i;
! 205: DL dl;
! 206: MP m;
! 207: DP d;
! 208:
! 209: index = QTOS((Q)ARG0(arg));
! 210: n = get_numberfield()->n;
! 211: NEWDL(dl,n);
! 212: for ( i = 0; i < n; i++ ) dl->d[i] = 0;
! 213: dl->d[index] = 1; dl->td = 1;
! 214: NEWMP(m); m->dl = dl; m->c = (Obj)ONE; NEXT(m) = 0;
! 215: MKDP(n,m,d);
! 216: MKDAlg(d,ONE,*r);
! 217: }
! 218:
! 219:
! 220: void Pinv_or_split_dalg(NODE arg,Obj *rp)
! 221: {
! 222: NODE gen,t,nd0,nd;
! 223: LIST list;
! 224: int l,i,j,n;
! 225: DP *ps,*ps1,*psw;
! 226: NumberField nf;
! 227: DAlg inv;
! 228: extern struct order_spec *dp_current_spec;
! 229: struct order_spec *current_spec;
! 230:
! 231: gen = inv_or_split_dalg((DAlg)ARG0(arg),&inv);
! 232: if ( !gen )
! 233: *rp = (Obj)inv;
! 234: else {
! 235: nf = get_numberfield();
! 236: current_spec = dp_current_spec; initd(nf->spec);
! 237: l = length(gen);
! 238: n = nf->n;
! 239: ps = nf->ps;
! 240: psw = (DP *)ALLOCA((n+l)*sizeof(DP));
! 241: for ( i = j = 0; i < n; i++ ) {
! 242: for ( t = gen; t; t = NEXT(t) )
! 243: if ( dp_redble(ps[i],(DP)BDY(t)) ) break;
! 244: if ( !t )
! 245: psw[j++] = ps[i];
! 246: }
! 247: nd0 = 0;
! 248: /* gen[0] < gen[1] < ... */
! 249: /* psw[0] > psw[1] > ... */
! 250: for ( i = j-1, t = gen; i >= 0 && t; ) {
! 251: NEXTNODE(nd0,nd);
! 252: if ( compd(CO,psw[i],(DP)BDY(t)) > 0 ) {
! 253: BDY(nd) = BDY(t); t = NEXT(t);
! 254: } else
! 255: BDY(nd) = (pointer)psw[i--];
! 256: }
! 257: for ( ; i >= 0; i-- ) {
! 258: NEXTNODE(nd0,nd); BDY(nd) = (pointer)psw[i];
! 259: }
! 260: for ( ; t; t = NEXT(t) ) {
! 261: NEXTNODE(nd0,nd); BDY(nd) = BDY(t);
! 262: }
! 263: NEXT(nd) = 0;
! 264: MKLIST(list,nd0);
! 265: initd(current_spec);
! 266: *rp = (Obj)list;
! 267: }
! 268: }
! 269:
! 270: void Pnewalg(arg,rp)
! 271: NODE arg;
! 272: Alg *rp;
! 273: {
! 274: P p;
! 275: VL vl;
! 276: P c;
! 277:
! 278: p = (P)ARG0(arg);
! 279: if ( !p || OID(p) != O_P )
! 280: error("newalg : invalid argument");
! 281: clctv(CO,p,&vl);
! 282: if ( NEXT(vl) )
! 283: error("newalg : invalid argument");
! 284: c = COEF(DC(p));
! 285: if ( !NUM(c) || !RATN(c) )
! 286: error("newalg : invalid argument");
! 287: mkalg(p,rp);
! 288: }
! 289:
! 290: void mkalg(p,r)
! 291: P p;
! 292: Alg *r;
! 293: {
! 294: VL vl,mvl,nvl;
! 295: V a,tv;
! 296: char buf[BUFSIZ];
! 297: char *name;
! 298: P x,t,s;
! 299: Num c;
! 300: DCP dc,dcr,dcr0;
! 301:
! 302: for ( vl = ALG; vl; vl = NEXT(vl) )
! 303: if ( !cmpalgp(p,(P)vl->v->attr) ) {
! 304: a = vl->v; break;
! 305: }
! 306: if ( !vl ) {
! 307: NEWVL(vl); NEXT(vl) = ALG; ALG = vl;
! 308: NEWV(a); vl->v = a;
! 309: sprintf(buf,"#%d",ACNT++);
! 310: name = (char *)MALLOC(strlen(buf)+1);
! 311: strcpy(name,buf); NAME(a) = name;
! 312:
! 313: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 314: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); c = (Num)COEF(dc);
! 315: if ( NID(c) != N_A )
! 316: COEF(dcr) = (P)c;
! 317: else
! 318: COEF(dcr) = (P)BDY(((Alg)c));
! 319: }
! 320: NEXT(dcr) = 0; MKP(a,dcr0,t); a->attr = (pointer)t;
! 321:
! 322: sprintf(buf,"t%s",name); makevar(buf,&s);
! 323:
! 324: if ( NEXT(ALG) ) {
! 325: tv = (V)NEXT(ALG)->v->priv;
! 326: for ( vl = CO; NEXT(NEXT(vl)); vl = NEXT(vl) );
! 327: nvl = NEXT(vl); NEXT(vl) = 0;
! 328: for ( vl = CO; NEXT(vl) && (NEXT(vl)->v != tv); vl = NEXT(vl) );
! 329: mvl = NEXT(vl); NEXT(vl) = nvl; NEXT(nvl) = mvl;
! 330: }
! 331:
! 332: a->priv = (pointer)VR(s); VR(s)->priv = (pointer)a;
! 333: }
! 334: MKV(a,x); MKAlg(x,*r);
! 335: }
! 336:
! 337: int cmpalgp(p,defp)
! 338: P p,defp;
! 339: {
! 340: DCP dc,dcd;
! 341: P t;
! 342:
! 343: for ( dc = DC(p), dcd = DC(defp); dc && dcd;
! 344: dc = NEXT(dc), dcd = NEXT(dcd) ) {
! 345: if ( cmpz(DEG(dc),DEG(dcd)) )
! 346: break;
! 347: t = NID((Num)COEF(dc)) == N_A ? (P)BDY((Alg)COEF(dc)) : COEF(dc);
! 348: if ( compp(ALG,t,COEF(dcd)) )
! 349: break;
! 350: }
! 351: if ( dc || dcd )
! 352: return 1;
! 353: else
! 354: return 0;
! 355: }
! 356:
! 357: void Pdefpoly(arg,rp)
! 358: NODE arg;
! 359: P *rp;
! 360: {
! 361: asir_assert(ARG0(arg),O_N,"defpoly");
! 362: algptop((P)VR((P)BDY((Alg)ARG0(arg)))->attr,rp);
! 363: }
! 364:
! 365: void Pmainalg(arg,r)
! 366: NODE arg;
! 367: Alg *r;
! 368: {
! 369: Num c;
! 370: V v;
! 371: P b;
! 372:
! 373: c = (Num)(ARG0(arg));
! 374: if ( NID(c) <= N_R )
! 375: *r = 0;
! 376: else {
! 377: v = VR((P)BDY((Alg)c)); MKV(v,b); MKAlg(b,*r);
! 378: }
! 379: }
! 380:
! 381: void Palgtorat(arg,rp)
! 382: NODE arg;
! 383: Obj *rp;
! 384: {
! 385: asir_assert(ARG0(arg),O_N,"algtorat");
! 386: algtorat((Num)ARG0(arg),rp);
! 387: }
! 388:
! 389: void Prattoalg(arg,rp)
! 390: NODE arg;
! 391: Alg *rp;
! 392: {
! 393: asir_assert(ARG0(arg),O_R,"rattoalg");
! 394: rattoalg((Obj)ARG0(arg),rp);
! 395: }
! 396:
! 397: void Pgetalg(arg,rp)
! 398: NODE arg;
! 399: LIST *rp;
! 400: {
! 401: Obj t;
! 402: P p;
! 403: VL vl;
! 404: Num a;
! 405: Alg b;
! 406: NODE n0,n;
! 407:
! 408: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
! 409: vl = 0;
! 410: else {
! 411: t = BDY((Alg)a);
! 412: switch ( OID(t) ) {
! 413: case O_P: case O_R:
! 414: clctvr(ALG,t,&vl); break;
! 415: default:
! 416: vl = 0; break;
! 417: }
! 418: }
! 419: for ( n0 = 0; vl; vl = NEXT(vl) ) {
! 420: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
! 421: }
! 422: if ( n0 )
! 423: NEXT(n) = 0;
! 424: MKLIST(*rp,n0);
! 425: }
! 426:
! 427: void Pgetalgtree(arg,rp)
! 428: NODE arg;
! 429: LIST *rp;
! 430: {
! 431: Obj t;
! 432: P p;
! 433: VL vl,vl1,vl2;
! 434: Num a;
! 435: Alg b;
! 436: NODE n0,n;
! 437:
! 438: #if 0
! 439: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
! 440: vl = 0;
! 441: else {
! 442: t = BDY((Alg)a);
! 443: switch ( OID(t) ) {
! 444: case O_P:
! 445: clctalg((P)t,&vl); break;
! 446: case O_R:
! 447: clctalg(NM((R)t),&vl1);
! 448: clctalg(DN((R)t),&vl2);
! 449: mergev(ALG,vl1,vl2,&vl); break;
! 450: default:
! 451: vl = 0; break;
! 452: }
! 453: }
! 454: #else
! 455: get_algtree((Obj)ARG0(arg),&vl);
! 456: #endif
! 457: for ( n0 = 0; vl; vl = NEXT(vl) ) {
! 458: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
! 459: }
! 460: if ( n0 )
! 461: NEXT(n) = 0;
! 462: MKLIST(*rp,n0);
! 463: }
! 464:
! 465: void clctalg(p,vl)
! 466: P p;
! 467: VL *vl;
! 468: {
! 469: int n,i;
! 470: VL tvl;
! 471: VN vn,vn1;
! 472: P d;
! 473: DCP dc;
! 474:
! 475: for ( n = 0, tvl = ALG; tvl; tvl = NEXT(tvl), n++ );
! 476: vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 477: for ( i = n-1, tvl = ALG; tvl; tvl = NEXT(tvl), i-- ) {
! 478: vn[i].v = tvl->v;
! 479: vn[i].n = 0;
! 480: }
! 481: markv(vn,n,p);
! 482: for ( i = n-1; i >= 0; i-- ) {
! 483: if ( !vn[i].n )
! 484: continue;
! 485: d = (P)vn[i].v->attr;
! 486: for ( dc = DC(d); dc; dc = NEXT(dc) )
! 487: markv(vn,i,COEF(dc));
! 488: }
! 489: vn1 = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 490: for ( i = 0; i < n; i++ ) {
! 491: vn1[i].v = vn[n-1-i].v; vn1[i].n = vn[n-1-i].n;
! 492: }
! 493: vntovl(vn1,n,vl);
! 494: }
! 495:
! 496: void Palg(arg,rp)
! 497: NODE arg;
! 498: Alg *rp;
! 499: {
! 500: Q a;
! 501: VL vl;
! 502: P x;
! 503: int n;
! 504:
! 505: a = (Q)ARG0(arg);
! 506: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
! 507: *rp = 0;
! 508: else {
! 509: n = ACNT-QTOS(a)-1;
! 510: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
! 511: if ( vl ) {
! 512: MKV(vl->v,x); MKAlg(x,*rp);
! 513: } else
! 514: *rp = 0;
! 515: }
! 516: }
! 517:
! 518: void Palgv(arg,rp)
! 519: NODE arg;
! 520: Obj *rp;
! 521: {
! 522: Q a;
! 523: VL vl;
! 524: P x;
! 525: int n;
! 526: Alg b;
! 527:
! 528: a = (Q)ARG0(arg);
! 529: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
! 530: *rp = 0;
! 531: else {
! 532: n = ACNT-QTOS(a)-1;
! 533: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
! 534: if ( vl ) {
! 535: MKV(vl->v,x); MKAlg(x,b); algtorat((Num)b,rp);
! 536: } else
! 537: *rp = 0;
! 538: }
! 539: }
! 540:
! 541: void algptop(p,r)
! 542: P p,*r;
! 543: {
! 544: DCP dc,dcr,dcr0;
! 545:
! 546: if ( NUM(p) )
! 547: *r = (P)p;
! 548: else {
! 549: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 550: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 551: algptop(COEF(dc),&COEF(dcr));
! 552: }
! 553: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
! 554: }
! 555: }
! 556:
! 557: void algtorat(n,r)
! 558: Num n;
! 559: Obj *r;
! 560: {
! 561: Obj obj;
! 562: P nm,dn;
! 563:
! 564: if ( !n || NID(n) <= N_R )
! 565: *r = (Obj)n;
! 566: else {
! 567: obj = BDY((Alg)n);
! 568: if ( ID(obj) <= O_P )
! 569: algptop((P)obj,(P *)r);
! 570: else {
! 571: algptop(NM((R)obj),&nm); algptop(DN((R)obj),&dn);
! 572: divr(CO,(Obj)nm,(Obj)dn,r);
! 573: }
! 574: }
! 575: }
! 576:
! 577: void rattoalg(obj,n)
! 578: Obj obj;
! 579: Alg *n;
! 580: {
! 581: P nm,dn;
! 582: Obj t;
! 583:
! 584: if ( !obj || ID(obj) == O_N )
! 585: *n = (Alg)obj;
! 586: else if ( ID(obj) == O_P ) {
! 587: ptoalgp((P)obj,(P *)&t); MKAlg(t,*n);
! 588: } else {
! 589: ptoalgp(NM((R)obj),&nm); ptoalgp(DN((R)obj),&dn);
! 590: divr(ALG,(Obj)nm,(Obj)dn,&t); MKAlg(t,*n);
! 591: }
! 592: }
! 593:
! 594: void ptoalgp(p,r)
! 595: P p,*r;
! 596: {
! 597: DCP dc,dcr,dcr0;
! 598:
! 599: if ( NUM(p) )
! 600: *r = (P)p;
! 601: else {
! 602: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 603: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 604: ptoalgp(COEF(dc),&COEF(dcr));
! 605: }
! 606: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
! 607: }
! 608: }
! 609:
! 610: void Pinvalg_chrem(NODE arg,LIST *r)
! 611: {
! 612: NODE n;
! 613:
! 614: inva_chrem((P)ARG0(arg),(P)ARG1(arg),&n);
! 615: MKLIST(*r,n);
! 616: }
! 617:
! 618: void invalg_le(Alg a,LIST *r);
! 619:
! 620: void Pinvalg_le(NODE arg,LIST *r)
! 621: {
! 622: invalg_le((Alg)ARG0(arg),r);
! 623: }
! 624:
! 625: typedef struct oMono_nf {
! 626: DP mono;
! 627: DP nf;
! 628: Z dn;
! 629: } *Mono_nf;
! 630:
! 631: void invalg_le(Alg a,LIST *r)
! 632: {
! 633: Alg inv;
! 634: MAT mobj,sol;
! 635: int *rinfo,*cinfo;
! 636: P p,dn,ap;
! 637: VL vl,tvl;
! 638: Q c1,c2,c3,cont,c,mul;
! 639: Z two,iq,dn0,dn1,dnsol;
! 640: int i,j,n,len,k;
! 641: MP mp,mp0;
! 642: DP dp,nm,nm1,m,d,u,u1;
! 643: NODE b,b1,hlist,mblist,t,s,rev0,rev,hist;
! 644: DP *ps;
! 645: struct order_spec *spec;
! 646: Mono_nf h,h1;
! 647: Z nq,nr,nl,ng;
! 648: Q **mat,**solmat;
! 649: Q *w;
! 650: int *wi;
! 651:
! 652: ap = (P)BDY(a);
! 653: asir_assert(ap,O_P,"invalg_le");
! 654:
! 655: /* collecting algebraic numbers */
! 656: clctalg(ap,&vl);
! 657:
! 658: /* setup */
! 659: ptozp(ap,1,&c,&p);
! 660: STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec);
! 661: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 662: ps = (DP *)ALLOCA(n*sizeof(DP));
! 663:
! 664: /* conversion to DP */
! 665: for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) {
! 666: ptod(ALG,vl,tvl->v->attr,&ps[i]);
! 667: }
! 668: ptod(ALG,vl,p,&dp);
! 669: /* index list */
! 670: for ( b = 0, i = 0; i < n; i++ ) {
! 671: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
! 672: }
! 673: /* simplification */
! 674: dp_true_nf(b,dp,ps,1,&nm,(P *)&dn);
! 675:
! 676: /* construction of NF table */
! 677:
! 678: /* stdmono: <<0,...,0>> < ... < max */
! 679: for ( hlist = 0, i = 0; i < n; i++ ) {
! 680: MKNODE(b1,(pointer)ps[i],hlist); hlist = b1;
! 681: }
! 682: dp_mbase(hlist,&rev0);
! 683: for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) {
! 684: MKNODE(b1,BDY(rev),mblist); mblist = b1;
! 685: }
! 686: dn0 = ONE;
! 687: for ( hist = 0, t = mblist; t; t = NEXT(t) ) {
! 688: /* searching a predecessor */
! 689: for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) {
! 690: h = (Mono_nf)BDY(s);
! 691: if ( dp_redble(m,h->mono) )
! 692: break;
! 693: }
! 694: h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf));
! 695: if ( s ) {
! 696: dp_subd(m,h->mono,&d);
! 697: muld(CO,d,h->nf,&u);
! 698: dp_true_nf(b,u,ps,1,&nm1,(P *)&dn1);
! 699: mulz(h->dn,dn1,&h1->dn);
! 700: } else {
! 701: muld(CO,m,nm,&u);
! 702: dp_true_nf(b,u,ps,1,&nm1,(P *)&dn1);
! 703: h1->dn = dn1;
! 704: }
! 705: h1->mono = m;
! 706: h1->nf = nm1;
! 707: MKNODE(b1,(pointer)h1,hist); hist = b1;
! 708:
! 709: /* dn0 = LCM(dn0,h1->dn) */
! 710: gcdz(dn0,h1->dn,&ng); divz(dn0,ng,&nq);
! 711: mulz(nq,h1->dn,&nl); absz(nl,&dn0);
! 712: }
! 713: /* create a matrix */
! 714: len = length(mblist);
! 715: MKMAT(mobj,len,len+1);
! 716: mat = (Q **)BDY(mobj);
! 717: mat[len-1][len] = (Q)dn0;
! 718: for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) {
! 719: h = (Mono_nf)BDY(t);
! 720: nm1 = h->nf;
! 721: divq((Q)dn0,(Q)h->dn,&mul);
! 722: for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) )
! 723: if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) {
! 724: mulq(mul,(Q)mp->c,&mat[i][j]);
! 725: mp = NEXT(mp);
! 726: }
! 727: }
! 728: generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
! 729: solmat = (Q **)BDY(sol);
! 730: for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) )
! 731: if ( solmat[i][0] ) {
! 732: NEXTMP(mp0,mp);
! 733: mp->c = (Obj)solmat[i][0];
! 734: mp->dl = BDY((DP)BDY(t))->dl;
! 735: }
! 736: NEXT(mp) = 0; MKDP(n,mp0,u);
! 737: dp_ptozp(u,&u1);
! 738: divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont);
! 739: dtop(ALG,vl,u1,(Obj *)&ap);
! 740: MKAlg(ap,inv);
! 741: mulq((Q)dnsol,(Q)dn,&c1);
! 742: mulq(c1,c,&c2);
! 743: divq(c2,cont,&c3);
! 744: b = mknode(2,inv,c3);
! 745: MKLIST(*r,b);
! 746: }
! 747:
! 748: void get_algtree(Obj f,VL *r)
! 749: {
! 750: VL vl1,vl2,vl3;
! 751: Obj t;
! 752: DCP dc;
! 753: NODE b;
! 754: pointer *a;
! 755: pointer **m;
! 756: int len,row,col,i,j,l;
! 757:
! 758: if ( !f ) *r = 0;
! 759: else
! 760: switch ( OID(f) ) {
! 761: case O_N:
! 762: if ( NID((Num)f) != N_A ) *r = 0;
! 763: else {
! 764: t = BDY((Alg)f);
! 765: switch ( OID(t) ) {
! 766: case O_P:
! 767: clctalg((P)t,r); break;
! 768: case O_R:
! 769: clctalg(NM((R)t),&vl1);
! 770: clctalg(DN((R)t),&vl2);
! 771: mergev(ALG,vl1,vl2,r); break;
! 772: default:
! 773: *r = 0; break;
! 774: }
! 775: }
! 776: break;
! 777: case O_P:
! 778: vl1 = 0;
! 779: for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
! 780: get_algtree((Obj)COEF(dc),&vl2);
! 781: mergev(ALG,vl1,vl2,&vl3);
! 782: vl1 = vl3;
! 783: }
! 784: *r = vl1;
! 785: break;
! 786: case O_R:
! 787: get_algtree((Obj)NM((R)f),&vl1);
! 788: get_algtree((Obj)DN((R)f),&vl2);
! 789: mergev(ALG,vl1,vl2,r);
! 790: break;
! 791: case O_LIST:
! 792: vl1 = 0;
! 793: for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
! 794: get_algtree((Obj)BDY(b),&vl2);
! 795: mergev(ALG,vl1,vl2,&vl3);
! 796: vl1 = vl3;
! 797: }
! 798: *r = vl1;
! 799: break;
! 800: case O_VECT:
! 801: vl1 = 0;
! 802: l = ((VECT)f)->len;
! 803: a = BDY((VECT)f);
! 804: for ( i = 0; i < l; i++ ) {
! 805: get_algtree((Obj)a[i],&vl2);
! 806: mergev(ALG,vl1,vl2,&vl3);
! 807: vl1 = vl3;
! 808: }
! 809: *r = vl1;
! 810: break;
! 811: case O_MAT:
! 812: vl1 = 0;
! 813: row = ((MAT)f)->row; col = ((MAT)f)->col;
! 814: m = BDY((MAT)f);
! 815: for ( i = 0; i < row; i++ )
! 816: for ( j = 0; j < col; j++ ) {
! 817: get_algtree((Obj)m[i][j],&vl2);
! 818: mergev(ALG,vl1,vl2,&vl3);
! 819: vl1 = vl3;
! 820: }
! 821: *r = vl1;
! 822: break;
! 823: default:
! 824: *r = 0;
! 825: break;
! 826: }
! 827: }
! 828:
! 829: void algobjtorat(Obj f,Obj *r)
! 830: {
! 831: Obj t;
! 832: DCP dc,dcr,dcr0;
! 833: P p,nm,dn;
! 834: R rat;
! 835: NODE b,s,s0;
! 836: VECT v;
! 837: MAT mat;
! 838: LIST list;
! 839: pointer *a;
! 840: pointer **m;
! 841: int len,row,col,i,j,l;
! 842:
! 843: if ( !f ) *r = 0;
! 844: else
! 845: switch ( OID(f) ) {
! 846: case O_N:
! 847: algtorat((Num)f,r);
! 848: break;
! 849: case O_P:
! 850: dcr0 = 0;
! 851: for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
! 852: NEXTDC(dcr0,dcr);
! 853: algobjtorat((Obj)COEF(dc),&t);
! 854: COEF(dcr) = (P)t;
! 855: DEG(dcr) = DEG(dc);
! 856: }
! 857: NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p;
! 858: break;
! 859: case O_R:
! 860: algobjtorat((Obj)NM((R)f),&t); nm = (P)t;
! 861: algobjtorat((Obj)DN((R)f),&t); dn = (P)t;
! 862: MKRAT(nm,dn,0,rat); *r = (Obj)rat;
! 863: break;
! 864: case O_LIST:
! 865: s0 = 0;
! 866: for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
! 867: NEXTNODE(s0,s);
! 868: algobjtorat((Obj)BDY(b),&t);
! 869: BDY(s) = (pointer)t;
! 870: }
! 871: NEXT(s) = 0;
! 872: MKLIST(list,s0);
! 873: *r = (Obj)list;
! 874: break;
! 875: case O_VECT:
! 876: l = ((VECT)f)->len;
! 877: a = BDY((VECT)f);
! 878: MKVECT(v,l);
! 879: for ( i = 0; i < l; i++ ) {
! 880: algobjtorat((Obj)a[i],&t);
! 881: BDY(v)[i] = (pointer)t;
! 882: }
! 883: *r = (Obj)v;
! 884: break;
! 885: case O_MAT:
! 886: row = ((MAT)f)->row; col = ((MAT)f)->col;
! 887: m = BDY((MAT)f);
! 888: MKMAT(mat,row,col);
! 889: for ( i = 0; i < row; i++ )
! 890: for ( j = 0; j < col; j++ ) {
! 891: algobjtorat((Obj)m[i][j],&t);
! 892: BDY(mat)[i][j] = (pointer)t;
! 893: }
! 894: *r = (Obj)mat;
! 895: break;
! 896: default:
! 897: *r = f;
! 898: break;
! 899: }
! 900: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to algnum.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin