Annotation of OpenXM_contrib2/asir2000/builtin/fctr.c, Revision 1.16
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.16 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/fctr.c,v 1.15 2002/10/23 07:54:57 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.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},
73: {"afctr",Pafctr,2},
74: {"agcd",Pagcd,3},
75: {"modsqfr",Pmodsqfr,2},
76: {"modfctr",Pmodfctr,2},
1.10 noro 77: {"sfsqfr",Psfsqfr,1},
1.15 noro 78: {"sffctr",Psffctr,1},
1.10 noro 79: {"sfufctr",Psfufctr,1},
80: {"sfbfctr",Psfbfctr,-4},
1.11 noro 81: {"sfmintdeg",Psfmintdeg,5},
1.13 noro 82: {"sfgcd",Psfgcd,2},
1.1 noro 83: #if 0
84: {"ddd",Pddd,2},
85: {"newddd",Pnewddd,2},
86: #endif
87: {"ddd_tab",Pddd_tab,2},
88: {"irred_check",Pirred_check,2},
89: {"nfctr_mod",Pnfctr_mod,2},
90: {0,0,0},
91: };
1.16 ! noro 92:
! 93: /* bivariate_hensel_special(f(x,y):monic in x,g0(x),h0(y),x,y,d) */
! 94:
! 95: void Pbivariate_hensel_special(arg,rp)
! 96: NODE arg;
! 97: LIST *rp;
! 98: {
! 99: DCP dc;
! 100: struct oVN vn[2];
! 101: P f,g0,h0,ak,bk,gk,hk;
! 102: V vx,vy;
! 103: VL nvl;
! 104: Q qk,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));
! 113: d = QTOS((Q)ARG5(arg));
! 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: }
1.1 noro 126:
127: void Pfctr(arg,rp)
128: NODE arg;
129: LIST *rp;
130: {
131: DCP dc;
132:
133: asir_assert(ARG0(arg),O_P,"fctr");
1.5 noro 134: if ( argc(arg) == 1 )
135: fctrp(CO,(P)ARG0(arg),&dc);
136: else {
137: asir_assert(ARG1(arg),O_P,"fctr");
138: fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc);
139: }
1.1 noro 140: dcptolist(dc,rp);
141: }
142:
143: void Pgcd(arg,rp)
144: NODE arg;
145: P *rp;
146: {
147: P p1,p2,g1,g2,g;
148: Num m;
149: int mod;
150:
151: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
152: asir_assert(p1,O_P,"gcd");
153: asir_assert(p2,O_P,"gcd");
154: if ( !p1 )
155: *rp = p2;
156: else if ( !p2 )
157: *rp = p1;
158: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
1.4 noro 159: gcdprsp(CO,p1,p2,rp);
1.1 noro 160: else if ( argc(arg) == 2 )
161: ezgcdp(CO,p1,p2,rp);
162: else {
163: m = (Num)ARG2(arg);
164: asir_assert(m,O_P,"gcd");
165: mod = QTOS((Q)m);
166: ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
167: gcdprsmp(CO,mod,g1,g2,&g);
168: mptop(g,rp);
169: }
170: }
171:
172: void Pgcdz(arg,rp)
173: NODE arg;
174: P *rp;
175: {
176: P p1,p2,t;
177: Q c1,c2;
178: N n;
179:
180: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
181: asir_assert(p1,O_P,"gcdz");
182: asir_assert(p2,O_P,"gcdz");
183: if ( !p1 )
184: *rp = p2;
185: else if ( !p2 )
186: *rp = p1;
187: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
188: error("gcdz : invalid argument");
189: else if ( NUM(p1) || NUM(p2) ) {
190: if ( NUM(p1) )
191: c1 = (Q)p1;
192: else
193: ptozp(p1,1,&c1,&t);
194: if ( NUM(p2) )
195: c2 = (Q)p2;
196: else
197: ptozp(p2,1,&c2,&t);
198: gcdn(NM(c1),NM(c2),&n); NTOQ(n,1,c1); *rp = (P)c1;
199: } else {
200: #if 0
201: w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
202: #endif
203: ezgcdpz(CO,p1,p2,rp);
204: }
205: }
206:
207: void Plcm(arg,rp)
208: NODE arg;
209: P *rp;
210: {
211: P t1,t2,p1,p2,g,q;
212: Q c;
213:
214: p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
215: asir_assert(p1,O_P,"lcm");
216: asir_assert(p2,O_P,"lcm");
217: if ( !p1 || !p2 )
218: *rp = 0;
219: else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
220: error("lcm : invalid argument");
221: else {
222: ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
223: ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
224: }
225: }
226:
227: void Psqfr(arg,rp)
228: NODE arg;
229: LIST *rp;
230: {
231: DCP dc;
232:
233: asir_assert(ARG0(arg),O_P,"sqfr");
234: sqfrp(CO,(P)ARG0(arg),&dc);
235: dcptolist(dc,rp);
236: }
237:
238: void Pufctrhint(arg,rp)
239: NODE arg;
240: LIST *rp;
241: {
242: DCP dc;
243:
244: asir_assert(ARG0(arg),O_P,"ufctrhint");
245: asir_assert(ARG1(arg),O_N,"ufctrhint");
246: ufctr((P)ARG0(arg),QTOS((Q)ARG1(arg)),&dc);
247: dcptolist(dc,rp);
248: }
249:
250: #if 0
251: Pmgcd(arg,rp)
252: NODE arg;
253: Obj *rp;
254: {
255: NODE node,tn;
256: int i,m;
257: P *l;
258:
259: node = BDY((LIST)ARG0(arg));
260: for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
261: m = i; l = (P *)ALLOCA(m*sizeof(P));
262: for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
263: l[i] = (P)BDY(tn);
264: nezgcdnpz(CO,l,m,rp);
265: }
266: #endif
267:
268: void Pcont(arg,rp)
269: NODE arg;
270: P *rp;
271: {
272: DCP dc;
273: int m;
274: P p,p1;
275: P *l;
276: V v;
277:
278: asir_assert(ARG0(arg),O_P,"cont");
279: p = (P)ARG0(arg);
280: if ( NUM(p) )
281: *rp = p;
282: else {
283: if ( argc(arg) == 2 ) {
284: v = VR((P)ARG1(arg));
285: change_mvar(CO,p,v,&p1);
286: if ( VR(p1) != v ) {
287: *rp = p1; return;
288: } else
289: p = p1;
290: }
291: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
292: l = (P *)ALLOCA(m*sizeof(P));
293: for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
294: l[m] = COEF(dc);
295: nezgcdnpz(CO,l,m,rp);
296: }
297: }
298:
299: void Pptozp(arg,rp)
300: NODE arg;
301: P *rp;
302: {
303: Q t;
304:
305: asir_assert(ARG0(arg),O_P,"ptozp");
306: ptozp((P)ARG0(arg),1,&t,rp);
307: }
308:
309: void Pafctr(arg,rp)
310: NODE arg;
311: LIST *rp;
312: {
313: DCP dc;
314:
315: asir_assert(ARG0(arg),O_P,"afctr");
316: asir_assert(ARG1(arg),O_P,"afctr");
317: afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
318: dcptolist(dc,rp);
319: }
320:
321: void Pagcd(arg,rp)
322: NODE arg;
323: P *rp;
324: {
325: asir_assert(ARG0(arg),O_P,"agcd");
326: asir_assert(ARG1(arg),O_P,"agcd");
327: asir_assert(ARG2(arg),O_P,"agcd");
328: gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
329: }
330:
331: #if 1
332: #define Mulum mulum
333: #define Divum divum
334: #define Mulsum mulsum
335: #define Gcdum gcdum
336: #endif
337:
338: void Mulum(), Mulsum(), Gcdum();
339: int Divum();
340:
341: #define FCTR 0 /* berlekamp */
342: #define SQFR 1
343: #define DDD 2 /* Cantor-Zassenhauss */
344: #define NEWDDD 3 /* berlekamp + root-finding by Cantor-Zassenhauss */
345:
346: UM *resberle();
347:
348: void Pmodfctr(arg,rp)
349: NODE arg;
350: LIST *rp;
351: {
352: DCP dc;
353: int mod;
354:
355: mod = QTOS((Q)ARG1(arg));
356: if ( mod < 0 )
357: error("modfctr : invalid modulus");
358: modfctrp(ARG0(arg),mod,NEWDDD,&dc);
1.6 noro 359: if ( !dc ) {
360: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
361: }
362: dcptolist(dc,rp);
1.13 noro 363: }
364:
365: void Psfgcd(arg,rp)
366: NODE arg;
367: LIST *rp;
368: {
369: P ps[2];
370:
371: ps[0] = (P)ARG0(arg);
372: ps[1] = (P)ARG1(arg);
373: gcdsf(CO,ps,2,rp);
1.6 noro 374: }
375:
1.15 noro 376: void Psffctr(arg,rp)
377: NODE arg;
378: LIST *rp;
379: {
380: DCP dc;
381:
382: mfctrsf(CO,ARG0(arg),&dc);
383: dcptolist(dc,rp);
384: }
385:
1.10 noro 386: void Psfsqfr(arg,rp)
387: NODE arg;
388: LIST *rp;
389: {
390: DCP dc;
391:
1.14 noro 392: sqfrsf(CO,ARG0(arg),&dc);
1.10 noro 393: dcptolist(dc,rp);
394: }
395:
396: void Psfufctr(arg,rp)
1.6 noro 397: NODE arg;
398: LIST *rp;
399: {
400: DCP dc;
401:
1.15 noro 402: ufctrsf(ARG0(arg),&dc);
1.1 noro 403: dcptolist(dc,rp);
1.7 noro 404: }
405:
406: void Psfbfctr(arg,rp)
407: NODE arg;
408: LIST *rp;
409: {
410: V x,y;
1.8 noro 411: DCP dc,dct;
412: P t;
413: struct oVL vl1,vl2;
414: VL vl;
1.10 noro 415: int degbound;
1.7 noro 416:
417: x = VR((P)ARG1(arg));
418: y = VR((P)ARG2(arg));
1.8 noro 419: vl1.v = x; vl1.next = &vl2;
420: vl2.v = y; vl2.next = 0;
421: vl = &vl1;
1.10 noro 422: if ( argc(arg) == 4 )
423: degbound = QTOS((Q)ARG3(arg));
424: else
425: degbound = -1;
1.8 noro 426:
1.10 noro 427: sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
1.8 noro 428: for ( dct = dc; dct; dct = NEXT(dct) ) {
429: reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
1.7 noro 430: }
1.8 noro 431: dcptolist(dc,rp);
1.1 noro 432: }
433:
1.11 noro 434: void Psfmintdeg(arg,rp)
435: NODE arg;
436: P *rp;
437: {
438: V x,y;
439: P r;
440: struct oVL vl1,vl2;
441: VL vl;
442: int dy,c;
443:
444: x = VR((P)ARG1(arg));
445: y = VR((P)ARG2(arg));
446: vl1.v = x; vl1.next = &vl2;
447: vl2.v = y; vl2.next = 0;
448: vl = &vl1;
449: dy = QTOS((Q)ARG3(arg));
450: c = QTOS((Q)ARG4(arg));
451: sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
452: reorderp(CO,vl,r,rp);
453: }
454:
1.1 noro 455: void Pmodsqfr(arg,rp)
456: NODE arg;
457: LIST *rp;
458: {
459: DCP dc;
460:
1.9 noro 461: if ( !ARG0(arg) ) {
1.1 noro 462: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 463: } else
464: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc);
1.1 noro 465: dcptolist(dc,rp);
466: }
467:
468: void Pddd(arg,rp)
469: NODE arg;
470: LIST *rp;
471: {
472: DCP dc;
473:
1.9 noro 474: if ( !ARG0(arg) ) {
1.1 noro 475: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 476: } else
477: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc);
1.1 noro 478: dcptolist(dc,rp);
479: }
480:
481: void Pnewddd(arg,rp)
482: NODE arg;
483: LIST *rp;
484: {
1.9 noro 485: DCP dc=0;
1.1 noro 486:
1.9 noro 487: if ( !ARG0(arg) ) {
1.1 noro 488: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 489: } else
490: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc);
1.1 noro 491: dcptolist(dc,rp);
492: }
493:
494: void Pirred_check(arg,rp)
495: NODE arg;
496: Q *rp;
497: {
498: P p;
499: UM mp;
500: int r,mod;
501:
502: p = (P)ARG0(arg);
503: if ( !p ) {
504: *rp = 0; return;
505: }
506: mp = W_UMALLOC(UDEG(p));
507: mod = QTOS((Q)ARG1(arg));
508: ptoum(mod,p,mp);
509: r = irred_check(mp,mod);
510: if ( r )
511: *rp = ONE;
512: else
513: *rp = 0;
514: }
515:
516: void Pnfctr_mod(arg,rp)
517: NODE arg;
518: Q *rp;
519: {
520: P p;
521: UM mp;
522: int r,mod;
523:
524: p = (P)ARG0(arg);
525: if ( !p ) {
526: *rp = 0; return;
527: }
528: mp = W_UMALLOC(UDEG(p));
529: mod = QTOS((Q)ARG1(arg));
530: ptoum(mod,p,mp);
531: r = nfctr_mod(mp,mod);
532: STOQ(r,*rp);
533: }
534:
535: void Pddd_tab(arg,rp)
536: NODE arg;
537: VECT *rp;
538: {
539: P p;
540: UM mp,t,q,r1,w,w1;
541: UM *r,*s;
542: int dr,mod,n,i;
543: VECT result;
544: V v;
545:
546: p = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
547: v = VR(p);
548: n = UDEG(p); mp = W_UMALLOC(n);
549: ptoum(mod,p,mp);
550: r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
551: r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
552: t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
553: DEG(t) = mod; COEF(t)[mod] = 1;
554: q = W_UMALLOC(mod);
555: dr = divum(mod,t,mp,q);
556: DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
557: s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
558: w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
559: w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
560: for ( i = 1; i < n; i++ ) {
561: DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
562: mulum(mod,r1,w,w1);
563: dr = divum(mod,w1,mp,q); DEG(w1) = dr;
564: s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
565: }
566: for ( i = 2; i < n; i++ ) {
567: mult_mod_tab(r[i-1],mod,s,w,n);
568: r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
569: }
570: MKVECT(result,n);
571: for ( i = 0; i < n; i++ )
572: umtop(v,r[i],(P *)&BDY(result)[i]);
573: *rp = result;
574: }
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