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