![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to D.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to D.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
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/engine/D.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $
! 49: */
1.1 noro 50: #include "ca.h"
51:
52: void dtest(f,list,hint,dcp)
53: P f;
54: ML list;
55: int hint;
56: DCP *dcp;
57: {
58: int n,np,bound,q;
59: int i,j,k;
60: int *win;
61: P g,factor,cofactor;
62: Q csum,csumt;
63: DCP dcf,dcf0;
64: LUM *c;
65: ML wlist;
66: int z;
67:
68: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
69: win = W_ALLOC(np+1);
70: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
71: wlist = W_MLALLOC(np); wlist->n = list->n;
72: wlist->mod = list->mod; wlist->bound = list->bound;
73: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
74: for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; ) {
75: #if 0
76: if ( !(++z % 10000) )
77: fprintf(stderr,"z=%d\n",z);
78: #endif
79: if ( degtest(k,win,wlist,hint) &&
80: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
81: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
82: g = cofactor;
83: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
84: for ( i = 0; i < k - 1; i++ )
85: for ( j = win[i] + 1; j < win[i + 1]; j++ )
86: c[j-i-1] = c[j];
87: for ( j = win[k-1] + 1; j <= np; j++ )
88: c[j-k] = c[j];
89: if ( ( np -= k ) < k )
90: break;
91: if ( np - win[0] + 1 < k )
92: if ( ++k > np )
93: break;
94: else
95: for ( i = 0; i < k; i++ )
96: win[i] = i + 1;
97: else
98: for ( i = 1; i < k; i++ )
99: win[i] = win[0] + i;
100: } else if ( !ncombi(1,np,k,win) )
101: if ( k == np )
102: break;
103: else
104: for ( i = 0, ++k; i < k; i++ )
105: win[i] = i + 1;
106: }
107: NEXTDC(dcf0,dcf); COEF(dcf) = g;
108: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
109: }
110:
111: void dtestsql(f,list,dc,dcp)
112: P f;
113: ML list;
114: struct oDUM *dc;
115: DCP *dcp;
116: {
117: int j,n,m,b;
118: P t,s,fq,fr;
119: P *true;
120: Q tq;
121: LUM *c;
122: DCP dcr,dcr0;
123:
124: n = list->n; m = list->mod; b = list->bound; c = (LUM *)list->c;
125: true = (P*)ALLOCA(n*sizeof(P));
126: for ( j = 0; j < n; j++ ) {
127: dtestsq(m,b,f,c[j],&t);
128: if ( t )
129: true[j] = t;
130: else {
131: *dcp = 0;
132: return;
133: }
134: }
135: for ( t = f, j = 0; j < n; j++ ) {
136: STOQ(dc[j].n,tq); pwrp(CO,true[j],tq,&s); udivpz(t,s,&fq,&fr);
137: if ( fq && !fr )
138: t = fq;
139: else {
140: *dcp = 0;
141: return;
142: }
143: }
144: for ( j = 0, dcr = dcr0 = 0; j < n; j++ ) {
145: NEXTDC(dcr0,dcr); STOQ(dc[j].n,DEG(dcr)); COEF(dcr) = true[j];
146: }
147: NEXT(dcr) = 0; *dcp = dcr0;
148: }
149:
150: void dtestsq(q,bound,f,fl,g)
151: int q,bound;
152: P f;
153: LUM fl;
154: P *g;
155: {
156: P lcf,t,fq,fr,s;
157: struct oML list;
158: int in = 0;
159:
160: list.n = 1;
161: list.mod = q;
162: list.bound = bound;
163: list.c[0] = (pointer)fl;
164:
165: mullumarray(f,&list,1,&in,&t); mulp(CO,f,COEF(DC(f)),&lcf);
166: udivpz(lcf,t,&fq,&fr);
167: if( fq && !fr )
168: ptozp(t,1,(Q *)&s,g);
169: else
170: *g = 0;
171: }
172:
173: void dtestroot(m,b,f,fl,dc,dcp)
174: int m,b;
175: P f;
176: LUM fl;
177: struct oDUM *dc;
178: DCP *dcp;
179: {
180: P t,s,u;
181: DCP dcr;
182: Q q;
183:
184: dtestroot1(m,b,f,fl,&t);
185: if ( !t ) {
186: *dcp = 0;
187: return;
188: }
189: STOQ(dc[0].n,q); pwrp(CO,t,q,&s); subp(CO,s,f,&u);
190: if ( u )
191: *dcp = 0;
192: else {
193: NEWDC(dcr); STOQ(dc[0].n,DEG(dcr));
194: COEF(dcr) = t; NEXT(dcr) = 0; *dcp = dcr;
195: }
196: }
197:
198: void dtestroot1(q,bound,f,fl,g)
199: int q,bound;
200: P f;
201: LUM fl;
202: P *g;
203: {
204: P fq,fr,t;
205:
206: lumtop(VR(f),q,bound,fl,&t); udivpz(f,t,&fq,&fr);
207: *g = (fq && !fr) ? t : 0;
208: }
209:
210: int dtestmain(g,csum,list,k,in,fp,cofp)
211: P g;
212: Q csum;
213: ML list;
214: int k;
215: int *in;
216: P *fp,*cofp;
217: {
218: int mod;
219: P fmul,lcg;
220: Q csumg;
221: N nq,nr;
222: P fq,fr,t;
223:
224: if (!ctest(g,list,k,in))
225: return 0;
226: mod = list->mod;
227: mullumarray(g,list,k,in,&fmul); mulp(CO,g,COEF(DC(g)),&lcg);
228: if ( csum ) {
229: ucsump(fmul,&csumg);
230: if ( csumg ) {
231: divn(NM(csum),NM(csumg),&nq,&nr);
232: if ( nr )
233: return 0;
234: }
235: }
236: udivpz(lcg,fmul,&fq,&fr);
237: if ( fq && !fr ) {
238: ptozp(fq,1,(Q *)&t,cofp); ptozp(fmul,1,(Q *)&t,fp);
239: return 1;
240: } else
241: return 0;
242: }
243:
244: int degtest(k,win,list,hint)
245: int k;
246: int *win;
247: ML list;
248: int hint;
249: {
250: register int i,d;
251: LUM *c;
252:
253: if ( hint == 1 )
254: return 1;
255: for ( c = (LUM*)list->c, i = 0, d = 0; i < k; i++ )
256: d += DEG(c[win[i]]);
257: return !(d % hint);
258: }
259:
260: int ctest(g,list,k,in)
261: P g;
262: ML list;
263: int k;
264: int *in;
265: {
266: register int i;
267: int q,bound;
268: int *wm,*wm1,*tmpp;
269: DCP dc;
270: Q dvr;
271: N lcn,cstn,dndn,dmyn,rn;
272: LUM *l;
273:
274: for ( dc = DC(g); dc && DEG(dc); dc = NEXT(dc) );
275: if ( dc )
276: cstn = NM((Q)COEF(dc));
277: else
278: return 1;
279: q = list->mod; bound = list->bound;
280: ntobn(q,NM((Q)COEF(DC(g))),&lcn);;
281: W_CALLOC(bound+1,int,wm); W_CALLOC(bound+1,int,wm1);
282: for ( i = 0; i < PL(lcn); i++ )
283: wm[i] = BD(lcn)[i];
284: for ( i = 0, l = (LUM *)list->c; i < k; i++ ) {
285: mulpadic(q,bound,wm,COEF(l[in[i]])[0],wm1);
286: tmpp = wm; wm = wm1; wm1 = tmpp;
287: }
288: padictoq(q,bound,wm,&dvr);
289: kmuln(NM((Q)COEF(DC(g))),cstn,&dndn); divn(dndn,NM(dvr),&dmyn,&rn);
290: return rn ? 0 : 1;
291: }
292:
293: /*
294: int ncombi(n0,n,k,c)
295: int n0,n,k,*c;
296: {
297: register int i,tmp;
298:
299: if ( !k )
300: return 0;
301: if ( !ncombi(c[1],n,k-1,c+1) ) {
302: if ( c[0] + k > n )
303: return 0;
304: else {
305: for ( i = 0, tmp = c[0]; i < k; i++ )
306: c[i] = tmp + i + 1;
307: return 1;
308: }
309: } else
310: return 1;
311: }
312: */
313:
314: int ncombi(n0,n,k,c)
315: int n0,n,k,*c;
316: {
317: register int i,t;
318:
319: if ( !k )
320: return 0;
321: for ( i = k-1; i >= 0 && c[i] == n+i-(k-1); i-- );
322: if ( i < 0 )
323: return 0;
324: t = ++c[i++];
325: for ( t++ ; i < k; i++, t++ )
326: c[i] = t;
327: return 1;
328: }
329:
330: void nthrootn(number,n,root)
331: int n;
332: N number,*root;
333: {
334: N s,t,u,pn,base,n1,n2,q,r,gcd,num;
335: int sgn,index,p,i,tmp,tp,mlr,num0;
336:
337: for ( i = 0; !(n % 2); n /= 2, i++ );
338: for ( index = 0, num = number; ; index++ ) {
339: if ( n == 1 )
340: goto TAIL;
341: p = lprime[index];
342: if ( !p )
343: error("nthrootn : lprime[] exhausted.");
344: if ( !(num0 = rem(num,p)) )
345: continue;
346: STON(n,n1); STON(p-1,n2); gcdn(n1,n2,&gcd);
347: if ( !UNIN(gcd) )
348: continue;
349: tp = pwrm(p,num0,invm(n,p-1)); STON(tp,s);
350: mlr = invm(dmb(p,n,pwrm(p,tp,n-1),&tmp),p);
351: STON(p,base); STON(p,pn);
352: while ( 1 ) {
353: pwrn(s,n,&t); sgn = subn(num,t,&u);
354: if ( !u ) {
355: num = s;
356: break;
357: }
358: if ( sgn < 0 ) {
359: *root = 0;
360: return;
361: }
362: divn(u,base,&q,&r);
363: if ( r ) {
364: *root = 0;
365: return;
366: }
367: STON(dmb(p,mlr,rem(q,p),&tmp),t);
368: kmuln(t,base,&u); addn(u,s,&t); s = t;
369: kmuln(base,pn,&t); base = t;
370: }
371: TAIL :
372: for ( ; i; i-- ) {
373: sqrtn(num,&t);
374: if ( !t ) {
375: *root = 0;
376: return;
377: }
378: num = t;
379: }
380: *root = num;
381: return;
382: }
383: }
384:
385: void sqrtn(number,root)
386: N number,*root;
387: {
388: N a,s,r,q;
389: int sgn;
390:
391: for ( a = ONEN; ; ) {
392: divn(number,a,&q,&r); sgn = subn(q,a,&s);
393: if ( !s ) {
394: *root = !r ? a : 0;
395: return;
396: } else if ( UNIN(s) ) {
397: *root = 0;
398: return;
399: } else {
400: divin(s,2,&q);
401: if ( sgn > 0 )
402: addn(a,q,&r);
403: else
404: subn(a,q,&r);
405: a = r;
406: }
407: }
408: }
409:
410: void lumtop(v,mod,bound,f,g)
411: V v;
412: int mod;
413: int bound;
414: LUM f;
415: P *g;
416: {
417: DCP dc,dc0;
418: int **l;
419: int i;
420: Q q;
421:
422: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
423: padictoq(mod,bound,l[i],&q);
424: if ( q ) {
425: NEXTDC(dc0,dc);
426: if ( i )
427: STOQ(i,DEG(dc));
428: else
429: DEG(dc) = 0;
430: COEF(dc) = (P)q;
431: }
432: }
433: if ( !dc0 )
434: *g = 0;
435: else {
436: NEXT(dc) = 0; MKP(v,dc0,*g);
437: }
438: }
439:
440: void padictoq(mod,bound,p,qp)
441: int mod,bound,*p;
442: Q *qp;
443: {
444: register int h,i,t;
445: int br,sgn;
446: unsigned int *ptr;
447: N n,tn;
448: int *c;
449:
450: c = W_ALLOC(bound);
451: for ( i = 0; i < bound; i++ )
452: c[i] = p[i];
453: h = (mod%2?(mod-1)/2:mod/2); i = bound - 1;
454: while ( i >= 0 && c[i] == h ) i--;
455: if ( i == -1 || c[i] > h ) {
456: for (i = 0, br = 0; i < bound; i++ )
457: if ( ( t = -(c[i] + br) ) < 0 ) {
458: c[i] = t + mod; br = 1;
459: } else {
460: c[i] = 0; br = 0;
461: }
462: sgn = -1;
463: } else
464: sgn = 1;
465: for ( i = bound - 1; ( i >= 0 ) && ( c[i] == 0 ); i--);
466: if ( i == -1 )
467: *qp = 0;
468: else {
469: n = NALLOC(i+1); PL(n) = i+1;
470: for ( i = 0, ptr = BD(n); i < PL(n); i++ )
471: ptr[i] = c[i];
472: bnton(mod,n,&tn); NTOQ(tn,sgn,*qp);
473: }
474: }
475:
476: void mullumarray(f,list,k,in,g)
477: P f;
478: ML list;
479: int k;
480: int *in;
481: P *g;
482: {
483: int np,bound,q,n,i,u;
484: int *tmpp;
485: LUM lclum,wb0,wb1,tlum;
486: LUM *l;
487: N lc;
488:
489: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
490: W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
491: W_LUMALLOC(0,bound,lclum);
492: ntobn(q,NM((Q)COEF(DC(f))),&lc);
493: for ( i = 0, tmpp = COEF(lclum)[0], u = MIN(bound,PL(lc));
494: i < u; i++ )
495: tmpp[i] = BD(lc)[i];
496: l = (LUM *)list->c;
497: mullum(q,bound,lclum,l[in[0]],wb0);
498: for ( i = 1; i < k; i++ ) {
499: mullum(q,bound,l[in[i]],wb0,wb1);
500: tlum = wb0; wb0 = wb1; wb1 = tlum;
501: }
502: lumtop(VR(f),q,bound,wb0,g);
503: }
504:
505: void ucsump(f,s)
506: P f;
507: Q *s;
508: {
509: Q t,u;
510: DCP dc;
511:
512: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
513: addq((Q)COEF(dc),t,&u); t = u;
514: }
515: *s = t;
516: }
517: