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