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