Annotation of OpenXM_contrib2/asir2000/builtin/fctr.c, Revision 1.13
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.13 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/fctr.c,v 1.12 2002/01/13 07:11:46 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.13 ! noro 57: void Psfsqfr(),Psfbfctr(),Psfufctr(),Psfmintdeg(),Psfgcd();
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:
1.1 noro 62: struct ftab fctr_tab[] = {
1.5 noro 63: {"fctr",Pfctr,-2},
1.1 noro 64: {"gcd",Pgcd,-3},
65: {"gcdz",Pgcdz,2},
66: {"lcm",Plcm,2},
67: {"sqfr",Psqfr,1},
68: {"ufctrhint",Pufctrhint,2},
69: {"ptozp",Pptozp,1},
70: {"cont",Pcont,-2},
71: {"afctr",Pafctr,2},
72: {"agcd",Pagcd,3},
73: {"modsqfr",Pmodsqfr,2},
74: {"modfctr",Pmodfctr,2},
1.10 noro 75: {"sfsqfr",Psfsqfr,1},
76: {"sfufctr",Psfufctr,1},
77: {"sfbfctr",Psfbfctr,-4},
1.11 noro 78: {"sfmintdeg",Psfmintdeg,5},
1.13 ! noro 79: {"sfgcd",Psfgcd,2},
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);
1.13 ! noro 326: }
! 327:
! 328: void Psfgcd(arg,rp)
! 329: NODE arg;
! 330: LIST *rp;
! 331: {
! 332: P ps[2];
! 333:
! 334: ps[0] = (P)ARG0(arg);
! 335: ps[1] = (P)ARG1(arg);
! 336: gcdsf(CO,ps,2,rp);
1.6 noro 337: }
338:
1.10 noro 339: void Psfsqfr(arg,rp)
340: NODE arg;
341: LIST *rp;
342: {
343: DCP dc;
344:
345: sfsqfr(ARG0(arg),&dc);
346: dcptolist(dc,rp);
347: }
348:
349: void Psfufctr(arg,rp)
1.6 noro 350: NODE arg;
351: LIST *rp;
352: {
353: DCP dc;
354:
355: fctrsf(ARG0(arg),&dc);
1.1 noro 356: dcptolist(dc,rp);
1.7 noro 357: }
358:
359: void Psfbfctr(arg,rp)
360: NODE arg;
361: LIST *rp;
362: {
363: V x,y;
1.8 noro 364: DCP dc,dct;
365: P t;
366: struct oVL vl1,vl2;
367: VL vl;
1.10 noro 368: int degbound;
1.7 noro 369:
370: x = VR((P)ARG1(arg));
371: y = VR((P)ARG2(arg));
1.8 noro 372: vl1.v = x; vl1.next = &vl2;
373: vl2.v = y; vl2.next = 0;
374: vl = &vl1;
1.10 noro 375: if ( argc(arg) == 4 )
376: degbound = QTOS((Q)ARG3(arg));
377: else
378: degbound = -1;
1.8 noro 379:
1.10 noro 380: sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
1.8 noro 381: for ( dct = dc; dct; dct = NEXT(dct) ) {
382: reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
1.7 noro 383: }
1.8 noro 384: dcptolist(dc,rp);
1.1 noro 385: }
386:
1.11 noro 387: void Psfmintdeg(arg,rp)
388: NODE arg;
389: P *rp;
390: {
391: V x,y;
392: P r;
393: struct oVL vl1,vl2;
394: VL vl;
395: int dy,c;
396:
397: x = VR((P)ARG1(arg));
398: y = VR((P)ARG2(arg));
399: vl1.v = x; vl1.next = &vl2;
400: vl2.v = y; vl2.next = 0;
401: vl = &vl1;
402: dy = QTOS((Q)ARG3(arg));
403: c = QTOS((Q)ARG4(arg));
404: sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
405: reorderp(CO,vl,r,rp);
406: }
407:
1.1 noro 408: void Pmodsqfr(arg,rp)
409: NODE arg;
410: LIST *rp;
411: {
412: DCP dc;
413:
1.9 noro 414: if ( !ARG0(arg) ) {
1.1 noro 415: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 416: } else
417: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc);
1.1 noro 418: dcptolist(dc,rp);
419: }
420:
421: void Pddd(arg,rp)
422: NODE arg;
423: LIST *rp;
424: {
425: DCP dc;
426:
1.9 noro 427: if ( !ARG0(arg) ) {
1.1 noro 428: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 429: } else
430: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc);
1.1 noro 431: dcptolist(dc,rp);
432: }
433:
434: void Pnewddd(arg,rp)
435: NODE arg;
436: LIST *rp;
437: {
1.9 noro 438: DCP dc=0;
1.1 noro 439:
1.9 noro 440: if ( !ARG0(arg) ) {
1.1 noro 441: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9 noro 442: } else
443: modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc);
1.1 noro 444: dcptolist(dc,rp);
445: }
446:
447: void Pirred_check(arg,rp)
448: NODE arg;
449: Q *rp;
450: {
451: P p;
452: UM mp;
453: int r,mod;
454:
455: p = (P)ARG0(arg);
456: if ( !p ) {
457: *rp = 0; return;
458: }
459: mp = W_UMALLOC(UDEG(p));
460: mod = QTOS((Q)ARG1(arg));
461: ptoum(mod,p,mp);
462: r = irred_check(mp,mod);
463: if ( r )
464: *rp = ONE;
465: else
466: *rp = 0;
467: }
468:
469: void Pnfctr_mod(arg,rp)
470: NODE arg;
471: Q *rp;
472: {
473: P p;
474: UM mp;
475: int r,mod;
476:
477: p = (P)ARG0(arg);
478: if ( !p ) {
479: *rp = 0; return;
480: }
481: mp = W_UMALLOC(UDEG(p));
482: mod = QTOS((Q)ARG1(arg));
483: ptoum(mod,p,mp);
484: r = nfctr_mod(mp,mod);
485: STOQ(r,*rp);
486: }
487:
488: void Pddd_tab(arg,rp)
489: NODE arg;
490: VECT *rp;
491: {
492: P p;
493: UM mp,t,q,r1,w,w1;
494: UM *r,*s;
495: int dr,mod,n,i;
496: VECT result;
497: V v;
498:
499: p = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
500: v = VR(p);
501: n = UDEG(p); mp = W_UMALLOC(n);
502: ptoum(mod,p,mp);
503: r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
504: r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
505: t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
506: DEG(t) = mod; COEF(t)[mod] = 1;
507: q = W_UMALLOC(mod);
508: dr = divum(mod,t,mp,q);
509: DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
510: s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
511: w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
512: w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
513: for ( i = 1; i < n; i++ ) {
514: DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
515: mulum(mod,r1,w,w1);
516: dr = divum(mod,w1,mp,q); DEG(w1) = dr;
517: s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
518: }
519: for ( i = 2; i < n; i++ ) {
520: mult_mod_tab(r[i-1],mod,s,w,n);
521: r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
522: }
523: MKVECT(result,n);
524: for ( i = 0; i < n; i++ )
525: umtop(v,r[i],(P *)&BDY(result)[i]);
526: *rp = result;
527: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>