![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to fctr.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 fctr.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/fctr.c,v 1.10 2001/11/19 00:57:10 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: void Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint();
54: void Pptozp(), Pcont();
55: void Pafctr(), Pagcd();
56: void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab();
1.11 ! noro 57: void Psfsqfr(),Psfbfctr(),Psfufctr(),Psfmintdeg();
1.1 noro 58: void Pirred_check(), Pnfctr_mod();
59:
1.11 ! noro 60: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr);
! 61: void create_bmono(P c,V x,int i,V y,int j,P *mono);
! 62:
1.1 noro 63: struct ftab fctr_tab[] = {
1.5 noro 64: {"fctr",Pfctr,-2},
1.1 noro 65: {"gcd",Pgcd,-3},
66: {"gcdz",Pgcdz,2},
67: {"lcm",Plcm,2},
68: {"sqfr",Psqfr,1},
69: {"ufctrhint",Pufctrhint,2},
70: {"ptozp",Pptozp,1},
71: {"cont",Pcont,-2},
72: {"afctr",Pafctr,2},
73: {"agcd",Pagcd,3},
74: {"modsqfr",Pmodsqfr,2},
75: {"modfctr",Pmodfctr,2},
1.10 noro 76: {"sfsqfr",Psfsqfr,1},
77: {"sfufctr",Psfufctr,1},
78: {"sfbfctr",Psfbfctr,-4},
1.11 ! noro 79: {"sfmintdeg",Psfmintdeg,5},
1.1 noro 80: #if 0
81: {"ddd",Pddd,2},
82: {"newddd",Pnewddd,2},
83: #endif
84: {"ddd_tab",Pddd_tab,2},
85: {"irred_check",Pirred_check,2},
86: {"nfctr_mod",Pnfctr_mod,2},
87: {0,0,0},
88: };
89:
90: void Pfctr(arg,rp)
91: NODE arg;
92: LIST *rp;
93: {
94: DCP dc;
95:
96: asir_assert(ARG0(arg),O_P,"fctr");
1.5 noro 97: if ( argc(arg) == 1 )
98: fctrp(CO,(P)ARG0(arg),&dc);
99: else {
100: asir_assert(ARG1(arg),O_P,"fctr");
101: fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc);
102: }
1.1 noro 103: dcptolist(dc,rp);
104: }
105:
106: void Pgcd(arg,rp)
107: NODE arg;
108: P *rp;
109: {
110: P p1,p2,g1,g2,g;
111: Num m;
112: int mod;
113:
114: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
115: asir_assert(p1,O_P,"gcd");
116: asir_assert(p2,O_P,"gcd");
117: if ( !p1 )
118: *rp = p2;
119: else if ( !p2 )
120: *rp = p1;
121: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
1.4 noro 122: gcdprsp(CO,p1,p2,rp);
1.1 noro 123: else if ( argc(arg) == 2 )
124: ezgcdp(CO,p1,p2,rp);
125: else {
126: m = (Num)ARG2(arg);
127: asir_assert(m,O_P,"gcd");
128: mod = QTOS((Q)m);
129: ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
130: gcdprsmp(CO,mod,g1,g2,&g);
131: mptop(g,rp);
132: }
133: }
134:
135: void Pgcdz(arg,rp)
136: NODE arg;
137: P *rp;
138: {
139: P p1,p2,t;
140: Q c1,c2;
141: N n;
142:
143: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
144: asir_assert(p1,O_P,"gcdz");
145: asir_assert(p2,O_P,"gcdz");
146: if ( !p1 )
147: *rp = p2;
148: else if ( !p2 )
149: *rp = p1;
150: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
151: error("gcdz : invalid argument");
152: else if ( NUM(p1) || NUM(p2) ) {
153: if ( NUM(p1) )
154: c1 = (Q)p1;
155: else
156: ptozp(p1,1,&c1,&t);
157: if ( NUM(p2) )
158: c2 = (Q)p2;
159: else
160: ptozp(p2,1,&c2,&t);
161: gcdn(NM(c1),NM(c2),&n); NTOQ(n,1,c1); *rp = (P)c1;
162: } else {
163: #if 0
164: w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
165: #endif
166: ezgcdpz(CO,p1,p2,rp);
167: }
168: }
169:
170: void Plcm(arg,rp)
171: NODE arg;
172: P *rp;
173: {
174: P t1,t2,p1,p2,g,q;
175: Q c;
176:
177: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
178: asir_assert(p1,O_P,"lcm");
179: asir_assert(p2,O_P,"lcm");
180: if ( !p1 || !p2 )
181: *rp = 0;
182: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
183: error("lcm : invalid argument");
184: else {
185: ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
186: ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
187: }
188: }
189:
190: void Psqfr(arg,rp)
191: NODE arg;
192: LIST *rp;
193: {
194: DCP dc;
195:
196: asir_assert(ARG0(arg),O_P,"sqfr");
197: sqfrp(CO,(P)ARG0(arg),&dc);
198: dcptolist(dc,rp);
199: }
200:
201: void Pufctrhint(arg,rp)
202: NODE arg;
203: LIST *rp;
204: {
205: DCP dc;
206:
207: asir_assert(ARG0(arg),O_P,"ufctrhint");
208: asir_assert(ARG1(arg),O_N,"ufctrhint");
209: ufctr((P)ARG0(arg),QTOS((Q)ARG1(arg)),&dc);
210: dcptolist(dc,rp);
211: }
212:
213: #if 0
214: Pmgcd(arg,rp)
215: NODE arg;
216: Obj *rp;
217: {
218: NODE node,tn;
219: int i,m;
220: P *l;
221:
222: node = BDY((LIST)ARG0(arg));
223: for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
224: m = i; l = (P *)ALLOCA(m*sizeof(P));
225: for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
226: l[i] = (P)BDY(tn);
227: nezgcdnpz(CO,l,m,rp);
228: }
229: #endif
230:
231: void Pcont(arg,rp)
232: NODE arg;
233: P *rp;
234: {
235: DCP dc;
236: int m;
237: P p,p1;
238: P *l;
239: V v;
240:
241: asir_assert(ARG0(arg),O_P,"cont");
242: p = (P)ARG0(arg);
243: if ( NUM(p) )
244: *rp = p;
245: else {
246: if ( argc(arg) == 2 ) {
247: v = VR((P)ARG1(arg));
248: change_mvar(CO,p,v,&p1);
249: if ( VR(p1) != v ) {
250: *rp = p1; return;
251: } else
252: p = p1;
253: }
254: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
255: l = (P *)ALLOCA(m*sizeof(P));
256: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
257: l[m] = COEF(dc);
258: nezgcdnpz(CO,l,m,rp);
259: }
260: }
261:
262: void Pptozp(arg,rp)
263: NODE arg;
264: P *rp;
265: {
266: Q t;
267:
268: asir_assert(ARG0(arg),O_P,"ptozp");
269: ptozp((P)ARG0(arg),1,&t,rp);
270: }
271:
272: void Pafctr(arg,rp)
273: NODE arg;
274: LIST *rp;
275: {
276: DCP dc;
277:
278: asir_assert(ARG0(arg),O_P,"afctr");
279: asir_assert(ARG1(arg),O_P,"afctr");
280: afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
281: dcptolist(dc,rp);
282: }
283:
284: void Pagcd(arg,rp)
285: NODE arg;
286: P *rp;
287: {
288: asir_assert(ARG0(arg),O_P,"agcd");
289: asir_assert(ARG1(arg),O_P,"agcd");
290: asir_assert(ARG2(arg),O_P,"agcd");
291: gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
292: }
293:
294: #if 1
295: #define Mulum mulum
296: #define Divum divum
297: #define Mulsum mulsum
298: #define Gcdum gcdum
299: #endif
300:
301: void Mulum(), Mulsum(), Gcdum();
302: int Divum();
303:
304: #define FCTR 0 /* berlekamp */
305: #define SQFR 1
306: #define DDD 2 /* Cantor-Zassenhauss */
307: #define NEWDDD 3 /* berlekamp + root-finding by Cantor-Zassenhauss */
308:
309: UM *resberle();
310:
311: void Pmodfctr(arg,rp)
312: NODE arg;
313: LIST *rp;
314: {
315: DCP dc;
316: int mod;
317:
318: mod = QTOS((Q)ARG1(arg));
319: if ( mod < 0 )
320: error("modfctr : invalid modulus");
321: modfctrp(ARG0(arg),mod,NEWDDD,&dc);
1.6 noro 322: if ( !dc ) {
323: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
324: }
325: dcptolist(dc,rp);
326: }
327:
1.10 noro 328: void Psfsqfr(arg,rp)
329: NODE arg;
330: LIST *rp;
331: {
332: DCP dc;
333:
334: sfsqfr(ARG0(arg),&dc);
335: dcptolist(dc,rp);
336: }
337:
338: void Psfufctr(arg,rp)
1.6 noro 339: NODE arg;
340: LIST *rp;
341: {
342: DCP dc;
343:
344: fctrsf(ARG0(arg),&dc);
1.1 noro 345: dcptolist(dc,rp);
1.7 noro 346: }
347:
348: void Psfbfctr(arg,rp)
349: NODE arg;
350: LIST *rp;
351: {
352: V x,y;
1.8 noro 353: DCP dc,dct;
354: P t;
355: struct oVL vl1,vl2;
356: VL vl;
1.10 noro 357: int degbound;
1.7 noro 358:
359: x = VR((P)ARG1(arg));
360: y = VR((P)ARG2(arg));
1.8 noro 361: vl1.v = x; vl1.next = &vl2;
362: vl2.v = y; vl2.next = 0;
363: vl = &vl1;
1.10 noro 364: if ( argc(arg) == 4 )
365: degbound = QTOS((Q)ARG3(arg));
366: else
367: degbound = -1;
1.8 noro 368:
1.10 noro 369: sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
1.8 noro 370: for ( dct = dc; dct; dct = NEXT(dct) ) {
371: reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
1.7 noro 372: }
1.8 noro 373: dcptolist(dc,rp);
1.1 noro 374: }
375:
1.11 ! noro 376: void Psfmintdeg(arg,rp)
! 377: NODE arg;
! 378: P *rp;
! 379: {
! 380: V x,y;
! 381: P r;
! 382: struct oVL vl1,vl2;
! 383: VL vl;
! 384: int dy,c;
! 385:
! 386: x = VR((P)ARG1(arg));
! 387: y = VR((P)ARG2(arg));
! 388: vl1.v = x; vl1.next = &vl2;
! 389: vl2.v = y; vl2.next = 0;
! 390: vl = &vl1;
! 391: dy = QTOS((Q)ARG3(arg));
! 392: c = QTOS((Q)ARG4(arg));
! 393: sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
! 394: reorderp(CO,vl,r,rp);
! 395: }
! 396:
1.1 noro 397: void Pmodsqfr(arg,rp)
398: NODE arg;
399: LIST *rp;
400: {
401: DCP dc;
402:
1.9 noro 403: if ( !ARG0(arg) ) {
1.1 noro 404: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 405: } else
406: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc);
1.1 noro 407: dcptolist(dc,rp);
408: }
409:
410: void Pddd(arg,rp)
411: NODE arg;
412: LIST *rp;
413: {
414: DCP dc;
415:
1.9 noro 416: if ( !ARG0(arg) ) {
1.1 noro 417: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 418: } else
419: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc);
1.1 noro 420: dcptolist(dc,rp);
421: }
422:
423: void Pnewddd(arg,rp)
424: NODE arg;
425: LIST *rp;
426: {
1.9 noro 427: DCP dc=0;
1.1 noro 428:
1.9 noro 429: if ( !ARG0(arg) ) {
1.1 noro 430: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 431: } else
432: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc);
1.1 noro 433: dcptolist(dc,rp);
434: }
435:
436: void Pirred_check(arg,rp)
437: NODE arg;
438: Q *rp;
439: {
440: P p;
441: UM mp;
442: int r,mod;
443:
444: p = (P)ARG0(arg);
445: if ( !p ) {
446: *rp = 0; return;
447: }
448: mp = W_UMALLOC(UDEG(p));
449: mod = QTOS((Q)ARG1(arg));
450: ptoum(mod,p,mp);
451: r = irred_check(mp,mod);
452: if ( r )
453: *rp = ONE;
454: else
455: *rp = 0;
456: }
457:
458: void Pnfctr_mod(arg,rp)
459: NODE arg;
460: Q *rp;
461: {
462: P p;
463: UM mp;
464: int r,mod;
465:
466: p = (P)ARG0(arg);
467: if ( !p ) {
468: *rp = 0; return;
469: }
470: mp = W_UMALLOC(UDEG(p));
471: mod = QTOS((Q)ARG1(arg));
472: ptoum(mod,p,mp);
473: r = nfctr_mod(mp,mod);
474: STOQ(r,*rp);
475: }
476:
477: void Pddd_tab(arg,rp)
478: NODE arg;
479: VECT *rp;
480: {
481: P p;
482: UM mp,t,q,r1,w,w1;
483: UM *r,*s;
484: int dr,mod,n,i;
485: VECT result;
486: V v;
487:
488: p = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
489: v = VR(p);
490: n = UDEG(p); mp = W_UMALLOC(n);
491: ptoum(mod,p,mp);
492: r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
493: r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
494: t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
495: DEG(t) = mod; COEF(t)[mod] = 1;
496: q = W_UMALLOC(mod);
497: dr = divum(mod,t,mp,q);
498: DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
499: s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
500: w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
501: w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
502: for ( i = 1; i < n; i++ ) {
503: DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
504: mulum(mod,r1,w,w1);
505: dr = divum(mod,w1,mp,q); DEG(w1) = dr;
506: s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
507: }
508: for ( i = 2; i < n; i++ ) {
509: mult_mod_tab(r[i-1],mod,s,w,n);
510: r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
511: }
512: MKVECT(result,n);
513: for ( i = 0; i < n; i++ )
514: umtop(v,r[i],(P *)&BDY(result)[i]);
515: *rp = result;
1.11 ! noro 516: }
! 517:
! 518: struct lb {
! 519: int pos,len;
! 520: int *r;
! 521: int *hist;
! 522: };
! 523:
! 524: static NODE insert_lb(NODE g,struct lb *a)
! 525: {
! 526: NODE prev,cur,n;
! 527:
! 528: prev = 0; cur = g;
! 529: while ( cur ) {
! 530: if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
! 531: MKNODE(n,a,cur);
! 532: if ( !prev )
! 533: return n;
! 534: else {
! 535: NEXT(prev) = n;
! 536: return g;
! 537: }
! 538: } else {
! 539: prev = cur;
! 540: cur = NEXT(cur);
! 541: }
! 542: }
! 543: MKNODE(n,a,0);
! 544: NEXT(prev) = n;
! 545: return g;
! 546: }
! 547:
! 548: static void lnf(int *r,int *h,int n,int len,NODE g)
! 549: {
! 550: struct lb *t;
! 551: int pos,i,j,len1,c;
! 552: int *r1,*h1;
! 553:
! 554: for ( ; g; g = NEXT(g) ) {
! 555: t = (struct lb *)BDY(g);
! 556: pos = t->pos;
! 557: if ( c = r[pos] ) {
! 558: r1 = t->r;
! 559: h1 = t->hist;
! 560: len1 = t->len;
! 561: for ( i = pos; i < n; i++ )
! 562: r[i] = _subsf(r[i],_mulsf(r1[i],c));
! 563: for ( i = 0; i < len1; i++ )
! 564: h[i] = _subsf(h[i],_mulsf(h1[i],c));
! 565: }
! 566: }
! 567: for ( i = 0; i < n && !r[i]; i++ );
! 568: if ( i < n ) {
! 569: c = _invsf(r[i]);
! 570: for ( j = i; j < n; j++ )
! 571: r[j] = _mulsf(r[j],c);
! 572: for ( j = i; j < len; j++ )
! 573: h[j] = _mulsf(h[j],c);
! 574: }
! 575: }
! 576:
! 577: void print_vect(int *r,int len)
! 578: {
! 579: int i;
! 580:
! 581: for ( i = 0; i < len; i++ )
! 582: if ( r[i] ) printf("(%d %d)",i,IFTOF(r[i]));
! 583: printf("\n");
! 584: }
! 585:
! 586: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
! 587: {
! 588: V x,y;
! 589: int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
! 590: UP *rx;
! 591: DCP dc;
! 592: P t,f,mono,f1;
! 593: UP ut,h;
! 594: int ***nf;
! 595: int *r,*hist,*prev,*r1;
! 596: struct lb *lb;
! 597: GFS s;
! 598: NODE g;
! 599:
! 600: x = vl->v;
! 601: y = NEXT(vl)->v;
! 602: dx = getdeg(x,fx);
! 603: dxdy = dx*dy;
! 604: /* rx = -(fx-x^dx) */
! 605: rx = (UP *)CALLOC(dx,sizeof(UP));
! 606: for ( dc = DC(fx); dc; dc = NEXT(dc)) {
! 607: chsgnp(COEF(dc),&t);
! 608: ptoup(t,&ut);
! 609: rx[QTOS(DEG(dc))] = ut;
! 610: }
! 611: /* nf[d] = normal form table of monomials with total degree d */
! 612: nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
! 613: nf[0] = (int **)CALLOC(1,sizeof(int *));
! 614:
! 615: /* nf[0][0] = 1 */
! 616: r = (int *)CALLOC(dxdy,sizeof(int));
! 617: r[0] = _onesf();
! 618: nf[0][0] = r;
! 619:
! 620: hist = (int *)CALLOC(1,sizeof(int));
! 621: r[0] = _onesf();
! 622:
! 623: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
! 624: lb->pos = 0;
! 625: lb->r = r;
! 626: lb->hist = hist;
! 627: lb->len = 1;
! 628:
! 629: /* g : table of normal form as linear form */
! 630: MKNODE(g,lb,0);
! 631:
! 632: len = 1;
! 633: h = UPALLOC(dy);
! 634: for ( d = 1; ; d++ ) {
! 635: if ( d > c ){
! 636: return;
! 637: }
! 638: nf[d] = (int **)CALLOC(d+1,sizeof(int *));
! 639: len0 = len;
! 640: len += d+1;
! 641:
! 642: for ( i = d; i >= 0; i-- ) {
! 643: /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
! 644: /* nf(x^d) = nf(nf(x^(d-1))*x) */
! 645: r = (int *)CALLOC(dxdy,sizeof(int));
! 646: if ( i == 0 ) {
! 647: prev = nf[d-1][0];
! 648: bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
! 649:
! 650: /* create the head coeff */
! 651: for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
! 652: if ( prev[k] ) {
! 653: u = IFTOF(prev[k]);
! 654: MKGFS(u,s);
! 655: } else
! 656: s = 0;
! 657: COEF(h)[l] = (Num)s;
! 658: }
! 659: for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
! 660: DEG(h) = l;
! 661:
! 662: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
! 663: tmulup(rx[k],h,dy,&ut);
! 664: if ( ut )
! 665: for ( l = 0; l < dy; l++ ) {
! 666: s = (GFS)COEF(ut)[l];
! 667: if ( s ) {
! 668: u = CONT(s);
! 669: r[dyk+l] = _addsf(r[dyk+l],FTOIF(u));
! 670: }
! 671: }
! 672: }
! 673: } else {
! 674: prev = nf[d-1][i-1];
! 675: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
! 676: for ( l = 1; l < dy; l++ )
! 677: r[dyk+l] = prev[dyk+l-1];
! 678: }
! 679: }
! 680: nf[d][i] = r;
! 681: hist = (int *)CALLOC(len,sizeof(int));
! 682: hist[len0+i] = _onesf();
! 683: r1 = (int *)CALLOC(dxdy,sizeof(int));
! 684: bcopy(r,r1,dxdy*sizeof(int));
! 685: lnf(r1,hist,dxdy,len,g);
! 686: for ( k = 0; k < dxdy && !r1[k]; k++ );
! 687: if ( k == dxdy ) {
! 688: f = 0;
! 689: for ( k = j = 0; k <= d; k++ )
! 690: for ( i = 0; i <= k; i++, j++ )
! 691: if ( hist[j] ) {
! 692: u = IFTOF(hist[j]);
! 693: MKGFS(u,s);
! 694: /* mono = s*x^(k-i)*y^i */
! 695: create_bmono((P)s,x,k-i,y,i,&mono);
! 696: addp(vl,f,mono,&f1);
! 697: f = f1;
! 698: }
! 699: *fr = f;
! 700: return;
! 701: } else {
! 702: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
! 703: lb->pos = k;
! 704: lb->r = r1;
! 705: lb->hist = hist;
! 706: lb->len = len;
! 707: g = insert_lb(g,lb);
! 708: }
! 709: }
! 710: }
! 711: }
! 712:
! 713: void create_bmono(P c,V x,int i,V y,int j,P *mono)
! 714: {
! 715: P t,s;
! 716:
! 717: if ( !i )
! 718: if ( !j )
! 719: t = c;
! 720: else {
! 721: /* c*y^j */
! 722: MKV(y,t);
! 723: COEF(DC(t)) = c;
! 724: STOQ(j,DEG(DC(t)));
! 725: }
! 726: else if ( !j ) {
! 727: /* c*x^i */
! 728: MKV(x,t);
! 729: COEF(DC(t)) = c;
! 730: STOQ(i,DEG(DC(t)));
! 731: } else {
! 732: MKV(y,s);
! 733: COEF(DC(s)) = c;
! 734: STOQ(j,DEG(DC(s)));
! 735: MKV(x,t);
! 736: COEF(DC(t)) = s;
! 737: STOQ(i,DEG(DC(t)));
! 738: }
! 739: *mono = t;
1.1 noro 740: }