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