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