Annotation of OpenXM_contrib2/asir2018/engine/D.c, Revision 1.2
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.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/D.c,v 1.1 2018/09/19 05:45:06 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;
299: int sgn,index,p,i,tmp,tp,mlr,num0;
300:
301: for ( i = 0; !(n % 2); n /= 2, i++ );
1.2 ! noro 302: STOZ(n,z);
1.1 noro 303: for ( index = 0, num = number; ; index++ ) {
304: if ( n == 1 )
305: goto TAIL;
306: p = get_lprime(index);
307: if ( !(num0 = remqi((Q)num,p)) )
308: continue;
1.2 ! noro 309: STOZ(n,n1); STOZ(p-1,n2); gcdz(n1,n2,&gcd);
1.1 noro 310: if ( !UNIQ(gcd) )
311: continue;
1.2 ! noro 312: tp = pwrm(p,num0,invm(n,p-1)); STOZ(tp,s);
1.1 noro 313: mlr = invm(dmb(p,n,pwrm(p,tp,n-1),&tmp),p);
1.2 ! noro 314: STOZ(p,base); STOZ(p,pn);
1.1 noro 315: while ( 1 ) {
316: pwrz(s,z,&t); subz(num,t,&u);
317: if ( !u ) {
318: num = s;
319: break;
320: }
321: if ( sgnz(u) < 0 ) {
322: *root = 0;
323: return;
324: }
325: divqrz(u,base,&q,&r);
326: if ( r ) {
327: *root = 0;
328: return;
329: }
1.2 ! noro 330: STOZ(dmb(p,mlr,remqi((Q)q,p),&tmp),t);
1.1 noro 331: mulz(t,base,&u); addz(u,s,&t); s = t;
332: mulz(base,pn,&t); base = t;
333: }
334: TAIL :
335: for ( ; i; i-- ) {
336: sqrtz(num,&t);
337: if ( !t ) {
338: *root = 0;
339: return;
340: }
341: num = t;
342: }
343: *root = num;
344: return;
345: }
346: }
347:
348: void sqrtz(Z number,Z *root)
349: {
350: Z a,s,r,q,sa,two;
351: int sgn;
352:
353: for ( a = ONE; ; ) {
354: divqrz(number,a,&q,&r);
355: subz(q,a,&s);
356: if ( !s ) {
357: *root = !r ? a : 0;
358: return;
359: } else {
360: if ( sgnz(s) < 0 ) {
361: chsgnz(s,&sa); sgn = -1;
362: } else {
363: sa = s; sgn = 1;
364: }
365: if ( UNIZ(sa) ) {
366: *root = 0;
367: return;
368: } else {
1.2 ! noro 369: STOZ(2,two);
1.1 noro 370: divqrz(sa,two,&q,&r);
371: if ( sgn > 0 )
372: addz(a,q,&r);
373: else
374: subz(a,q,&r);
375: a = r;
376: }
377: }
378: }
379: }
380:
381: void lumtop(V v,int mod,int bound,LUM f,P *g)
382: {
383: DCP dc,dc0;
384: int **l;
385: int i;
386: Z q;
387:
388: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
389: padictoq(mod,bound,l[i],&q);
390: if ( q ) {
391: NEXTDC(dc0,dc);
392: if ( i )
1.2 ! noro 393: STOZ(i,DEG(dc));
1.1 noro 394: else
395: DEG(dc) = 0;
396: COEF(dc) = (P)q;
397: }
398: }
399: if ( !dc0 )
400: *g = 0;
401: else {
402: NEXT(dc) = 0; MKP(v,dc0,*g);
403: }
404: }
405:
406: void padictoq(int mod,int bound,int *p,Z *qp)
407: {
408: int h,i,t;
409: int br,sgn;
410: int *c;
411: Z z;
412:
413: c = W_ALLOC(bound);
414: for ( i = 0; i < bound; i++ )
415: c[i] = p[i];
416: h = (mod%2?(mod-1)/2:mod/2); i = bound - 1;
417: while ( i >= 0 && c[i] == h ) i--;
418: if ( i == -1 || c[i] > h ) {
419: for (i = 0, br = 0; i < bound; i++ )
420: if ( ( t = -(c[i] + br) ) < 0 ) {
421: c[i] = t + mod; br = 1;
422: } else {
423: c[i] = 0; br = 0;
424: }
425: sgn = -1;
426: } else
427: sgn = 1;
428: for ( i = bound - 1; ( i >= 0 ) && ( c[i] == 0 ); i--);
429: if ( i == -1 )
430: *qp = 0;
431: else {
432: nadictoz(mod,i+1,(unsigned int *)c,&z);
433: if ( sgn < 0 ) chsgnz(z,qp);
434: else *qp = z;
435: }
436: }
437:
438: void padictoq_unsigned(int,int,int *,Z *);
439:
440: void lumtop_unsigned(V v,int mod,int bound,LUM f,P *g)
441: {
442: DCP dc,dc0;
443: int **l;
444: int i;
445: Z q;
446:
447: for ( dc0 = NULL, i = DEG(f), l = COEF(f); i >= 0; i-- ) {
448: padictoq_unsigned(mod,bound,l[i],&q);
449: if ( q ) {
450: NEXTDC(dc0,dc);
451: if ( i )
1.2 ! noro 452: STOZ(i,DEG(dc));
1.1 noro 453: else
454: DEG(dc) = 0;
455: COEF(dc) = (P)q;
456: }
457: }
458: if ( !dc0 )
459: *g = 0;
460: else {
461: NEXT(dc) = 0; MKP(v,dc0,*g);
462: }
463: }
464:
465: void padictoq_unsigned(int mod,int bound,int *p,Z *qp)
466: {
467: nadictoz(mod,bound,(unsigned int *)p,qp);
468: }
469:
470: void mullumarray(P f,ML list,int k,int *in,P *g)
471: {
472: int np,bound,q,n,i,u;
473: int *tmpp;
474: LUM lclum,wb0,wb1,tlum;
475: LUM *l;
476: int len;
477: unsigned int *lc;
478:
479: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
480: W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
481: W_LUMALLOC(0,bound,lclum);
482: len = ztonadic(q,(Z)COEF(DC(f)),&lc);
483: for ( i = 0, tmpp = COEF(lclum)[0], u = MIN(bound,len); i < u; i++ )
484: tmpp[i] = lc[i];
485: l = (LUM *)list->c;
486: mullum(q,bound,lclum,l[in[0]],wb0);
487: for ( i = 1; i < k; i++ ) {
488: mullum(q,bound,l[in[i]],wb0,wb1);
489: tlum = wb0; wb0 = wb1; wb1 = tlum;
490: }
491: lumtop(VR(f),q,bound,wb0,g);
492: }
493:
494: void ucsump(P f,Q *s)
495: {
496: Q t,u;
497: DCP dc;
498:
499: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
500: addq((Q)COEF(dc),t,&u); t = u;
501: }
502: *s = t;
503: }
504:
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