Annotation of OpenXM_contrib2/asir2000/engine/D.c, Revision 1.6
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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/D.c,v 1.5 2002/03/15 02:52:10 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: {
1.6 ! noro 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; ; ) {
1.1 noro 75: #if 0
1.6 ! noro 76: if ( !(++z % 10000) )
! 77: fprintf(stderr,"z=%d\n",z);
1.1 noro 78: #endif
1.6 ! noro 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;
1.1 noro 109: }
110:
111: void dtestsql(f,list,dc,dcp)
112: P f;
113: ML list;
114: struct oDUM *dc;
115: DCP *dcp;
116: {
1.6 ! noro 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;
1.1 noro 148: }
149:
150: void dtestsq(q,bound,f,fl,g)
151: int q,bound;
152: P f;
153: LUM fl;
154: P *g;
155: {
1.6 ! noro 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;
1.1 noro 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: {
1.6 ! noro 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: }
1.1 noro 196: }
197:
198: void dtestroot1(q,bound,f,fl,g)
199: int q,bound;
200: P f;
201: LUM fl;
202: P *g;
203: {
1.6 ! noro 204: P fq,fr,t;
1.1 noro 205:
1.6 ! noro 206: lumtop(VR(f),q,bound,fl,&t); udivpz(f,t,&fq,&fr);
! 207: *g = (fq && !fr) ? t : 0;
1.1 noro 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: {
1.6 ! noro 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;
1.1 noro 242: }
243:
244: int degtest(k,win,list,hint)
245: int k;
246: int *win;
247: ML list;
248: int hint;
249: {
1.6 ! noro 250: register int i,d;
! 251: LUM *c;
1.1 noro 252:
1.6 ! noro 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);
1.1 noro 258: }
259:
260: int ctest(g,list,k,in)
261: P g;
262: ML list;
263: int k;
264: int *in;
265: {
1.6 ! noro 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;
1.1 noro 291: }
292:
293: /*
294: int ncombi(n0,n,k,c)
295: int n0,n,k,*c;
296: {
1.6 ! noro 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;
1.1 noro 311: }
312: */
313:
314: int ncombi(n0,n,k,c)
315: int n0,n,k,*c;
316: {
1.6 ! noro 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;
1.1 noro 328: }
329:
330: void nthrootn(number,n,root)
331: int n;
332: N number,*root;
333: {
1.6 ! noro 334: N s,t,u,pn,base,n1,n2,q,r,gcd,num;
! 335: int sgn,index,p,i,tmp,tp,mlr,num0;
1.1 noro 336:
1.6 ! noro 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 = get_lprime(index);
! 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: }
1.1 noro 381: }
382:
383: void sqrtn(number,root)
384: N number,*root;
385: {
1.6 ! noro 386: N a,s,r,q;
! 387: int sgn;
1.1 noro 388:
1.6 ! noro 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: }
1.1 noro 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: {
1.6 ! noro 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: }
1.1 noro 436: }
1.6 ! noro 437:
1.1 noro 438: void padictoq(mod,bound,p,qp)
439: int mod,bound,*p;
440: Q *qp;
441: {
1.6 ! noro 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);
! 471: }
1.5 noro 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: {
1.6 ! noro 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: }
1.5 noro 504: }
1.6 ! noro 505:
1.5 noro 506: void padictoq_unsigned(mod,bound,p,qp)
507: int mod,bound,*p;
508: Q *qp;
509: {
1.6 ! noro 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);
! 527: }
1.1 noro 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: {
1.6 ! noro 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);
1.1 noro 557: }
558:
559: void ucsump(f,s)
560: P f;
561: Q *s;
562: {
1.6 ! noro 563: Q t,u;
! 564: DCP dc;
1.1 noro 565:
1.6 ! noro 566: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
! 567: addq((Q)COEF(dc),t,&u); t = u;
! 568: }
! 569: *s = t;
1.1 noro 570: }
571:
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 / asir2000 / engine