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