Annotation of OpenXM_contrib2/asir2018/lib/nf, 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: extern IDIVV_REM,IDIVV_HIST;
! 51:
! 52: def nf_mod(B,G,MOD,DIR,PS)
! 53: {
! 54: setmod(MOD);
! 55: for ( D = 0, H = []; G; ) {
! 56: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 57: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
! 58: R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
! 59: CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
! 60: U = G-CR*R;
! 61: print(".",2);
! 62: if ( !U )
! 63: return D;
! 64: break;
! 65: }
! 66: }
! 67: if ( U )
! 68: G = U;
! 69: else {
! 70: D += dp_hm(G); G = dp_rest(G); print(!D,2);
! 71: }
! 72: }
! 73: return D;
! 74: }
! 75:
! 76: def nf_mod_ext(B,G,MOD,DIR,PS)
! 77: {
! 78: setmod(MOD);
! 79: for ( D = 0, H = []; G; ) {
! 80: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 81: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
! 82: R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
! 83: CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
! 84: U = G-CR*R;
! 85: H = cons([CR,strtov("t"+rtostr(car(L)))],H);
! 86: print([length(H)]);
! 87: if ( !U )
! 88: return [D,reverse(H)];
! 89: break;
! 90: }
! 91: }
! 92: if ( U )
! 93: G = U;
! 94: else {
! 95: D += dp_hm(G); G = dp_rest(G);
! 96: }
! 97: }
! 98: return [D,reverse(H)];
! 99: }
! 100:
! 101: def nf(B,G,M,PS)
! 102: {
! 103: for ( D = 0, H = []; G; ) {
! 104: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 105: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
! 106: GCD = igcd(dp_hc(G),dp_hc(R));
! 107: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
! 108: U = CG*G-dp_subd(G,R)*CR*R;
! 109: H = cons(car(L),H);
! 110: if ( !U )
! 111: return [D,M,reverse(H)];
! 112: D *= CG; M *= CG;
! 113: break;
! 114: }
! 115: }
! 116: if ( U )
! 117: G = U;
! 118: else {
! 119: D += dp_hm(G); G = dp_rest(G);
! 120: }
! 121: }
! 122: return [D,M,reverse(H)];
! 123: }
! 124: extern M_NM,M_DN;
! 125:
! 126: def nf_demand(B,G,DIR,PSH)
! 127: {
! 128: T1 = T2 = T3 = T4 = 0;
! 129: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
! 130: for ( D = 0, H = []; G; ) {
! 131: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 132: if ( dp_redble(G,PSH[car(L)]) > 0 ) {
! 133: R = bload(DIR+rtostr(car(L)));
! 134: H = cons(car(L),H);
! 135: print([length(H),dp_mag(dp_hm(G))]);
! 136: GCD = igcd(dp_hc(G),dp_hc(R));
! 137: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
! 138: T0 = time()[0]; A1 = CG*G; T1 += time()[0]-T0;
! 139: T0 = time()[0]; A2 = dp_subd(G,R)*CR*R; T2 += time()[0]-T0;
! 140: T0 = time()[0]; U = A1-A2; T3 += time()[0]-T0;
! 141: if ( !U )
! 142: return [D,reverse(H),T1,T2,T3,T4];
! 143: D *= CG;
! 144: break;
! 145: }
! 146: }
! 147: if ( U ) {
! 148: G = U;
! 149: if ( D ) {
! 150: if ( dp_mag(dp_hm(D)) > Mag ) {
! 151: T0 = time()[0];
! 152: print("ptozp2");
! 153: T = dp_ptozp2(D,G); D = T[0]; G = T[1];
! 154: T4 += time()[0]-T0;
! 155: Mag = idiv(dp_mag(dp_hm(D))*M_NM,M_DN);
! 156: }
! 157: } else {
! 158: if ( dp_mag(dp_hm(G)) > Mag ) {
! 159: T0 = time()[0];
! 160: print("ptozp");
! 161: G = dp_ptozp(G);
! 162: T4 += time()[0]-T0;
! 163: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
! 164: }
! 165: }
! 166: } else {
! 167: D += dp_hm(G); G = dp_rest(G);
! 168: }
! 169: }
! 170: return [D,reverse(H),T1,T2,T3,T4];
! 171: }
! 172:
! 173: def nf_demand_d(B,G,DIR,NDIR,PDIR,PSH,PROC)
! 174: {
! 175: INDEX = 0;
! 176: T1 = T2 = 0;
! 177: /* dp_gr_flags(["Dist",PROC]); */
! 178: if ( PROC != [] ) {
! 179: PROC = newvect(length(PROC),PROC);
! 180: NPROC = size(PROC)[0];
! 181: } else
! 182: NPROC = 0;
! 183: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
! 184: Kara = dp_set_kara()*27; /* XXX */
! 185: D = [0,0]; R = [1,G];
! 186: print("");
! 187: for ( H = []; R[1]; ) {
! 188: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 189: if ( dp_redble(R[1],PSH[car(L)]) > 0 ) {
! 190: Red = bload(DIR+rtostr(car(L)));
! 191: H = cons(car(L),H);
! 192: /* D0=dp_mag(D[0]*<<0>>); D1=dp_mag(dp_hm(D[1]));
! 193: R0=dp_mag(R[0]*<<0>>); R1=dp_mag(dp_hm(R[1]));
! 194: print([car(L),length(H),[D0,D1,dp_nt(D[1])],[R0,R1,dp_nt(R[1])]],2); */
! 195:
! 196: GCD = igcd(dp_hc(R[1]),dp_hc(Red));
! 197: CR = idiv(dp_hc(Red),GCD);
! 198: CRed = idiv(dp_hc(R[1]),GCD);
! 199:
! 200: T0 = time()[3];
! 201: if ( (PROC != []) && (dp_mag(dp_hm(Red)) > Kara) ) {
! 202: print("d",2);
! 203: rpc(PROC[0],"dp_imul_index",CRed,car(L));
! 204: U = CR*R[1] - dp_subd(R[1],Red)*rpcrecv(PROC[0]);
! 205: } else {
! 206: print("l",2);
! 207: U = CR*R[1] - dp_subd(R[1],Red)*CRed*Red;
! 208: }
! 209: TT = time()[3]-T0; T1 += TT; /* print(TT); */
! 210: if ( !U )
! 211: return [D,reverse(H),T1,T2];
! 212: break;
! 213: }
! 214: }
! 215: if ( U ) {
! 216: if ( (HMAG=dp_mag(dp_hm(U))) > Mag ) {
! 217: T0 = time()[3];
! 218: if ( HMAG > Kara ) {
! 219: print("D",2);
! 220: T = dp_ptozp_d(U,PROC,NPROC);
! 221: } else {
! 222: print("L",2);
! 223: T = dp_ptozp(U);
! 224: }
! 225: TT = time()[3]-T0; T2 += TT; /* print(TT); */
! 226: Cont = idiv(dp_hc(U),dp_hc(T));
! 227: D0 = kmul(CR,D[0]);
! 228: R0 = kmul(Cont,R[0]);
! 229: S = igcd(D0,R0);
! 230: D = [idiv(D0,S),D[1]];
! 231: R = [idiv(R0,S),T];
! 232: Mag = idiv(dp_mag(dp_hm(R[1]))*M_NM,M_DN);
! 233: } else {
! 234: D = [kmul(CR,D[0]),D[1]];
! 235: R = [R[0],U];
! 236: }
! 237: } else {
! 238: C = kmul(dp_hc(R[1]),R[0]);
! 239: T = igcd(D[0],C);
! 240: D = [T,idiv(D[0],T)*D[1]+idiv(C,T)*dp_ht(R[1])];
! 241: R = [R[0],dp_rest(R[1])];
! 242: }
! 243: }
! 244: return [D[1],reverse(H),T1,T2];
! 245: }
! 246:
! 247: extern PTOZP,DPCV,NFSDIR;
! 248:
! 249: def imulv(A,B)
! 250: {
! 251: return A*B;
! 252: }
! 253:
! 254: def dp_imul(P,N)
! 255: {
! 256: return N*P;
! 257: }
! 258:
! 259: def mul(A,B)
! 260: {
! 261: return A*B;
! 262: }
! 263:
! 264: def dp_imul_index(C,INDEX)
! 265: {
! 266: T = C*bload(NFSDIR+rtostr(INDEX));
! 267: return T;
! 268: }
! 269:
! 270: def reg_dp_dtov(P)
! 271: {
! 272: PTOZP = P;
! 273: DPCV = dp_dtov(PTOZP);
! 274: }
! 275:
! 276: def reg_dp_iqr(G)
! 277: {
! 278: N = size(DPCV)[0];
! 279: for ( R = [], I = 0; I < N; I++ ) {
! 280: DPCV[I] = T = irem(DPCV[I],G);
! 281: if ( T )
! 282: R = cons(T,R);
! 283: }
! 284: return R;
! 285: }
! 286:
! 287: def reg_dp_cont(P)
! 288: {
! 289: PTOZP = P;
! 290: C = dp_dtov(P);
! 291: return igcd(C);
! 292: }
! 293:
! 294: def reg_dp_idiv(GCD)
! 295: {
! 296: return dp_idiv(PTOZP,GCD);
! 297: }
! 298:
! 299: def dp_cont_d(G,PROC,NPROC)
! 300: {
! 301: C = dp_sep(G,NPROC);
! 302: N = size(C)[0];
! 303: for ( I = 0; I < N; I++ )
! 304: rpc(PROC[I],"reg_dp_cont",C[I]);
! 305: R = map(rpcrecv,PROC);
! 306: GCD = igcd(R);
! 307: return GCD;
! 308: }
! 309:
! 310: /*
! 311: def iqrv(C,D)
! 312: {
! 313: return map(iqr,C,D);
! 314: }
! 315: */
! 316:
! 317: #if 1
! 318: def dp_ptozp_d(G,PROC,NPROC)
! 319: {
! 320: C = dp_dtov(G); L = size(C)[0];
! 321: T0 = time()[3];
! 322: D0 = igcd_estimate(C);
! 323: TT = time()[3]-T0; TG1 += TT;
! 324: CS = sepvect(C,NPROC+1);
! 325: N = size(CS)[0]; N1 = N-1;
! 326: QR = newvect(N); Q = newvect(L); R = newvect(L);
! 327: T0 = time()[3];
! 328: for ( I = 0; I < N1; I++ )
! 329: rpc0(PROC[I],"iqrv",CS[I],D0);
! 330: QR[N1] = iqrv(CS[N1],D0);
! 331: for ( I = 0; I < N1; I++ )
! 332: QR[I] = rpcrecv(PROC[I]);
! 333: TT = time()[3]-T0; TD += TT;
! 334: for ( I = J = 0; I < N; I++ ) {
! 335: T = QR[I]; M = size(T)[0];
! 336: for ( K = 0; K < M; K++, J++ ) {
! 337: Q[J] = T[K][0]; R[J] = T[K][1];
! 338: }
! 339: }
! 340: T0 = time()[3];
! 341: D1 = igcd(R);
! 342: GCD = igcd(D0,D1);
! 343: A = idiv(D0,GCD);
! 344: TT = time()[3]-T0; TG2 += TT;
! 345: T0 = time()[3];
! 346: for ( I = 0; I < L; I++ )
! 347: Q[I] = kmul(A,Q[I])+idiv(R[I],GCD);
! 348: TT = time()[3]-T0; TM += TT;
! 349: print([TG1,TD,TG2,TM,dp_mag(D0*<<0>>),dp_mag(D1*<<0>>),dp_mag(GCD*<<0>>)],2);
! 350: return dp_vtod(Q,G);
! 351: }
! 352: #endif
! 353:
! 354: #if 0
! 355: def dp_ptozp_d(G,PROC,NPROC)
! 356: {
! 357: C = dp_sep(G,NPROC+1);
! 358: N = size(C)[0]; N1 = N-1;
! 359: R = newvect(N);
! 360: for ( I = 0; I < N1; I++ )
! 361: rpc(PROC[I],"reg_igcdv",dp_dtov(C[I]));
! 362: R[N1] = reg_igcdv(dp_dtov(C[N1]));
! 363: for ( I = 0; I < N1; I++ )
! 364: R[I] = rpcrecv(PROC[I]);
! 365: GCD = igcd(R);
! 366: for ( I = 0; I < N1; I++ )
! 367: rpc(PROC[I],"idivv_final",GCD);
! 368: S = dp_vtod(idivv_final(GCD),C[N1]);
! 369: for ( I = 0; I < N1; I++ )
! 370: S += dp_vtod(rpcrecv(PROC[I]),C[I]);
! 371: return S;
! 372: }
! 373: #endif
! 374:
! 375: #if 0
! 376: def dp_ptozp_d(G,PROC,NPROC)
! 377: {
! 378: C = dp_sep(G,NPROC+1);
! 379: N = size(C)[0]; N1 = N-1;
! 380: R = newvect(N);
! 381: for ( I = 0; I < N1; I++ )
! 382: rpc(PROC[I],"reg_dp_cont",C[I]);
! 383: R[N1] = igcd(dp_dtov(C[N1]));
! 384: for ( I = 0; I < N1; I++ )
! 385: R[I] = rpcrecv(PROC[I]);
! 386: GCD = igcd(R);
! 387: for ( I = 0; I < N1; I++ )
! 388: rpc(PROC[I],"reg_dp_idiv",GCD);
! 389: S = dp_idiv(C[N1],GCD);
! 390: for ( I = 0; I < N1; I++ )
! 391: S += rpcrecv(PROC[I]);
! 392: return S;
! 393: }
! 394: #endif
! 395:
! 396: #if 0
! 397: def dp_ptozp_d(G,PROC,NPROC)
! 398: {
! 399: C = dp_sep(G,NPROC);
! 400: N = size(C)[0]; T = newvect(N);
! 401: for ( I = 0; I < N; I++ )
! 402: T[I] = PROC[I];
! 403: PROC = T;
! 404: for ( I = 0; I < N; I++ )
! 405: rpc(PROC[I],"reg_dp_dtov",C[I]);
! 406: A = dp_dtov(G); A = isort(A); L = size(A)[0];
! 407: map(rpcrecv,PROC);
! 408: while ( 1 ) {
! 409: GCD = igcd(A[0],A[1]);
! 410: for ( I = 2; I < L; I++ ) {
! 411: GT = igcd(GCD,A[I]);
! 412: if ( GCD == GT )
! 413: break;
! 414: else
! 415: GCD = GT;
! 416: }
! 417: if ( I == L )
! 418: break;
! 419: else {
! 420: map(rpc,PROC,"reg_dp_iqr",GCD);
! 421: R = map(rpcrecv,PROC);
! 422: for ( J = 0, U = []; J < N; J++ )
! 423: U = append(R[J],U);
! 424: L = length(U);
! 425: if ( L == 0 )
! 426: break;
! 427: else
! 428: U = cons(GCD,U);
! 429: A = isort(newvect(L+1,U));
! 430: print(["length=",L,"minmag=",dp_mag(A[0]*<<0>>)]);
! 431: }
! 432: }
! 433: for ( I = 0; I < N; I++ )
! 434: rpc(PROC[I],"reg_dp_idiv",GCD);
! 435: R = map(rpcrecv,PROC);
! 436: for ( I = 0, S = 0; I < N; I++ )
! 437: S += R[I];
! 438: return S;
! 439: }
! 440: #endif
! 441:
! 442: def dp_ptozp2_d(D,G,PROC,NPROC)
! 443: {
! 444: T = dp_ptozp_d(D+G,PROC,NPROC);
! 445: for ( S = 0; D; D = dp_rest(D), T = dp_rest(T) )
! 446: S += dp_hm(T);
! 447: return [S,T];
! 448: }
! 449:
! 450: def genindex(N)
! 451: {
! 452: R = [];
! 453: for ( I = 0; I < N; I++ )
! 454: R = cons(I,R);
! 455: return reverse(R);
! 456: }
! 457:
! 458: def nftest_ext_mod(N1,N2,N,MOD,DIR,PSH)
! 459: {
! 460: S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
! 461: R = nf_mod_ext(genindex(N-1),S,MOD,DIR,PSH);
! 462: return R;
! 463: }
! 464:
! 465: def nftest_mod(N1,N2,N,MOD,DIR,PSH)
! 466: {
! 467: S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
! 468: R = nf_mod(genindex(N-1),S,MOD,DIR,PSH);
! 469: return dp_mdtod(R);
! 470: }
! 471:
! 472: def nftest(N1,N2,N,DIR,PSH)
! 473: {
! 474: S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
! 475: R = nf_demand(genindex(N-1),S,DIR,PSH);
! 476: return R;
! 477: }
! 478:
! 479: def nftest_d(N1,N2,N,DIR,NDIR,PDIR,PSH,PROC)
! 480: {
! 481: S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
! 482: R = nf_demand_d(genindex(N-1),S,DIR,NDIR,PDIR,PSH,PROC);
! 483: return R;
! 484: }
! 485:
! 486: def abs(A)
! 487: {
! 488: return A >= 0 ? A : -A;
! 489: }
! 490:
! 491: def sort(A)
! 492: {
! 493: N = size(A)[0];
! 494: for ( I = 0; I < N; I++ ) {
! 495: for ( M = I, J = I + 1; J < N; J++ )
! 496: if ( A[J] < A[M] )
! 497: M = J;
! 498: if ( I != M ) {
! 499: T = A[I]; A[I] = A[M]; A[M] = T;
! 500: }
! 501: }
! 502: return A;
! 503: }
! 504:
! 505: def igcdv(C)
! 506: {
! 507: C = sort(C); N = size(C)[0];
! 508: G = igcd(C[0],C[1]);
! 509: for ( I = 2; I < N; I++ ) {
! 510: T = time();
! 511: G = igcd(G,C[I]);
! 512: /* print([I,dp_mag(G*<<0>>),time()[0]-T[0]]); */
! 513: }
! 514: return G;
! 515: }
! 516:
! 517: def igcdv2(C)
! 518: {
! 519: C = sort(C); N = size(C)[0];
! 520: G = igcd(C[0],C[1]);
! 521: for ( I = 2; I < N; I++ ) {
! 522: T = time();
! 523: C[I] %= G;
! 524: GT = igcd(G,C[I]);
! 525: if ( G == GT ) {
! 526: for ( J = I+1, NZ=0; J < N; J++ ) {
! 527: C[J] %= G;
! 528: if ( C[J] )
! 529: NZ = 1;
! 530: }
! 531: if ( !NZ )
! 532: break;
! 533: } else
! 534: G = GT;
! 535: print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
! 536: }
! 537: return G;
! 538: }
! 539:
! 540: def igcdv_d(C,P,NP)
! 541: {
! 542: C = isort(C); N = size(C)[0];
! 543: G = igcd(C[0],C[1]);
! 544: for ( I = 2; I < N; I++ ) {
! 545: T = time();
! 546: C[I] %= G;
! 547: GT = igcd(G,C[I]);
! 548: if ( G == GT ) {
! 549: for ( J = I+1, NZ=0; J < N; J++ ) {
! 550: C[J] %= G;
! 551: if ( C[J] )
! 552: NZ = 1;
! 553: }
! 554: if ( !NZ )
! 555: break;
! 556: } else
! 557: G = GT;
! 558: print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
! 559: }
! 560: return G;
! 561: }
! 562:
! 563: def dup(A)
! 564: {
! 565: N = size(A)[0]; V = newvect(N);
! 566: for ( I = 0; I < N; I++ )
! 567: V[I] = A[I];
! 568: return V;
! 569: }
! 570:
! 571: def idivv_init(A)
! 572: {
! 573: IDIVV_REM = dup(A);
! 574: }
! 575:
! 576: def igcd_cofactor(D)
! 577: {
! 578: T = map(iqr,IDIVV_REM,D);
! 579: N = size(T)[0];
! 580: Q = newvect(N); R = newvect(N);
! 581: for ( I = 0; I < N; I++ ) {
! 582: Q[I] = T[I][0]; R[I] = T[I][1];
! 583: }
! 584: D1 = igcdv(dup(R));
! 585: if ( !D1 ) {
! 586: IDIVV_HIST = [[Q,D]];
! 587: return D;
! 588: } else {
! 589: Q1 = map(idiv,R,D1);
! 590: IDIVV_HIST = [[Q,D],[Q1,D1]];
! 591: return D1;
! 592: }
! 593: }
! 594:
! 595: def idivv_step(D)
! 596: {
! 597: T = map(iqr,IDIVV_REM,D);
! 598: N = size(T)[0];
! 599: Q = newvect(N); R = newvect(N);
! 600: for ( I = 0, NZ = 0; I < N; I++ ) {
! 601: Q[I] = T[I][0]; R[I] = T[I][1];
! 602: if ( R[I] )
! 603: NZ = 1;
! 604: }
! 605: IDIVV_REM = R;
! 606: IDIVV_HIST = cons([Q,D],IDIVV_HIST);
! 607: return NZ;
! 608: }
! 609:
! 610: def igcdv3(C)
! 611: {
! 612: idivv_init(C);
! 613: C = isort(C); N = size(C)[0];
! 614: D = igcd(C[0],C[1]);
! 615: for ( I = 2; I < N; I++ ) {
! 616: T = time();
! 617: C[I] %= G;
! 618: GT = igcd(G,C[I]);
! 619: if ( G == GT ) {
! 620: NZ = idivv_step(G);
! 621: if ( !NZ )
! 622: break;
! 623: } else
! 624: G = GT;
! 625: }
! 626: CF = idivv_final(G);
! 627: return [G,CF];
! 628: }
! 629:
! 630: def igcdv4(C)
! 631: {
! 632: idivv_init(C);
! 633: C = isort(C); N = size(C)[0];
! 634: D = igcd_estimate(C);
! 635: G = igcd_cofactor(D);
! 636: G = igcd(G,D);
! 637: CF = idivv_final(G);
! 638: return [G,CF];
! 639: }
! 640:
! 641: def igcd_estimate(C)
! 642: {
! 643: C = isort(C); N = size(C)[0];
! 644: M = idiv(N,2);
! 645: for ( I = 0, S = 0; I < M; I++ )
! 646: S += C[I]>0?C[I]:-C[I];
! 647: for ( T = 0; I < N; I++ )
! 648: T += C[I]>0?C[I]:-C[I];
! 649: if ( !S && !T )
! 650: return igcdv(C);
! 651: else
! 652: return igcd(S,T);
! 653: }
! 654:
! 655: def reg_igcdv(C)
! 656: {
! 657: D = igcd_estimate(C);
! 658: T = map(iqr,C,D); N = size(T)[0];
! 659: Q = newvect(N); R = newvect(N);
! 660: for ( I = 0; I < N; I++ ) {
! 661: Q[I] = T[I][0]; R[I] = T[I][1];
! 662: }
! 663: D1 = igcdv(dup(R));
! 664: if ( !D1 ) {
! 665: IDIVV_HIST = [[Q,D]];
! 666: return D;
! 667: } else {
! 668: Q1 = map(idiv,R,D1);
! 669: IDIVV_HIST = [[Q,D],[Q1,D1]];
! 670: return igcd(D,D1);
! 671: }
! 672: }
! 673:
! 674: def idivv_final(D)
! 675: {
! 676: for ( T = IDIVV_HIST, S = 0; T != []; T = cdr(T) ) {
! 677: H = car(T); S += map(kmul,H[0],idiv(H[1],D));
! 678: }
! 679: return S;
! 680: }
! 681:
! 682: def dp_vtod(C,P)
! 683: {
! 684: for ( R = 0, I = 0; P; P = dp_rest(P), I++ )
! 685: R += C[I]*dp_ht(P);
! 686: return R;
! 687: }
! 688:
! 689: def dp_nt(P)
! 690: {
! 691: for ( I = 0; P; P = dp_rest(P), I++ );
! 692: return I;
! 693: }
! 694: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>