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