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