Annotation of OpenXM_contrib2/asir2018/engine/D.c, Revision 1.1
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: *
! 48: * $OpenXM$
! 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++ ) {
! 129: STOQ(dc[j].n,tq); pwrp(CO,true[j],tq,&s); udivpz(t,s,&fq,&fr);
! 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++ ) {
! 138: NEXTDC(dcr0,dcr); STOQ(dc[j].n,DEG(dcr)); COEF(dcr) = true[j];
! 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: }
! 173: STOQ(dc[0].n,q); pwrp(CO,t,q,&s); subp(CO,s,f,&u);
! 174: if ( u )
! 175: *dcp = 0;
! 176: else {
! 177: NEWDC(dcr); STOQ(dc[0].n,DEG(dcr));
! 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++ );
! 302: STOQ(n,z);
! 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;
! 309: STOQ(n,n1); STOQ(p-1,n2); gcdz(n1,n2,&gcd);
! 310: if ( !UNIQ(gcd) )
! 311: continue;
! 312: tp = pwrm(p,num0,invm(n,p-1)); STOQ(tp,s);
! 313: mlr = invm(dmb(p,n,pwrm(p,tp,n-1),&tmp),p);
! 314: STOQ(p,base); STOQ(p,pn);
! 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: }
! 330: STOQ(dmb(p,mlr,remqi((Q)q,p),&tmp),t);
! 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 {
! 369: STOQ(2,two);
! 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 )
! 393: STOQ(i,DEG(dc));
! 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 )
! 452: STOQ(i,DEG(dc));
! 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