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