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