![[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|
Return to algnum.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.5 2001/10/09 01:36:05 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.1 noro 56:
57: void mkalg(P,Alg *);
58: int cmpalgp(P,P);
59: void algptop(P,P *);
60: void algtorat(Num,Obj *);
61: void rattoalg(Obj,Alg *);
62: void ptoalgp(P,P *);
1.4 noro 63: void clctalg(P,VL *);
1.1 noro 64:
65: struct ftab alg_tab[] = {
1.6 ! noro 66: {"invalg_le",Pinvalg_le,1},
1.1 noro 67: {"defpoly",Pdefpoly,1},
68: {"newalg",Pnewalg,1},
69: {"mainalg",Pmainalg,1},
70: {"algtorat",Palgtorat,1},
71: {"rattoalg",Prattoalg,1},
72: {"getalg",Pgetalg,1},
73: {"getalgtree",Pgetalgtree,1},
74: {"alg",Palg,1},
75: {"algv",Palgv,1},
76: {0,0,0},
77: };
78:
79: static int UCN,ACNT;
80:
81: void Pnewalg(arg,rp)
82: NODE arg;
83: Alg *rp;
84: {
85: P p;
86: VL vl;
87: P c;
88:
89: p = (P)ARG0(arg);
90: if ( !p || OID(p) != O_P )
91: error("newalg : invalid argument");
92: clctv(CO,p,&vl);
93: if ( NEXT(vl) )
94: error("newalg : invalid argument");
95: c = COEF(DC(p));
96: if ( !NUM(c) || !RATN(c) )
97: error("newalg : invalid argument");
98: mkalg(p,rp);
99: }
100:
101: void mkalg(p,r)
102: P p;
103: Alg *r;
104: {
105: VL vl,mvl,nvl;
106: V a,tv;
107: char buf[BUFSIZ];
108: char *name;
109: P x,t,s;
110: Num c;
111: DCP dc,dcr,dcr0;
112:
113: for ( vl = ALG; vl; vl = NEXT(vl) )
114: if ( !cmpalgp(p,(P)vl->v->attr) ) {
115: a = vl->v; break;
116: }
117: if ( !vl ) {
118: NEWVL(vl); NEXT(vl) = ALG; ALG = vl;
119: NEWV(a); vl->v = a;
120: sprintf(buf,"#%d",ACNT++);
121: name = (char *)MALLOC(strlen(buf)+1);
122: strcpy(name,buf); NAME(a) = name;
123:
124: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
125: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); c = (Num)COEF(dc);
126: if ( NID(c) != N_A )
127: COEF(dcr) = (P)c;
128: else
129: COEF(dcr) = (P)BDY(((Alg)c));
130: }
131: NEXT(dcr) = 0; MKP(a,dcr0,t); a->attr = (pointer)t;
132:
133: sprintf(buf,"t%s",name); makevar(buf,&s);
134:
135: if ( NEXT(ALG) ) {
136: tv = (V)NEXT(ALG)->v->priv;
137: for ( vl = CO; NEXT(NEXT(vl)); vl = NEXT(vl) );
138: nvl = NEXT(vl); NEXT(vl) = 0;
139: for ( vl = CO; NEXT(vl) && (NEXT(vl)->v != tv); vl = NEXT(vl) );
140: mvl = NEXT(vl); NEXT(vl) = nvl; NEXT(nvl) = mvl;
141: }
142:
143: a->priv = (pointer)VR(s); VR(s)->priv = (pointer)a;
144: }
145: MKV(a,x); MKAlg(x,*r);
146: }
147:
148: int cmpalgp(p,defp)
149: P p,defp;
150: {
151: DCP dc,dcd;
152: P t;
153:
154: for ( dc = DC(p), dcd = DC(defp); dc && dcd;
155: dc = NEXT(dc), dcd = NEXT(dcd) ) {
156: if ( cmpq(DEG(dc),DEG(dcd)) )
157: break;
158: t = NID((Num)COEF(dc)) == N_A ? (P)BDY((Alg)COEF(dc)) : COEF(dc);
159: if ( compp(ALG,t,COEF(dcd)) )
160: break;
161: }
162: if ( dc || dcd )
163: return 1;
164: else
165: return 0;
166: }
167:
168: void Pdefpoly(arg,rp)
169: NODE arg;
170: P *rp;
171: {
172: asir_assert(ARG0(arg),O_N,"defpoly");
173: algptop((P)VR((P)BDY((Alg)ARG0(arg)))->attr,rp);
174: }
175:
176: void Pmainalg(arg,r)
177: NODE arg;
178: Alg *r;
179: {
180: Num c;
181: V v;
182: P b;
183:
184: c = (Num)(ARG0(arg));
185: if ( NID(c) <= N_R )
186: *r = 0;
187: else {
188: v = VR((P)BDY((Alg)c)); MKV(v,b); MKAlg(b,*r);
189: }
190: }
191:
192: void Palgtorat(arg,rp)
193: NODE arg;
194: Obj *rp;
195: {
196: asir_assert(ARG0(arg),O_N,"algtorat");
197: algtorat((Num)ARG0(arg),rp);
198: }
199:
200: void Prattoalg(arg,rp)
201: NODE arg;
202: Alg *rp;
203: {
204: asir_assert(ARG0(arg),O_R,"rattoalg");
205: rattoalg((Obj)ARG0(arg),rp);
206: }
207:
208: void Pgetalg(arg,rp)
209: NODE arg;
210: LIST *rp;
211: {
212: Obj t;
213: P p;
214: VL vl;
215: Num a;
216: Alg b;
217: NODE n0,n;
218:
219: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
220: vl = 0;
221: else {
222: t = BDY((Alg)a);
223: switch ( OID(t) ) {
224: case O_P: case O_R:
225: clctvr(ALG,t,&vl); break;
226: default:
227: vl = 0; break;
228: }
229: }
230: for ( n0 = 0; vl; vl = NEXT(vl) ) {
231: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
232: }
233: if ( n0 )
234: NEXT(n) = 0;
235: MKLIST(*rp,n0);
236: }
237:
238: void Pgetalgtree(arg,rp)
239: NODE arg;
240: LIST *rp;
241: {
242: Obj t;
243: P p;
244: VL vl,vl1,vl2;
245: Num a;
246: Alg b;
247: NODE n0,n;
248:
249: if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
250: vl = 0;
251: else {
252: t = BDY((Alg)a);
253: switch ( OID(t) ) {
254: case O_P:
1.5 noro 255: clctalg((P)t,&vl); break;
1.1 noro 256: case O_R:
257: clctalg(NM((R)t),&vl1);
258: clctalg(DN((R)t),&vl2);
259: mergev(ALG,vl1,vl2,&vl); break;
260: default:
261: vl = 0; break;
262: }
263: }
264: for ( n0 = 0; vl; vl = NEXT(vl) ) {
265: NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
266: }
267: if ( n0 )
268: NEXT(n) = 0;
269: MKLIST(*rp,n0);
270: }
271:
272: void clctalg(p,vl)
273: P p;
274: VL *vl;
275: {
276: int n,i;
277: VL tvl;
278: VN vn,vn1;
279: P d;
280: DCP dc;
281:
282: for ( n = 0, tvl = ALG; tvl; tvl = NEXT(tvl), n++ );
283: vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
284: for ( i = n-1, tvl = ALG; tvl; tvl = NEXT(tvl), i-- ) {
285: vn[i].v = tvl->v;
286: vn[i].n = 0;
287: }
288: markv(vn,n,p);
289: for ( i = n-1; i >= 0; i-- ) {
290: if ( !vn[i].n )
291: continue;
292: d = (P)vn[i].v->attr;
293: for ( dc = DC(d); dc; dc = NEXT(dc) )
294: markv(vn,i,COEF(dc));
295: }
296: vn1 = (VN) ALLOCA((n+1)*sizeof(struct oVN));
297: for ( i = 0; i < n; i++ ) {
298: vn1[i].v = vn[n-1-i].v; vn1[i].n = vn[n-1-i].n;
299: }
300: vntovl(vn1,n,vl);
301: }
302:
303: void Palg(arg,rp)
304: NODE arg;
305: Alg *rp;
306: {
307: Q a;
308: VL vl;
309: P x;
310: int n;
311:
312: a = (Q)ARG0(arg);
313: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
314: *rp = 0;
315: else {
316: n = ACNT-QTOS(a)-1;
317: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
318: if ( vl ) {
319: MKV(vl->v,x); MKAlg(x,*rp);
320: } else
321: *rp = 0;
322: }
323: }
324:
325: void Palgv(arg,rp)
326: NODE arg;
327: Obj *rp;
328: {
329: Q a;
330: VL vl;
331: P x;
332: int n;
333: Alg b;
334:
335: a = (Q)ARG0(arg);
336: if ( a && (OID(a) != O_N || NID(a) != N_Q || !INT(a)) )
337: *rp = 0;
338: else {
339: n = ACNT-QTOS(a)-1;
340: for ( vl = ALG; vl && n; vl = NEXT(vl), n-- );
341: if ( vl ) {
342: MKV(vl->v,x); MKAlg(x,b); algtorat((Num)b,rp);
343: } else
344: *rp = 0;
345: }
346: }
347:
348: void algptop(p,r)
349: P p,*r;
350: {
351: DCP dc,dcr,dcr0;
352:
353: if ( NUM(p) )
354: *r = (P)p;
355: else {
356: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
357: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
358: algptop(COEF(dc),&COEF(dcr));
359: }
360: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
361: }
362: }
363:
364: void algtorat(n,r)
365: Num n;
366: Obj *r;
367: {
368: Obj obj;
369: P nm,dn;
370:
371: if ( !n || NID(n) <= N_R )
372: *r = (Obj)n;
373: else {
374: obj = BDY((Alg)n);
375: if ( ID(obj) <= O_P )
376: algptop((P)obj,(P *)r);
377: else {
378: algptop(NM((R)obj),&nm); algptop(DN((R)obj),&dn);
379: divr(CO,(Obj)nm,(Obj)dn,r);
380: }
381: }
382: }
383:
384: void rattoalg(obj,n)
385: Obj obj;
386: Alg *n;
387: {
388: P nm,dn;
389: Obj t;
390:
391: if ( !obj || ID(obj) == O_N )
392: *n = (Alg)obj;
393: else if ( ID(obj) == O_P ) {
394: ptoalgp((P)obj,(P *)&t); MKAlg(t,*n);
395: } else {
396: ptoalgp(NM((R)obj),&nm); ptoalgp(DN((R)obj),&dn);
397: divr(ALG,(Obj)nm,(Obj)dn,&t); MKAlg(t,*n);
398: }
399: }
400:
401: void ptoalgp(p,r)
402: P p,*r;
403: {
404: DCP dc,dcr,dcr0;
405:
406: if ( NUM(p) )
407: *r = (P)p;
408: else {
409: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
410: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
411: ptoalgp(COEF(dc),&COEF(dcr));
412: }
413: NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
414: }
1.6 ! noro 415: }
! 416:
! 417: void invalg_le(Alg a,LIST *r);
! 418:
! 419: void Pinvalg_le(NODE arg,LIST *r)
! 420: {
! 421: invalg_le((Alg)ARG0(arg),r);
! 422: }
! 423:
! 424: typedef struct oMono_nf {
! 425: DP mono;
! 426: DP nf;
! 427: Q dn;
! 428: } *Mono_nf;
! 429:
! 430: void invalg_le(Alg a,LIST *r)
! 431: {
! 432: Alg inv;
! 433: MAT mobj,sol;
! 434: int *rinfo,*cinfo;
! 435: P p,dn,dn1,ap;
! 436: VL vl,tvl;
! 437: Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol;
! 438: int i,j,n,len,k;
! 439: MP mp,mp0;
! 440: DP dp,nm,nm1,m,d,u,u1;
! 441: NODE b,b1,hlist,mblist,t,s,rev0,rev,hist;
! 442: DP *ps;
! 443: struct order_spec *spec;
! 444: Mono_nf h,h1;
! 445: N nq,nr,nl,ng;
! 446: Q **mat,**solmat;
! 447: Q *w;
! 448: int *wi;
! 449:
! 450: ap = (P)BDY(a);
! 451: asir_assert(ap,O_P,"invalg_le");
! 452:
! 453: /* collecting algebraic numbers */
! 454: clctalg(ap,&vl);
! 455:
! 456: /* setup */
! 457: ptozp(ap,1,&c,&p);
! 458: STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec);
! 459: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 460: ps = (DP *)ALLOCA(n*sizeof(DP));
! 461:
! 462: /* conversion to DP */
! 463: for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) {
! 464: ptod(ALG,vl,tvl->v->attr,&ps[i]);
! 465: }
! 466: ptod(ALG,vl,p,&dp);
! 467: /* index list */
! 468: for ( b = 0, i = 0; i < n; i++ ) {
! 469: STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
! 470: }
! 471: /* simplification */
! 472: dp_true_nf(b,dp,ps,1,&nm,&dn);
! 473:
! 474: /* construction of NF table */
! 475:
! 476: /* stdmono: <<0,...,0>> < ... < max */
! 477: for ( hlist = 0, i = 0; i < n; i++ ) {
! 478: MKNODE(b1,(pointer)ps[i],hlist); hlist = b1;
! 479: }
! 480: dp_mbase(hlist,&rev0);
! 481: for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) {
! 482: MKNODE(b1,BDY(rev),mblist); mblist = b1;
! 483: }
! 484: dn0 = ONE;
! 485: for ( hist = 0, t = mblist; t; t = NEXT(t) ) {
! 486: /* searching a predecessor */
! 487: for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) {
! 488: h = (Mono_nf)BDY(s);
! 489: if ( dp_redble(m,h->mono) )
! 490: break;
! 491: }
! 492: h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf));
! 493: if ( s ) {
! 494: dp_subd(m,h->mono,&d);
! 495: muld(CO,d,h->nf,&u);
! 496: dp_true_nf(b,u,ps,1,&nm1,&dn1);
! 497: mulq(h->dn,(Q)dn1,&h1->dn);
! 498: } else {
! 499: muld(CO,m,nm,&u);
! 500: dp_true_nf(b,u,ps,1,&nm1,&dn1);
! 501: h1->dn = (Q)dn1;
! 502: }
! 503: h1->mono = m;
! 504: h1->nf = nm1;
! 505: MKNODE(b1,(pointer)h1,hist); hist = b1;
! 506:
! 507: /* dn0 = LCM(dn0,h1->dn) */
! 508: gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr);
! 509: muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0);
! 510: }
! 511: /* create a matrix */
! 512: len = length(mblist);
! 513: MKMAT(mobj,len,len+1);
! 514: mat = (Q **)BDY(mobj);
! 515: mat[len-1][len] = dn0;
! 516: for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) {
! 517: h = (Mono_nf)BDY(t);
! 518: nm1 = h->nf;
! 519: divq((Q)dn0,h->dn,&mul);
! 520: for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) )
! 521: if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) {
! 522: mulq(mul,(Q)mp->c,&mat[i][j]);
! 523: mp = NEXT(mp);
! 524: }
! 525: }
! 526: #if 0
! 527: w = (Q *)ALLOCA((len+1)*sizeof(Q));
! 528: wi = (int *)ALLOCA((len+1)*sizeof(int));
! 529: for ( i = 0; i < len; i++ ) {
! 530: for ( j = 0, k = 0; j <= len; j++ )
! 531: if ( mat[i][j] ) {
! 532: w[k] = mat[i][j];
! 533: wi[k] = j;
! 534: k++;
! 535: }
! 536: removecont_array(w,k);
! 537: for ( j = 0; j < k; j++ )
! 538: mat[i][wi[j]] = w[j];
! 539: }
! 540: #endif
! 541: generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
! 542: solmat = (Q **)BDY(sol);
! 543: for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) )
! 544: if ( solmat[i][0] ) {
! 545: NEXTMP(mp0,mp);
! 546: mp->c = (P)solmat[i][0];
! 547: mp->dl = BDY((DP)BDY(t))->dl;
! 548: }
! 549: NEXT(mp) = 0; MKDP(n,mp0,u);
! 550: dp_ptozp(u,&u1);
! 551: divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont);
! 552: dtop(ALG,vl,u1,&ap);
! 553: MKAlg(ap,inv);
! 554: mulq(dnsol,(Q)dn,&c1);
! 555: mulq(c1,c,&c2);
! 556: divq(c2,cont,&c3);
! 557: b = mknode(2,inv,c3);
! 558: MKLIST(*r,b);
1.1 noro 559: }