Annotation of OpenXM_contrib2/asir2018/builtin/fctr.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/fctr.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 Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint();
54: void Pptozp(), Pcont(), Psfcont();
55: void Pafctr(), Pagcd();
56: void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab();
57: void Psfsqfr(),Psffctr(),Psfbfctr(),Psfufctr(),Psfmintdeg(),Psfgcd();
58: void Pirred_check(), Pnfctr_mod();
59: void Pbivariate_hensel_special();
60:
61: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr);
62:
63: struct ftab fctr_tab[] = {
64: {"bivariate_hensel_special",Pbivariate_hensel_special,6},
65: {"fctr",Pfctr,-2},
66: {"gcd",Pgcd,-3},
67: {"gcdz",Pgcdz,2},
68: {"lcm",Plcm,2},
69: {"sqfr",Psqfr,1},
70: {"ufctrhint",Pufctrhint,2},
71: {"ptozp",Pptozp,1},
72: {"cont",Pcont,-2},
73: {"sfcont",Psfcont,-2},
74: {"afctr",Pafctr,2},
75: {"agcd",Pagcd,3},
76: {"modsqfr",Pmodsqfr,2},
77: {"modfctr",Pmodfctr,2},
78: {"sfsqfr",Psfsqfr,1},
79: {"sffctr",Psffctr,1},
80: {"sfufctr",Psfufctr,1},
81: {"sfbfctr",Psfbfctr,-4},
82: {"sfmintdeg",Psfmintdeg,5},
83: {"sfgcd",Psfgcd,2},
84: #if 0
85: {"ddd",Pddd,2},
86: {"newddd",Pnewddd,2},
87: #endif
88: {"ddd_tab",Pddd_tab,2},
89: {"irred_check",Pirred_check,2},
90: {"nfctr_mod",Pnfctr_mod,2},
91: {0,0,0},
92: };
93:
94: /* bivariate_hensel_special(f(x,y):monic in x,g0(x),h0(y),x,y,d) */
95:
96: void Pbivariate_hensel_special(NODE arg,LIST *rp)
97: {
98: DCP dc;
99: struct oVN vn[2];
100: P f,g0,h0,ak,bk,gk,hk;
101: V vx,vy;
102: VL nvl;
103: Z qk;
104: Q cbd,bb;
105: int d;
106: NODE n;
107:
108: f = (P)ARG0(arg);
109: g0 = (P)ARG1(arg);
110: h0 = (P)ARG2(arg);
111: vx = VR((P)ARG3(arg));
112: vy = VR((P)ARG4(arg));
1.2 ! noro 113: d = ZTOS((Q)ARG5(arg));
1.1 noro 114: NEWVL(nvl); nvl->v = vx;
115: NEWVL(NEXT(nvl)); NEXT(nvl)->v = vy;
116: NEXT(NEXT(nvl)) = 0;
117: vn[0].v = vy; vn[0].n = 0;
118: vn[1].v = 0; vn[1].n = 0;
119: cbound(nvl,f,&cbd);
120: addq(cbd,cbd,&bb);
121: henzq1(g0,h0,bb,&bk,&ak,&qk);
122: henmv(nvl,vn,f,g0,h0,ak,bk,(P)ONE,(P)ONE,(P)ONE,(P)ONE,qk,d,&gk,&hk);
123: n = mknode(2,gk,hk);
124: MKLIST(*rp,n);
125: }
126:
127: void Pfctr(NODE arg,LIST *rp)
128: {
129: DCP dc;
130:
131: asir_assert(ARG0(arg),O_P,"fctr");
132: if ( argc(arg) == 1 )
133: fctrp(CO,(P)ARG0(arg),&dc);
134: else {
135: asir_assert(ARG1(arg),O_P,"fctr");
136: fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc);
137: }
138: dcptolist(dc,rp);
139: }
140:
141: void Pgcd(NODE arg,P *rp)
142: {
143: P p1,p2,g1,g2,g;
144: Num m;
145: int mod;
146:
147: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
148: asir_assert(p1,O_P,"gcd");
149: asir_assert(p2,O_P,"gcd");
150: if ( !p1 )
151: *rp = p2;
152: else if ( !p2 )
153: *rp = p1;
154: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
155: gcdprsp(CO,p1,p2,rp);
156: else if ( argc(arg) == 2 )
157: ezgcdp(CO,p1,p2,rp);
158: else {
159: m = (Num)ARG2(arg);
160: asir_assert(m,O_P,"gcd");
1.2 ! noro 161: mod = ZTOS((Q)m);
1.1 noro 162: ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
163: gcdprsmp(CO,mod,g1,g2,&g);
164: mptop(g,rp);
165: }
166: }
167:
168: void Pgcdz(NODE arg,P *rp)
169: {
170: P p1,p2,t;
171: Q c1,c2;
172: Z n1,n2,n;
173:
174: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
175: asir_assert(p1,O_P,"gcdz");
176: asir_assert(p2,O_P,"gcdz");
177: if ( !p1 )
178: *rp = p2;
179: else if ( !p2 )
180: *rp = p1;
181: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
182: error("gcdz : invalid argument");
183: else if ( NUM(p1) || NUM(p2) ) {
184: if ( NUM(p1) )
185: c1 = (Q)p1;
186: else
187: ptozp(p1,1,&c1,&t);
188: if ( NUM(p2) )
189: c2 = (Q)p2;
190: else
191: ptozp(p2,1,&c2,&t);
192: /* XXX */
193: nmq(c1,&n1); nmq(c2,&n2);
194: gcdz(n1,n2,&n); *rp = (P)n;
195: } else {
196: #if 0
197: w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
198: #endif
199: ezgcdpz(CO,p1,p2,rp);
200: }
201: }
202:
203: void Plcm(NODE arg,P *rp)
204: {
205: P t1,t2,p1,p2,g,q;
206: Q c;
207:
208: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
209: asir_assert(p1,O_P,"lcm");
210: asir_assert(p2,O_P,"lcm");
211: if ( !p1 || !p2 )
212: *rp = 0;
213: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
214: error("lcm : invalid argument");
215: else {
216: ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
217: ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
218: }
219: }
220:
221: void Psqfr(NODE arg,LIST *rp)
222: {
223: DCP dc;
224:
225: asir_assert(ARG0(arg),O_P,"sqfr");
226: sqfrp(CO,(P)ARG0(arg),&dc);
227: dcptolist(dc,rp);
228: }
229:
230: void Pufctrhint(NODE arg,LIST *rp)
231: {
232: DCP dc;
233:
234: asir_assert(ARG0(arg),O_P,"ufctrhint");
235: asir_assert(ARG1(arg),O_N,"ufctrhint");
1.2 ! noro 236: ufctr((P)ARG0(arg),ZTOS((Q)ARG1(arg)),&dc);
1.1 noro 237: dcptolist(dc,rp);
238: }
239:
240: #if 0
241: Pmgcd(arg,rp)
242: NODE arg;
243: Obj *rp;
244: {
245: NODE node,tn;
246: int i,m;
247: P *l;
248:
249: node = BDY((LIST)ARG0(arg));
250: for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
251: m = i; l = (P *)ALLOCA(m*sizeof(P));
252: for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
253: l[i] = (P)BDY(tn);
254: nezgcdnpz(CO,l,m,rp);
255: }
256: #endif
257:
258: void Pcont(NODE arg,P *rp)
259: {
260: DCP dc;
261: int m;
262: P p,p1;
263: P *l;
264: V v;
265:
266: asir_assert(ARG0(arg),O_P,"cont");
267: p = (P)ARG0(arg);
268: if ( NUM(p) )
269: *rp = p;
270: else {
271: if ( argc(arg) == 2 ) {
272: v = VR((P)ARG1(arg));
273: change_mvar(CO,p,v,&p1);
274: if ( VR(p1) != v ) {
275: *rp = p1; return;
276: } else
277: p = p1;
278: }
279: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
280: l = (P *)ALLOCA(m*sizeof(P));
281: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
282: l[m] = COEF(dc);
283: nezgcdnpz(CO,l,m,rp);
284: }
285: }
286:
287: void Psfcont(NODE arg,P *rp)
288: {
289: DCP dc;
290: MP mp;
291: int m;
292: Obj obj;
293: P p,p1;
294: P *l;
295: V v;
296:
297: obj = (Obj)ARG0(arg);
298: if ( !obj || NUM(obj) )
299: *rp = (P)obj;
300: else if ( OID(obj) == O_P ) {
301: p = (P)obj;
302: if ( argc(arg) == 2 ) {
303: v = VR((P)ARG1(arg));
304: change_mvar(CO,p,v,&p1);
305: if ( VR(p1) != v ) {
306: *rp = p1; return;
307: } else
308: p = p1;
309: }
310: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
311: l = (P *)ALLOCA(m*sizeof(P));
312: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
313: l[m] = COEF(dc);
314: gcdsf(CO,l,m,rp);
315: } else if ( OID(obj) == O_DP ) {
316: for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++ );
317: l = (P *)ALLOCA(m*sizeof(P));
318: for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++)
319: l[m] = (P)mp->c;
320: gcdsf(CO,l,m,rp);
321: }
322: }
323:
324: void Pptozp(NODE arg,Obj *rp)
325: {
326: Q t;
327: NODE tt,p;
328: NODE n,n0;
329: char *key;
330: P pp;
331: LIST list;
332: int get_factor=0;
333:
334: asir_assert(ARG0(arg),O_P,"ptozp");
335:
336: /* analyze the option */
337: if ( current_option ) {
338: for ( tt = current_option; tt; tt = NEXT(tt) ) {
339: p = BDY((LIST)BDY(tt));
340: key = BDY((STRING)BDY(p));
341: /* value = (Obj)BDY(NEXT(p)); */
342: if ( !strcmp(key,"factor") ) get_factor=1;
343: else {
344: error("ptozp: unknown option.");
345: }
346: }
347: }
348:
349: ptozp((P)ARG0(arg),1,&t,&pp);
350:
351: /* printexpr(NULL,t); */
352: /* if the option factor is given, then it returns the answer
353: in the format [zpoly, num] where num*zpoly is equal to the argument.*/
354: if (get_factor) {
355: n0 = mknode(2,pp,t);
356: MKLIST(list,n0);
357: *rp = (Obj)list;
358: } else
359: *rp = (Obj)pp;
360: }
361:
362: void Pafctr(NODE arg,LIST *rp)
363: {
364: DCP dc;
365:
366: asir_assert(ARG0(arg),O_P,"afctr");
367: asir_assert(ARG1(arg),O_P,"afctr");
368: afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
369: dcptolist(dc,rp);
370: }
371:
372: void Pagcd(NODE arg,P *rp)
373: {
374: asir_assert(ARG0(arg),O_P,"agcd");
375: asir_assert(ARG1(arg),O_P,"agcd");
376: asir_assert(ARG2(arg),O_P,"agcd");
377: gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
378: }
379:
380: #if 1
381: #define Mulum mulum
382: #define Divum divum
383: #define Mulsum mulsum
384: #define Gcdum gcdum
385: #endif
386:
387: void Mulum(), Mulsum(), Gcdum();
388: int Divum();
389:
390: #define FCTR 0 /* berlekamp */
391: #define SQFR 1
392: #define DDD 2 /* Cantor-Zassenhauss */
393: #define NEWDDD 3 /* berlekamp + root-finding by Cantor-Zassenhauss */
394:
395: UM *resberle();
396:
397: void reduce_sfdc(DCP sfdc, DCP *dc);
398:
399: void Pmodfctr(NODE arg,LIST *rp)
400: {
401: DCP dc,dcu;
402: int mod,i,t;
403: P p;
404: Obj u;
405: VL vl;
406:
1.2 ! noro 407: mod = ZTOS((Q)ARG1(arg));
1.1 noro 408: if ( mod < 0 )
409: error("modfctr : invalid modulus");
410: p = (P)ARG0(arg);
411: clctv(CO,p,&vl);
412: if ( !vl ) {
413: NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE; NEXT(dc) = 0;
414: } else if ( !NEXT(vl) )
415: modfctrp(ARG0(arg),mod,NEWDDD,&dc);
416: else {
417: /* XXX 16384 should be replaced by a macro */
418: for ( i = 1, t = mod; t*mod < 16384; t *= mod, i++ );
419: current_ff = FF_GFS;
420: setmod_sf(mod,i);
421: simp_ff((Obj)p,&u);
422: mfctrsf(CO,(P)u,&dcu);
423: reduce_sfdc(dcu,&dc);
424: }
425: if ( !dc ) {
426: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
427: }
428: dcptolist(dc,rp);
429: }
430:
431: void Psfgcd(NODE arg,P *rp)
432: {
433: P ps[2];
434:
435: ps[0] = (P)ARG0(arg);
436: ps[1] = (P)ARG1(arg);
437: gcdsf(CO,ps,2,rp);
438: }
439:
440: void Psffctr(NODE arg,LIST *rp)
441: {
442: DCP dc;
443:
444: mfctrsf(CO,ARG0(arg),&dc);
445: dcptolist(dc,rp);
446: }
447:
448: void Psfsqfr(NODE arg,LIST *rp)
449: {
450: DCP dc;
451:
452: sqfrsf(CO,ARG0(arg),&dc);
453: dcptolist(dc,rp);
454: }
455:
456: void Psfufctr(NODE arg,LIST *rp)
457: {
458: DCP dc;
459:
460: ufctrsf(ARG0(arg),&dc);
461: dcptolist(dc,rp);
462: }
463:
464: void Psfbfctr(NODE arg,LIST *rp)
465: {
466: V x,y;
467: DCP dc,dct;
468: P t;
469: struct oVL vl1,vl2;
470: VL vl;
471: int degbound;
472:
473: x = VR((P)ARG1(arg));
474: y = VR((P)ARG2(arg));
475: vl1.v = x; vl1.next = &vl2;
476: vl2.v = y; vl2.next = 0;
477: vl = &vl1;
478: if ( argc(arg) == 4 )
1.2 ! noro 479: degbound = ZTOS((Q)ARG3(arg));
1.1 noro 480: else
481: degbound = -1;
482:
483: sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
484: for ( dct = dc; dct; dct = NEXT(dct) ) {
485: reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
486: }
487: dcptolist(dc,rp);
488: }
489:
490: void Psfmintdeg(NODE arg,P *rp)
491: {
492: V x,y;
493: P r;
494: struct oVL vl1,vl2;
495: VL vl;
496: int dy,c;
497:
498: x = VR((P)ARG1(arg));
499: y = VR((P)ARG2(arg));
500: vl1.v = x; vl1.next = &vl2;
501: vl2.v = y; vl2.next = 0;
502: vl = &vl1;
1.2 ! noro 503: dy = ZTOS((Q)ARG3(arg));
! 504: c = ZTOS((Q)ARG4(arg));
1.1 noro 505: sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
506: reorderp(CO,vl,r,rp);
507: }
508:
509: void Pmodsqfr(NODE arg,LIST *rp)
510: {
511: DCP dc;
512:
513: if ( !ARG0(arg) ) {
514: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
515: } else
1.2 ! noro 516: modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),SQFR,&dc);
1.1 noro 517: dcptolist(dc,rp);
518: }
519:
520: void Pddd(NODE arg,LIST *rp)
521: {
522: DCP dc;
523:
524: if ( !ARG0(arg) ) {
525: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
526: } else
1.2 ! noro 527: modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),DDD,&dc);
1.1 noro 528: dcptolist(dc,rp);
529: }
530:
531: void Pnewddd(NODE arg,LIST *rp)
532: {
533: DCP dc=0;
534:
535: if ( !ARG0(arg) ) {
536: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
537: } else
1.2 ! noro 538: modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),NEWDDD,&dc);
1.1 noro 539: dcptolist(dc,rp);
540: }
541:
542: void Pirred_check(NODE arg,Z *rp)
543: {
544: P p;
545: UM mp;
546: int r,mod;
547:
548: p = (P)ARG0(arg);
549: if ( !p ) {
550: *rp = 0; return;
551: }
552: mp = W_UMALLOC(UDEG(p));
1.2 ! noro 553: mod = ZTOS((Q)ARG1(arg));
1.1 noro 554: ptoum(mod,p,mp);
555: r = irred_check(mp,mod);
556: if ( r )
557: *rp = ONE;
558: else
559: *rp = 0;
560: }
561:
562: void Pnfctr_mod(NODE arg,Z *rp)
563: {
564: P p;
565: UM mp;
566: int r,mod;
567:
568: p = (P)ARG0(arg);
569: if ( !p ) {
570: *rp = 0; return;
571: }
572: mp = W_UMALLOC(UDEG(p));
1.2 ! noro 573: mod = ZTOS((Q)ARG1(arg));
1.1 noro 574: ptoum(mod,p,mp);
575: r = nfctr_mod(mp,mod);
1.2 ! noro 576: STOZ(r,*rp);
1.1 noro 577: }
578:
579: void Pddd_tab(NODE arg,VECT *rp)
580: {
581: P p;
582: UM mp,t,q,r1,w,w1;
583: UM *r,*s;
584: int dr,mod,n,i;
585: VECT result;
586: V v;
587:
1.2 ! noro 588: p = (P)ARG0(arg); mod = ZTOS((Q)ARG1(arg));
1.1 noro 589: v = VR(p);
590: n = UDEG(p); mp = W_UMALLOC(n);
591: ptoum(mod,p,mp);
592: r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
593: r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
594: t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
595: DEG(t) = mod; COEF(t)[mod] = 1;
596: q = W_UMALLOC(mod);
597: dr = divum(mod,t,mp,q);
598: DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
599: s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
600: w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
601: w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
602: for ( i = 1; i < n; i++ ) {
603: DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
604: mulum(mod,r1,w,w1);
605: dr = divum(mod,w1,mp,q); DEG(w1) = dr;
606: s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
607: }
608: for ( i = 2; i < n; i++ ) {
609: mult_mod_tab(r[i-1],mod,s,w,n);
610: r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
611: }
612: MKVECT(result,n);
613: for ( i = 0; i < n; i++ )
614: umtop(v,r[i],(P *)&BDY(result)[i]);
615: *rp = result;
616: }
617:
618: void reduce_sfdc(DCP sfdc,DCP *dcr)
619: {
620: P c,t,s,u,f;
621: DCP dc0,dc,tdc;
622: DCP *a;
623: int i,j,n;
624:
625: if ( !current_gfs_ext ) {
626: /* we simply apply sfptop() */
627: for ( dc0 = 0; sfdc; sfdc = NEXT(sfdc) ) {
628: NEXTDC(dc0,dc);
629: DEG(dc) = DEG(sfdc);
630: sfptop(COEF(sfdc),&COEF(dc));
631: }
632: NEXT(dc) = 0;
633: *dcr = dc0;
634: return;
635: }
636:
637: if ( NUM(COEF(sfdc)) ) {
638: sfptop(COEF(sfdc),&c);
639: sfdc = NEXT(sfdc);
640: } else
641: c = (P)ONE;
642:
643: for ( n = 0, tdc = sfdc; tdc; tdc = NEXT(tdc), n++ );
644: a = (DCP *)ALLOCA(n*sizeof(DCP));
645: for ( i = 0, tdc = sfdc; i < n; tdc = NEXT(tdc), i++ )
646: a[i] = tdc;
647:
648: dc0 = 0; NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = c;
649: for ( i = 0; i < n; i++ ) {
650: if ( !a[i] )
651: continue;
652: t = COEF(a[i]);
653: f = t;
654: while ( 1 ) {
655: sf_galois_action(t,(Q)ONE,&s);
656: for ( j = i; j < n; j++ )
657: if ( a[j] && !compp(CO,s,COEF(a[j])) )
658: break;
659: if ( j == n )
660: error("reduce_sfdc : cannot happen");
661: if ( j == i ) {
662: NEXTDC(dc0,dc); DEG(dc) = DEG(a[i]);
663: sfptop(f,&COEF(dc));
664: break;
665: } else {
666: mulp(CO,f,s,&u); f = u;
667: t = s;
668: a[j] = 0;
669: }
670: }
671: }
672: *dcr = dc0;
673: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / builtin