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