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