Annotation of OpenXM_contrib2/asir2000/builtin/algnum.c, Revision 1.16
1.2 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
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 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: *
1.16 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.15 2017/08/31 02:36:20 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();
54: void Palg(), Palgv(), Pgetalgtree();
1.6 noro 55: void Pinvalg_le();
1.7 noro 56: void Pset_field(),Palgtodalg(),Pdalgtoalg();
1.10 noro 57: void Pinv_or_split_dalg();
1.11 noro 58: void Pdalgtoup();
59: void Pget_field_defpoly();
60: void Pget_field_generator();
1.1 noro 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 *);
1.4 noro 68: void clctalg(P,VL *);
1.8 noro 69: void get_algtree(Obj f,VL *r);
1.11 noro 70: void Pinvalg_chrem();
71: void Pdalgtodp();
72: void Pdptodalg();
1.1 noro 73:
74: struct ftab alg_tab[] = {
1.16 ! noro 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},
1.1 noro 96: };
97:
98: static int UCN,ACNT;
1.7 noro 99:
100: void Pset_field(NODE arg,Q *rp)
101: {
1.16 ! noro 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;
1.7 noro 121: }
122:
123: void Palgtodalg(NODE arg,DAlg *rp)
124: {
1.16 ! noro 125: algtodalg((Alg)ARG0(arg),rp);
1.7 noro 126: }
127:
128: void Pdalgtoalg(NODE arg,Alg *rp)
129: {
1.16 ! noro 130: dalgtoalg((DAlg)ARG0(arg),rp);
1.10 noro 131: }
132:
1.11 noro 133: void Pdalgtodp(NODE arg,LIST *r)
134: {
1.16 ! noro 135: NODE b;
! 136: DP nm;
! 137: Q 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);
1.11 noro 145: }
146:
147: void Pdptodalg(NODE arg,DAlg *r)
148: {
1.16 ! noro 149: DP d,nm,nm1;
! 150: MP m;
! 151: Q c,a;
! 152: DAlg t;
! 153:
! 154: d = (DP)ARG0(arg);
! 155: if ( !d ) *r = 0;
! 156: else {
! 157: for ( m = BDY(d); m; m = NEXT(m) )
! 158: if ( !INT((Q)m->c) ) break;
! 159: if ( !m ) {
! 160: MKDAlg(d,(Q)ONE,t);
! 161: } else {
! 162: dp_ptozp(d,&nm);
! 163: divq((Q)BDY(d)->c,(Q)BDY(nm)->c,&c);
! 164: NTOQ(NM(c),SGN(c),a);
! 165: muldc(CO,nm,(Obj)a,&nm1);
! 166: NTOQ(DN(c),1,a);
! 167: MKDAlg(nm1,a,t);
! 168: }
! 169: simpdalg(t,r);
! 170: }
1.11 noro 171: }
172:
173: void Pdalgtoup(NODE arg,LIST *r)
174: {
1.16 ! noro 175: NODE b;
! 176: int pos;
! 177: P up;
! 178: DP nm;
! 179: Q dn,q;
! 180:
! 181: pos = dalgtoup((DAlg)ARG0(arg),&up,&dn);
! 182: STOQ(pos,q);
! 183: b = mknode(3,up,dn,q);
! 184: MKLIST(*r,b);
1.11 noro 185: }
186:
1.10 noro 187: NODE inv_or_split_dalg(DAlg,DAlg *);
1.16 ! noro 188: NumberField get_numberfield();
1.10 noro 189:
1.11 noro 190: void Pget_field_defpoly(NODE arg,DAlg *r)
191: {
1.16 ! noro 192: NumberField nf;
! 193: DP d;
1.11 noro 194:
1.16 ! noro 195: nf = get_numberfield();
! 196: d = nf->ps[QTOS((Q)ARG0(arg))];
! 197: MKDAlg(d,ONE,*r);
1.11 noro 198: }
199:
200: void Pget_field_generator(NODE arg,DAlg *r)
201: {
1.16 ! noro 202: int index,n,i;
! 203: DL dl;
! 204: MP m;
! 205: DP d;
! 206:
! 207: index = QTOS((Q)ARG0(arg));
! 208: n = get_numberfield()->n;
! 209: NEWDL(dl,n);
! 210: for ( i = 0; i < n; i++ ) dl->d[i] = 0;
! 211: dl->d[index] = 1; dl->td = 1;
! 212: NEWMP(m); m->dl = dl; m->c = (Obj)ONE; NEXT(m) = 0;
! 213: MKDP(n,m,d);
! 214: MKDAlg(d,ONE,*r);
1.11 noro 215: }
216:
217:
1.10 noro 218: void Pinv_or_split_dalg(NODE arg,Obj *rp)
219: {
1.16 ! noro 220: NODE gen,t,nd0,nd;
! 221: LIST list;
! 222: int l,i,j,n;
! 223: DP *ps,*ps1,*psw;
! 224: NumberField nf;
! 225: DAlg inv;
! 226: extern struct order_spec *dp_current_spec;
! 227: struct order_spec *current_spec;
! 228:
! 229: gen = inv_or_split_dalg((DAlg)ARG0(arg),&inv);
! 230: if ( !gen )
! 231: *rp = (Obj)inv;
! 232: else {
! 233: nf = get_numberfield();
! 234: current_spec = dp_current_spec; initd(nf->spec);
! 235: l = length(gen);
! 236: n = nf->n;
! 237: ps = nf->ps;
! 238: psw = (DP *)ALLOCA((n+l)*sizeof(DP));
! 239: for ( i = j = 0; i < n; i++ ) {
! 240: for ( t = gen; t; t = NEXT(t) )
! 241: if ( dp_redble(ps[i],(DP)BDY(t)) ) break;
! 242: if ( !t )
! 243: psw[j++] = ps[i];
! 244: }
! 245: nd0 = 0;
! 246: /* gen[0] < gen[1] < ... */
! 247: /* psw[0] > psw[1] > ... */
! 248: for ( i = j-1, t = gen; i >= 0 && t; ) {
! 249: NEXTNODE(nd0,nd);
! 250: if ( compd(CO,psw[i],(DP)BDY(t)) > 0 ) {
! 251: BDY(nd) = BDY(t); t = NEXT(t);
! 252: } else
! 253: BDY(nd) = (pointer)psw[i--];
! 254: }
! 255: for ( ; i >= 0; i-- ) {
! 256: NEXTNODE(nd0,nd); BDY(nd) = (pointer)psw[i];
! 257: }
! 258: for ( ; t; t = NEXT(t) ) {
! 259: NEXTNODE(nd0,nd); BDY(nd) = BDY(t);
! 260: }
! 261: NEXT(nd) = 0;
! 262: MKLIST(list,nd0);
! 263: initd(current_spec);
! 264: *rp = (Obj)list;
! 265: }
1.7 noro 266: }
1.1 noro 267:
268: void Pnewalg(arg,rp)
269: NODE arg;
270: Alg *rp;
271: {
1.16 ! noro 272: P p;
! 273: VL vl;
! 274: P c;
! 275:
! 276: p = (P)ARG0(arg);
! 277: if ( !p || OID(p) != O_P )
! 278: error("newalg : invalid argument");
! 279: clctv(CO,p,&vl);
! 280: if ( NEXT(vl) )
! 281: error("newalg : invalid argument");
! 282: c = COEF(DC(p));
! 283: if ( !NUM(c) || !RATN(c) )
! 284: error("newalg : invalid argument");
! 285: mkalg(p,rp);
1.1 noro 286: }
287:
288: void mkalg(p,r)
289: P p;
290: Alg *r;
291: {
1.16 ! noro 292: VL vl,mvl,nvl;
! 293: V a,tv;
! 294: char buf[BUFSIZ];
! 295: char *name;
! 296: P x,t,s;
! 297: Num c;
! 298: DCP dc,dcr,dcr0;
! 299:
! 300: for ( vl = ALG; vl; vl = NEXT(vl) )
! 301: if ( !cmpalgp(p,(P)vl->v->attr) ) {
! 302: a = vl->v; break;
! 303: }
! 304: if ( !vl ) {
! 305: NEWVL(vl); NEXT(vl) = ALG; ALG = vl;
! 306: NEWV(a); vl->v = a;
! 307: sprintf(buf,"#%d",ACNT++);
! 308: name = (char *)MALLOC(strlen(buf)+1);
! 309: strcpy(name,buf); NAME(a) = name;
! 310:
! 311: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 312: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); c = (Num)COEF(dc);
! 313: if ( NID(c) != N_A )
! 314: COEF(dcr) = (P)c;
! 315: else
! 316: COEF(dcr) = (P)BDY(((Alg)c));
! 317: }
! 318: NEXT(dcr) = 0; MKP(a,dcr0,t); a->attr = (pointer)t;
! 319:
! 320: sprintf(buf,"t%s",name); makevar(buf,&s);
! 321:
! 322: if ( NEXT(ALG) ) {
! 323: tv = (V)NEXT(ALG)->v->priv;
! 324: for ( vl = CO; NEXT(NEXT(vl)); vl = NEXT(vl) );
! 325: nvl = NEXT(vl); NEXT(vl) = 0;
! 326: for ( vl = CO; NEXT(vl) && (NEXT(vl)->v != tv); vl = NEXT(vl) );
! 327: mvl = NEXT(vl); NEXT(vl) = nvl; NEXT(nvl) = mvl;
! 328: }
! 329:
! 330: a->priv = (pointer)VR(s); VR(s)->priv = (pointer)a;
! 331: }
! 332: MKV(a,x); MKAlg(x,*r);
1.1 noro 333: }
334:
335: int cmpalgp(p,defp)
336: P p,defp;
337: {
1.16 ! noro 338: DCP dc,dcd;
! 339: P t;
1.1 noro 340:
1.16 ! noro 341: for ( dc = DC(p), dcd = DC(defp); dc && dcd;
! 342: dc = NEXT(dc), dcd = NEXT(dcd) ) {
! 343: if ( cmpq(DEG(dc),DEG(dcd)) )
! 344: break;
! 345: t = NID((Num)COEF(dc)) == N_A ? (P)BDY((Alg)COEF(dc)) : COEF(dc);
! 346: if ( compp(ALG,t,COEF(dcd)) )
! 347: break;
! 348: }
! 349: if ( dc || dcd )
! 350: return 1;
! 351: else
! 352: return 0;
1.1 noro 353: }
354:
355: void Pdefpoly(arg,rp)
356: NODE arg;
357: P *rp;
358: {
1.16 ! noro 359: asir_assert(ARG0(arg),O_N,"defpoly");
! 360: algptop((P)VR((P)BDY((Alg)ARG0(arg)))->attr,rp);
1.1 noro 361: }
362:
363: void Pmainalg(arg,r)
364: NODE arg;
365: Alg *r;
366: {
1.16 ! noro 367: Num c;
! 368: V v;
! 369: P b;
! 370:
! 371: c = (Num)(ARG0(arg));
! 372: if ( NID(c) <= N_R )
! 373: *r = 0;
! 374: else {
! 375: v = VR((P)BDY((Alg)c)); MKV(v,b); MKAlg(b,*r);
! 376: }
1.1 noro 377: }
378:
379: void Palgtorat(arg,rp)
380: NODE arg;
381: Obj *rp;
382: {
1.16 ! noro 383: asir_assert(ARG0(arg),O_N,"algtorat");
! 384: algtorat((Num)ARG0(arg),rp);
1.1 noro 385: }
386:
387: void Prattoalg(arg,rp)
388: NODE arg;
389: Alg *rp;
390: {
1.16 ! noro 391: asir_assert(ARG0(arg),O_R,"rattoalg");
! 392: rattoalg((Obj)ARG0(arg),rp);
1.1 noro 393: }
394:
395: void Pgetalg(arg,rp)
396: NODE arg;
397: LIST *rp;
398: {
1.16 ! noro 399: Obj t;
! 400: P p;
! 401: VL vl;
! 402: Num a;
! 403: Alg b;
! 404: NODE n0,n;
! 405:
! 406: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
! 407: vl = 0;
! 408: else {
! 409: t = BDY((Alg)a);
! 410: switch ( OID(t) ) {
! 411: case O_P: case O_R:
! 412: clctvr(ALG,t,&vl); break;
! 413: default:
! 414: vl = 0; break;
! 415: }
! 416: }
! 417: for ( n0 = 0; vl; vl = NEXT(vl) ) {
! 418: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
! 419: }
! 420: if ( n0 )
! 421: NEXT(n) = 0;
! 422: MKLIST(*rp,n0);
1.1 noro 423: }
424:
425: void Pgetalgtree(arg,rp)
426: NODE arg;
427: LIST *rp;
428: {
1.16 ! noro 429: Obj t;
! 430: P p;
! 431: VL vl,vl1,vl2;
! 432: Num a;
! 433: Alg b;
! 434: NODE n0,n;
1.1 noro 435:
1.8 noro 436: #if 0
1.16 ! noro 437: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
! 438: vl = 0;
! 439: else {
! 440: t = BDY((Alg)a);
! 441: switch ( OID(t) ) {
! 442: case O_P:
! 443: clctalg((P)t,&vl); break;
! 444: case O_R:
! 445: clctalg(NM((R)t),&vl1);
! 446: clctalg(DN((R)t),&vl2);
! 447: mergev(ALG,vl1,vl2,&vl); break;
! 448: default:
! 449: vl = 0; break;
! 450: }
! 451: }
1.8 noro 452: #else
1.16 ! noro 453: get_algtree((Obj)ARG0(arg),&vl);
1.8 noro 454: #endif
1.16 ! noro 455: for ( n0 = 0; vl; vl = NEXT(vl) ) {
! 456: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
! 457: }
! 458: if ( n0 )
! 459: NEXT(n) = 0;
! 460: MKLIST(*rp,n0);
1.1 noro 461: }
462:
463: void clctalg(p,vl)
464: P p;
465: VL *vl;
466: {
1.16 ! noro 467: int n,i;
! 468: VL tvl;
! 469: VN vn,vn1;
! 470: P d;
! 471: DCP dc;
! 472:
! 473: for ( n = 0, tvl = ALG; tvl; tvl = NEXT(tvl), n++ );
! 474: vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 475: for ( i = n-1, tvl = ALG; tvl; tvl = NEXT(tvl), i-- ) {
! 476: vn[i].v = tvl->v;
! 477: vn[i].n = 0;
! 478: }
! 479: markv(vn,n,p);
! 480: for ( i = n-1; i >= 0; i-- ) {
! 481: if ( !vn[i].n )
! 482: continue;
! 483: d = (P)vn[i].v->attr;
! 484: for ( dc = DC(d); dc; dc = NEXT(dc) )
! 485: markv(vn,i,COEF(dc));
! 486: }
! 487: vn1 = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 488: for ( i = 0; i < n; i++ ) {
! 489: vn1[i].v = vn[n-1-i].v; vn1[i].n = vn[n-1-i].n;
! 490: }
! 491: vntovl(vn1,n,vl);
1.1 noro 492: }
493:
494: void Palg(arg,rp)
495: NODE arg;
496: Alg *rp;
497: {
1.16 ! noro 498: Q a;
! 499: VL vl;
! 500: P x;
! 501: int n;
! 502:
! 503: a = (Q)ARG0(arg);
! 504: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
! 505: *rp = 0;
! 506: else {
! 507: n = ACNT-QTOS(a)-1;
! 508: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
! 509: if ( vl ) {
! 510: MKV(vl->v,x); MKAlg(x,*rp);
! 511: } else
! 512: *rp = 0;
! 513: }
1.1 noro 514: }
515:
516: void Palgv(arg,rp)
517: NODE arg;
518: Obj *rp;
519: {
1.16 ! noro 520: Q a;
! 521: VL vl;
! 522: P x;
! 523: int n;
! 524: Alg b;
! 525:
! 526: a = (Q)ARG0(arg);
! 527: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
! 528: *rp = 0;
! 529: else {
! 530: n = ACNT-QTOS(a)-1;
! 531: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
! 532: if ( vl ) {
! 533: MKV(vl->v,x); MKAlg(x,b); algtorat((Num)b,rp);
! 534: } else
! 535: *rp = 0;
! 536: }
1.1 noro 537: }
538:
539: void algptop(p,r)
540: P p,*r;
541: {
1.16 ! noro 542: DCP dc,dcr,dcr0;
1.1 noro 543:
1.16 ! noro 544: if ( NUM(p) )
! 545: *r = (P)p;
! 546: else {
! 547: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 548: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 549: algptop(COEF(dc),&COEF(dcr));
! 550: }
! 551: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
! 552: }
1.1 noro 553: }
554:
555: void algtorat(n,r)
556: Num n;
557: Obj *r;
558: {
1.16 ! noro 559: Obj obj;
! 560: P nm,dn;
1.1 noro 561:
1.16 ! noro 562: if ( !n || NID(n) <= N_R )
! 563: *r = (Obj)n;
! 564: else {
! 565: obj = BDY((Alg)n);
! 566: if ( ID(obj) <= O_P )
! 567: algptop((P)obj,(P *)r);
! 568: else {
! 569: algptop(NM((R)obj),&nm); algptop(DN((R)obj),&dn);
! 570: divr(CO,(Obj)nm,(Obj)dn,r);
! 571: }
! 572: }
1.1 noro 573: }
574:
575: void rattoalg(obj,n)
576: Obj obj;
577: Alg *n;
578: {
1.16 ! noro 579: P nm,dn;
! 580: Obj t;
1.1 noro 581:
1.16 ! noro 582: if ( !obj || ID(obj) == O_N )
! 583: *n = (Alg)obj;
! 584: else if ( ID(obj) == O_P ) {
! 585: ptoalgp((P)obj,(P *)&t); MKAlg(t,*n);
! 586: } else {
! 587: ptoalgp(NM((R)obj),&nm); ptoalgp(DN((R)obj),&dn);
! 588: divr(ALG,(Obj)nm,(Obj)dn,&t); MKAlg(t,*n);
! 589: }
1.1 noro 590: }
591:
592: void ptoalgp(p,r)
593: P p,*r;
594: {
1.16 ! noro 595: DCP dc,dcr,dcr0;
1.1 noro 596:
1.16 ! noro 597: if ( NUM(p) )
! 598: *r = (P)p;
! 599: else {
! 600: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 601: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 602: ptoalgp(COEF(dc),&COEF(dcr));
! 603: }
! 604: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
! 605: }
1.11 noro 606: }
607:
608: void Pinvalg_chrem(NODE arg,LIST *r)
609: {
1.16 ! noro 610: NODE n;
1.11 noro 611:
1.16 ! noro 612: inva_chrem((P)ARG0(arg),(P)ARG1(arg),&n);
! 613: MKLIST(*r,n);
1.6 noro 614: }
615:
616: void invalg_le(Alg a,LIST *r);
617:
618: void Pinvalg_le(NODE arg,LIST *r)
619: {
1.16 ! noro 620: invalg_le((Alg)ARG0(arg),r);
1.6 noro 621: }
622:
623: typedef struct oMono_nf {
1.16 ! noro 624: DP mono;
! 625: DP nf;
! 626: Q dn;
1.6 noro 627: } *Mono_nf;
628:
629: void invalg_le(Alg a,LIST *r)
630: {
1.16 ! noro 631: Alg inv;
! 632: MAT mobj,sol;
! 633: int *rinfo,*cinfo;
! 634: P p,dn,dn1,ap;
! 635: VL vl,tvl;
! 636: Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol;
! 637: int i,j,n,len,k;
! 638: MP mp,mp0;
! 639: DP dp,nm,nm1,m,d,u,u1;
! 640: NODE b,b1,hlist,mblist,t,s,rev0,rev,hist;
! 641: DP *ps;
! 642: struct order_spec *spec;
! 643: Mono_nf h,h1;
! 644: N nq,nr,nl,ng;
! 645: Q **mat,**solmat;
! 646: Q *w;
! 647: int *wi;
! 648:
! 649: ap = (P)BDY(a);
! 650: asir_assert(ap,O_P,"invalg_le");
! 651:
! 652: /* collecting algebraic numbers */
! 653: clctalg(ap,&vl);
! 654:
! 655: /* setup */
! 656: ptozp(ap,1,&c,&p);
! 657: STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec);
! 658: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 659: ps = (DP *)ALLOCA(n*sizeof(DP));
! 660:
! 661: /* conversion to DP */
! 662: for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) {
! 663: ptod(ALG,vl,tvl->v->attr,&ps[i]);
! 664: }
! 665: ptod(ALG,vl,p,&dp);
! 666: /* index list */
! 667: for ( b = 0, i = 0; i < n; i++ ) {
! 668: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
! 669: }
! 670: /* simplification */
! 671: dp_true_nf(b,dp,ps,1,&nm,&dn);
! 672:
! 673: /* construction of NF table */
! 674:
! 675: /* stdmono: <<0,...,0>> < ... < max */
! 676: for ( hlist = 0, i = 0; i < n; i++ ) {
! 677: MKNODE(b1,(pointer)ps[i],hlist); hlist = b1;
! 678: }
! 679: dp_mbase(hlist,&rev0);
! 680: for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) {
! 681: MKNODE(b1,BDY(rev),mblist); mblist = b1;
! 682: }
! 683: dn0 = ONE;
! 684: for ( hist = 0, t = mblist; t; t = NEXT(t) ) {
! 685: /* searching a predecessor */
! 686: for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) {
! 687: h = (Mono_nf)BDY(s);
! 688: if ( dp_redble(m,h->mono) )
! 689: break;
! 690: }
! 691: h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf));
! 692: if ( s ) {
! 693: dp_subd(m,h->mono,&d);
! 694: muld(CO,d,h->nf,&u);
! 695: dp_true_nf(b,u,ps,1,&nm1,&dn1);
! 696: mulq(h->dn,(Q)dn1,&h1->dn);
! 697: } else {
! 698: muld(CO,m,nm,&u);
! 699: dp_true_nf(b,u,ps,1,&nm1,&dn1);
! 700: h1->dn = (Q)dn1;
! 701: }
! 702: h1->mono = m;
! 703: h1->nf = nm1;
! 704: MKNODE(b1,(pointer)h1,hist); hist = b1;
! 705:
! 706: /* dn0 = LCM(dn0,h1->dn) */
! 707: gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr);
! 708: muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0);
! 709: }
! 710: /* create a matrix */
! 711: len = length(mblist);
! 712: MKMAT(mobj,len,len+1);
! 713: mat = (Q **)BDY(mobj);
! 714: mat[len-1][len] = dn0;
! 715: for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) {
! 716: h = (Mono_nf)BDY(t);
! 717: nm1 = h->nf;
! 718: divq((Q)dn0,h->dn,&mul);
! 719: for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) )
! 720: if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) {
! 721: mulq(mul,(Q)mp->c,&mat[i][j]);
! 722: mp = NEXT(mp);
! 723: }
! 724: }
1.6 noro 725: #if 0
1.16 ! noro 726: w = (Q *)ALLOCA((len+1)*sizeof(Q));
! 727: wi = (int *)ALLOCA((len+1)*sizeof(int));
! 728: for ( i = 0; i < len; i++ ) {
! 729: for ( j = 0, k = 0; j <= len; j++ )
! 730: if ( mat[i][j] ) {
! 731: w[k] = mat[i][j];
! 732: wi[k] = j;
! 733: k++;
! 734: }
! 735: removecont_array(w,k);
! 736: for ( j = 0; j < k; j++ )
! 737: mat[i][wi[j]] = w[j];
! 738: }
1.6 noro 739: #endif
1.16 ! noro 740: generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
! 741: solmat = (Q **)BDY(sol);
! 742: for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) )
! 743: if ( solmat[i][0] ) {
! 744: NEXTMP(mp0,mp);
! 745: mp->c = (Obj)solmat[i][0];
! 746: mp->dl = BDY((DP)BDY(t))->dl;
! 747: }
! 748: NEXT(mp) = 0; MKDP(n,mp0,u);
! 749: dp_ptozp(u,&u1);
! 750: divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont);
! 751: dtop(ALG,vl,u1,(Obj *)&ap);
! 752: MKAlg(ap,inv);
! 753: mulq(dnsol,(Q)dn,&c1);
! 754: mulq(c1,c,&c2);
! 755: divq(c2,cont,&c3);
! 756: b = mknode(2,inv,c3);
! 757: MKLIST(*r,b);
1.8 noro 758: }
759:
760: void get_algtree(Obj f,VL *r)
761: {
1.16 ! noro 762: VL vl1,vl2,vl3;
! 763: Obj t;
! 764: DCP dc;
! 765: NODE b;
! 766: pointer *a;
! 767: pointer **m;
! 768: int len,row,col,i,j,l;
! 769:
! 770: if ( !f ) *r = 0;
! 771: else
! 772: switch ( OID(f) ) {
! 773: case O_N:
! 774: if ( NID((Num)f) != N_A ) *r = 0;
! 775: else {
! 776: t = BDY((Alg)f);
! 777: switch ( OID(t) ) {
! 778: case O_P:
! 779: clctalg((P)t,r); break;
! 780: case O_R:
! 781: clctalg(NM((R)t),&vl1);
! 782: clctalg(DN((R)t),&vl2);
! 783: mergev(ALG,vl1,vl2,r); break;
! 784: default:
! 785: *r = 0; break;
! 786: }
! 787: }
! 788: break;
! 789: case O_P:
! 790: vl1 = 0;
! 791: for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
! 792: get_algtree((Obj)COEF(dc),&vl2);
! 793: mergev(ALG,vl1,vl2,&vl3);
! 794: vl1 = vl3;
! 795: }
! 796: *r = vl1;
! 797: break;
! 798: case O_R:
! 799: get_algtree((Obj)NM((R)f),&vl1);
! 800: get_algtree((Obj)DN((R)f),&vl2);
! 801: mergev(ALG,vl1,vl2,r);
! 802: break;
! 803: case O_LIST:
! 804: vl1 = 0;
! 805: for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
! 806: get_algtree((Obj)BDY(b),&vl2);
! 807: mergev(ALG,vl1,vl2,&vl3);
! 808: vl1 = vl3;
! 809: }
! 810: *r = vl1;
! 811: break;
! 812: case O_VECT:
! 813: vl1 = 0;
! 814: l = ((VECT)f)->len;
! 815: a = BDY((VECT)f);
! 816: for ( i = 0; i < l; i++ ) {
! 817: get_algtree((Obj)a[i],&vl2);
! 818: mergev(ALG,vl1,vl2,&vl3);
! 819: vl1 = vl3;
! 820: }
! 821: *r = vl1;
! 822: break;
! 823: case O_MAT:
! 824: vl1 = 0;
! 825: row = ((MAT)f)->row; col = ((MAT)f)->col;
! 826: m = BDY((MAT)f);
! 827: for ( i = 0; i < row; i++ )
! 828: for ( j = 0; j < col; j++ ) {
! 829: get_algtree((Obj)m[i][j],&vl2);
! 830: mergev(ALG,vl1,vl2,&vl3);
! 831: vl1 = vl3;
! 832: }
! 833: *r = vl1;
! 834: break;
! 835: default:
! 836: *r = 0;
! 837: break;
! 838: }
1.9 noro 839: }
840:
841: void algobjtorat(Obj f,Obj *r)
842: {
1.16 ! noro 843: Obj t;
! 844: DCP dc,dcr,dcr0;
! 845: P p,nm,dn;
! 846: R rat;
! 847: NODE b,s,s0;
! 848: VECT v;
! 849: MAT mat;
! 850: LIST list;
! 851: pointer *a;
! 852: pointer **m;
! 853: int len,row,col,i,j,l;
! 854:
! 855: if ( !f ) *r = 0;
! 856: else
! 857: switch ( OID(f) ) {
! 858: case O_N:
! 859: algtorat((Num)f,r);
! 860: break;
! 861: case O_P:
! 862: dcr0 = 0;
! 863: for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
! 864: NEXTDC(dcr0,dcr);
! 865: algobjtorat((Obj)COEF(dc),&t);
! 866: COEF(dcr) = (P)t;
! 867: DEG(dcr) = DEG(dc);
! 868: }
! 869: NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p;
! 870: break;
! 871: case O_R:
! 872: algobjtorat((Obj)NM((R)f),&t); nm = (P)t;
! 873: algobjtorat((Obj)DN((R)f),&t); dn = (P)t;
! 874: MKRAT(nm,dn,0,rat); *r = (Obj)rat;
! 875: break;
! 876: case O_LIST:
! 877: s0 = 0;
! 878: for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
! 879: NEXTNODE(s0,s);
! 880: algobjtorat((Obj)BDY(b),&t);
! 881: BDY(s) = (pointer)t;
! 882: }
! 883: NEXT(s) = 0;
! 884: MKLIST(list,s0);
! 885: *r = (Obj)list;
! 886: break;
! 887: case O_VECT:
! 888: l = ((VECT)f)->len;
! 889: a = BDY((VECT)f);
! 890: MKVECT(v,l);
! 891: for ( i = 0; i < l; i++ ) {
! 892: algobjtorat((Obj)a[i],&t);
! 893: BDY(v)[i] = (pointer)t;
! 894: }
! 895: *r = (Obj)v;
! 896: break;
! 897: case O_MAT:
! 898: row = ((MAT)f)->row; col = ((MAT)f)->col;
! 899: m = BDY((MAT)f);
! 900: MKMAT(mat,row,col);
! 901: for ( i = 0; i < row; i++ )
! 902: for ( j = 0; j < col; j++ ) {
! 903: algobjtorat((Obj)m[i][j],&t);
! 904: BDY(mat)[i][j] = (pointer)t;
! 905: }
! 906: *r = (Obj)mat;
! 907: break;
! 908: default:
! 909: *r = f;
! 910: break;
! 911: }
1.1 noro 912: }
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 / asir2000 / builtin