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