![[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
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.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/D.c,v 1.4 2001/06/07 04:54:39 noro Exp $
1.2 noro 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;
1.4 noro 341: p = get_lprime(index);
1.1 noro 342: if ( !(num0 = rem(num,p)) )
343: continue;
344: STON(n,n1); STON(p-1,n2); gcdn(n1,n2,&gcd);
345: if ( !UNIN(gcd) )
346: continue;
347: tp = pwrm(p,num0,invm(n,p-1)); STON(tp,s);
348: mlr = invm(dmb(p,n,pwrm(p,tp,n-1),&tmp),p);
349: STON(p,base); STON(p,pn);
350: while ( 1 ) {
351: pwrn(s,n,&t); sgn = subn(num,t,&u);
352: if ( !u ) {
353: num = s;
354: break;
355: }
356: if ( sgn < 0 ) {
357: *root = 0;
358: return;
359: }
360: divn(u,base,&q,&r);
361: if ( r ) {
362: *root = 0;
363: return;
364: }
365: STON(dmb(p,mlr,rem(q,p),&tmp),t);
366: kmuln(t,base,&u); addn(u,s,&t); s = t;
367: kmuln(base,pn,&t); base = t;
368: }
369: TAIL :
370: for ( ; i; i-- ) {
371: sqrtn(num,&t);
372: if ( !t ) {
373: *root = 0;
374: return;
375: }
376: num = t;
377: }
378: *root = num;
379: return;
380: }
381: }
382:
383: void sqrtn(number,root)
384: N number,*root;
385: {
386: N a,s,r,q;
387: int sgn;
388:
389: for ( a = ONEN; ; ) {
390: divn(number,a,&q,&r); sgn = subn(q,a,&s);
391: if ( !s ) {
392: *root = !r ? a : 0;
393: return;
394: } else if ( UNIN(s) ) {
395: *root = 0;
396: return;
397: } else {
398: divin(s,2,&q);
399: if ( sgn > 0 )
400: addn(a,q,&r);
401: else
402: subn(a,q,&r);
403: a = r;
404: }
405: }
406: }
407:
408: void lumtop(v,mod,bound,f,g)
409: V v;
410: int mod;
411: int bound;
412: LUM f;
413: P *g;
414: {
415: DCP dc,dc0;
416: int **l;
417: int i;
418: Q q;
419:
420: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
421: padictoq(mod,bound,l[i],&q);
422: if ( q ) {
423: NEXTDC(dc0,dc);
424: if ( i )
425: STOQ(i,DEG(dc));
426: else
427: DEG(dc) = 0;
428: COEF(dc) = (P)q;
429: }
430: }
431: if ( !dc0 )
432: *g = 0;
433: else {
434: NEXT(dc) = 0; MKP(v,dc0,*g);
435: }
436: }
437:
438: void padictoq(mod,bound,p,qp)
439: int mod,bound,*p;
440: Q *qp;
441: {
442: register int h,i,t;
443: int br,sgn;
444: unsigned int *ptr;
445: N n,tn;
446: int *c;
447:
448: c = W_ALLOC(bound);
449: for ( i = 0; i < bound; i++ )
450: c[i] = p[i];
451: h = (mod%2?(mod-1)/2:mod/2); i = bound - 1;
452: while ( i >= 0 && c[i] == h ) i--;
453: if ( i == -1 || c[i] > h ) {
454: for (i = 0, br = 0; i < bound; i++ )
455: if ( ( t = -(c[i] + br) ) < 0 ) {
456: c[i] = t + mod; br = 1;
457: } else {
458: c[i] = 0; br = 0;
459: }
460: sgn = -1;
461: } else
462: sgn = 1;
463: for ( i = bound - 1; ( i >= 0 ) && ( c[i] == 0 ); i--);
464: if ( i == -1 )
465: *qp = 0;
466: else {
467: n = NALLOC(i+1); PL(n) = i+1;
468: for ( i = 0, ptr = BD(n); i < PL(n); i++ )
469: ptr[i] = c[i];
470: bnton(mod,n,&tn); NTOQ(tn,sgn,*qp);
1.5 ! noro 471: }
! 472: }
! 473:
! 474: void padictoq_unsigned(int,int,int *,Q *);
! 475:
! 476: void lumtop_unsigned(v,mod,bound,f,g)
! 477: V v;
! 478: int mod;
! 479: int bound;
! 480: LUM f;
! 481: P *g;
! 482: {
! 483: DCP dc,dc0;
! 484: int **l;
! 485: int i;
! 486: Q q;
! 487:
! 488: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
! 489: padictoq_unsigned(mod,bound,l[i],&q);
! 490: if ( q ) {
! 491: NEXTDC(dc0,dc);
! 492: if ( i )
! 493: STOQ(i,DEG(dc));
! 494: else
! 495: DEG(dc) = 0;
! 496: COEF(dc) = (P)q;
! 497: }
! 498: }
! 499: if ( !dc0 )
! 500: *g = 0;
! 501: else {
! 502: NEXT(dc) = 0; MKP(v,dc0,*g);
! 503: }
! 504: }
! 505:
! 506: void padictoq_unsigned(mod,bound,p,qp)
! 507: int mod,bound,*p;
! 508: Q *qp;
! 509: {
! 510: register int h,i,t;
! 511: int br,sgn;
! 512: unsigned int *ptr;
! 513: N n,tn;
! 514: int *c;
! 515:
! 516: c = W_ALLOC(bound);
! 517: for ( i = 0; i < bound; i++ )
! 518: c[i] = p[i];
! 519: for ( i = bound - 1; ( i >= 0 ) && ( c[i] == 0 ); i--);
! 520: if ( i == -1 )
! 521: *qp = 0;
! 522: else {
! 523: n = NALLOC(i+1); PL(n) = i+1;
! 524: for ( i = 0, ptr = BD(n); i < PL(n); i++ )
! 525: ptr[i] = c[i];
! 526: bnton(mod,n,&tn); NTOQ(tn,1,*qp);
1.1 noro 527: }
528: }
529:
530: void mullumarray(f,list,k,in,g)
531: P f;
532: ML list;
533: int k;
534: int *in;
535: P *g;
536: {
537: int np,bound,q,n,i,u;
538: int *tmpp;
539: LUM lclum,wb0,wb1,tlum;
540: LUM *l;
541: N lc;
542:
543: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
544: W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
545: W_LUMALLOC(0,bound,lclum);
546: ntobn(q,NM((Q)COEF(DC(f))),&lc);
547: for ( i = 0, tmpp = COEF(lclum)[0], u = MIN(bound,PL(lc));
548: i < u; i++ )
549: tmpp[i] = BD(lc)[i];
550: l = (LUM *)list->c;
551: mullum(q,bound,lclum,l[in[0]],wb0);
552: for ( i = 1; i < k; i++ ) {
553: mullum(q,bound,l[in[i]],wb0,wb1);
554: tlum = wb0; wb0 = wb1; wb1 = tlum;
555: }
556: lumtop(VR(f),q,bound,wb0,g);
557: }
558:
559: void ucsump(f,s)
560: P f;
561: Q *s;
562: {
563: Q t,u;
564: DCP dc;
565:
566: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
567: addq((Q)COEF(dc),t,&u); t = u;
568: }
569: *s = t;
570: }
571: