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