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