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