![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gw.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / testing / noro|
Return to gw.rr CVS log | Up to [local] / OpenXM / src / asir-contrib / testing / noro |
1.1 ! noro 1: module noro_gw$
! 2:
! 3: localf set_order_mat, create_order_mat, inner_product, comp_order$
! 4: localf in_c12, comp_facet_preorder, compute_df, dp_boundary$
! 5: localf compute_compat_weight, compute_compat_weight_gmp, compute_last_w$
! 6: localf highest_w, mat_to_vlist, porder, nf_marked, appear, nf_marked2$
! 7: localf nf_marked_rec, comp_second, sort_by_order, generic_walk$
! 8: localf generic_walk_mod, gw, gw_mod, tolex_gw$
! 9:
! 10: static Cdd_Init, Cdd_Proc$
! 11: static Cddgmp_Init, Cddgmp_Proc$
! 12:
! 13: def set_order_mat(M,S,E,O) {
! 14: if ( O == 0 ) {
! 15: for ( J = S; J < E; J++ ) M[S][J] = 1;
! 16: for ( I = S+1; I < E; I++ ) M[I][E-(I-S)] = -1;
! 17: } else if ( O == 1 ) {
! 18: for ( J = S; J < E; J++ ) M[S][J] = 1;
! 19: for ( I = S+1; I < E; I++ ) M[I][I-1] = 1;
! 20: } else if ( O == 2 )
! 21: for ( I = S; I < E; I++ ) M[I][I] = 1;
! 22: }
! 23:
! 24: def create_order_mat(N,O)
! 25: {
! 26: M = matrix(N,N);
! 27: if ( O <= 2 )
! 28: set_order_mat(M,0,N,O);
! 29: else {
! 30: for ( T = O, S = 0; T != []; T = cdr(T) ) {
! 31: OT = car(T)[0]; ON = car(T)[1];
! 32: set_order_mat(M,S,S+ON,OT);
! 33: S += ON;
! 34: }
! 35: }
! 36: return M;
! 37: }
! 38:
! 39: def inner_product(V1,V2)
! 40: {
! 41: if ( V1 == 0 || V2 == 0 ) return 0;
! 42: N = size(V1)[0];
! 43: for ( S = 0, I = 0; I < N; I++ )
! 44: S += V1[I]*V2[I];
! 45: return S;
! 46: }
! 47:
! 48: def comp_order(V1,V2,W)
! 49: {
! 50: N = size(W)[0];
! 51: for ( I = 0; I < N; I++ ) {
! 52: T1 = inner_product(V1,W[I]);
! 53: T2 = inner_product(V2,W[I]);
! 54: if ( T1 < T2 ) return -1;
! 55: else if ( T1 > T2 ) return 1;
! 56: }
! 57: return 0;
! 58: }
! 59:
! 60: def in_c12(V,W1,W2)
! 61: {
! 62: T1 = comp_order(V,0,W1);
! 63: T2 = comp_order(V,0,W2);
! 64: if ( T1 > 0 && T2 < 0 ) return 1;
! 65: else return 0;
! 66: }
! 67:
! 68: /* U <= V according to facet preorder */
! 69:
! 70: def comp_facet_preorder(U,V,W1,W2)
! 71: {
! 72: /* U = 0 <=> U = -infty, U = 1 <=> U = infty */
! 73: if ( U == 0 ) return 1;
! 74: if ( U == 1 ) return 0;
! 75: N = size(V)[0];
! 76: for ( I = 0; I < N; I++ ) {
! 77: Left = inner_product(W2[I],U)*V;
! 78: Right = inner_product(W2[I],V)*U;
! 79: R = comp_order(Left,Right,W1);
! 80: if ( R < 0 ) return 1;
! 81: else if ( R > 0 ) return 0;
! 82: }
! 83: return 1;
! 84: }
! 85:
! 86: def compute_df(F,H)
! 87: {
! 88: H = dp_etov(H);
! 89: for ( R = [], T = F; T; T = dp_rest(T) )
! 90: R = cons(H-dp_etov(T),R);
! 91: return reverse(R);
! 92: }
! 93:
! 94: def dp_boundary(F)
! 95: {
! 96: for ( M = [], T = F; T; T = dp_rest(T) )
! 97: M = cons(dp_hm(T),M);
! 98: M = ltov(M); N = length(M);
! 99: for ( I = 0; I < N; I++ ) {
! 100: for ( MI = M[I], J = I+1; J < N; J++ ) {
! 101: if ( dp_redble(M[J],MI) ) break;
! 102: }
! 103: if ( J < N ) M[I] = 0;
! 104: }
! 105: for ( R = 0, I = 0; I < N; I++ )
! 106: R += M[I];
! 107: return R;
! 108: }
! 109:
! 110: def compute_compat_weight(F,H)
! 111: {
! 112: D = dp_compute_essential_df(F,H);
! 113: C = length(D[0]);
! 114: for ( I = 1; I < C; I++ ) {
! 115: V = vector(C);
! 116: V[I] = 1;
! 117: D = cons(vtol(V),D);
! 118: }
! 119: R = length(D);
! 120: if ( !Cdd_Init ) {
! 121: Cdd_Proc = ox_launch_nox(0,"ox_cdd");
! 122: Cdd_Init = 1;
! 123: }
! 124: ox_cmo_rpc(Cdd_Proc,"intpt",R,C,D);
! 125: Ret = ox_pop_cmo(Cdd_Proc);
! 126: V = vector(C-1);
! 127: for ( I = 1, D = 1; I < C; I++ ) {
! 128: V[I-1] = rint(Ret[I]/Ret[0]*100);
! 129: }
! 130: return V;
! 131: }
! 132:
! 133: def compute_compat_weight_gmp(F,H)
! 134: {
! 135: D = dp_compute_essential_df(F,H);
! 136: C = length(D[0]);
! 137: for ( I = 1; I < C; I++ ) {
! 138: V = vector(C);
! 139: V[I] = 1;
! 140: D = cons(vtol(V),D);
! 141: }
! 142: R = length(D);
! 143: if ( !Cddgmp_Init ) {
! 144: Cddgmp_Proc = ox_launch_nox(0,"ox_cddgmp");
! 145: Cddgmp_Init = 1;
! 146: }
! 147: ox_cmo_rpc(Cddgmp_Proc,"intpt",R,C,D);
! 148: Ret = ox_pop_cmo(Cdd_Procgmp);
! 149: V = vector(C-1);
! 150: for ( I = 1, D = 1; I < C; I++ ) {
! 151: V[I-1] = Ret[I];
! 152: D = ilcm(D,dn(Ret[I]));
! 153: }
! 154: V *= D;
! 155: return V;
! 156: }
! 157:
! 158: def compute_last_w(G,GH,W1,W2,W)
! 159: {
! 160: N = length(G);
! 161: for ( Dg = [], I = 0; I < N; I++ )
! 162: Dg = append(compute_df(G[I],GH[I]),Dg);
! 163: for ( V = [], T = Dg; T != []; T = cdr(T) ) {
! 164: S = car(T);
! 165: if ( in_c12(S,W1,W2) && comp_facet_preorder(W,S,W1,W2) )
! 166: V = cons(S,V);
! 167: }
! 168: if ( V == [] ) return 1;
! 169: for ( T = V; T != []; T = cdr(T) ) {
! 170: S = car(T);
! 171: for ( U = V; U != []; U = cdr(U) ) {
! 172: if ( !comp_facet_preorder(S,car(U),W1,W2) )
! 173: break;
! 174: }
! 175: if ( U == [] )
! 176: return S;
! 177: }
! 178: error("compute_last_w : cannot happen");
! 179: }
! 180:
! 181: def highest_w(G,GH,W1,W2,W)
! 182: {
! 183: N = length(G);
! 184: In = vector(N);
! 185: for ( I = 0; I < N; I++ ) {
! 186: F = G[I];
! 187: H = dp_etov(GH[I]);
! 188: for ( R = 0, T = F; T; T = dp_rest(T) ) {
! 189: S = H-dp_etov(T);
! 190: if ( comp_facet_preorder(S,W,W1,W2) && comp_facet_preorder(W,S,W1,W2) )
! 191: R += dp_hm(T);
! 192: }
! 193: In[I] = R;
! 194: }
! 195: return In;
! 196: }
! 197:
! 198: def mat_to_vlist(M)
! 199: {
! 200: N = size(M)[0];
! 201: R = vector(N);
! 202: for ( I = 0; I < N; I++ )
! 203: R[I] = M[I];
! 204: return R;
! 205: }
! 206:
! 207: def porder(F,G)
! 208: {
! 209: while ( 1 ) {
! 210: if ( !F ) {
! 211: if ( !G ) return 0;
! 212: else return -1;
! 213: } else if ( !G ) return 1;
! 214: else {
! 215: HF = dp_ht(F); HG = dp_ht(G);
! 216: if ( HF > HG ) return 1;
! 217: else if ( HF < HG ) return -1;
! 218: F = dp_rest(F); G = dp_rest(G);
! 219: }
! 220: }
! 221: }
! 222:
! 223: def nf_marked(B,G,PS,HPS)
! 224: {
! 225: M = 1;
! 226: for ( D = 0; G; ) {
! 227: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 228: if ( dp_redble(G,H=HPS[car(L)]) > 0 ) {
! 229: R = PS[car(L)];
! 230: GCD = igcd(dp_hc(G),dp_hc(H));
! 231: CG = idiv(dp_hc(H),GCD); CR = idiv(dp_hc(G),GCD);
! 232: U = CG*G-dp_subd(G,H)*CR*R;
! 233: if ( !U )
! 234: return [D,M];
! 235: D *= CG; M *= CG;
! 236: break;
! 237: }
! 238: }
! 239: if ( U ) {
! 240: G = U;
! 241: } else {
! 242: D += dp_hm(G); G = dp_rest(G);
! 243: }
! 244: }
! 245: return [D,M];
! 246: }
! 247:
! 248: def appear(H,F)
! 249: {
! 250: for ( ; F != []; F = cdr(F) )
! 251: if ( car(F) == H ) return 1;
! 252: return 0;
! 253: }
! 254:
! 255: def nf_marked2(B,G,PS,HPS)
! 256: {
! 257: M = 1;
! 258: Hist = [];
! 259: Post = 0;
! 260: for ( D = 0; G; ) {
! 261: for ( U = 0, Post1 = 0, L = B; L != []; L = cdr(L) ) {
! 262: if ( dp_redble(G,H=HPS[car(L)]) > 0 ) {
! 263: if ( appear(dp_ht(G),Hist) ) {
! 264: Post1 = dp_hm(G);
! 265: break;
! 266: }
! 267: R = PS[car(L)];
! 268: Hist = cons(dp_ht(G),Hist);
! 269: GCD = igcd(dp_hc(G),dp_hc(H));
! 270: CG = idiv(dp_hc(H),GCD); CR = idiv(dp_hc(G),GCD);
! 271: U = CG*G-dp_subd(G,H)*CR*R;
! 272: if ( !U )
! 273: return [D,Post,M];
! 274: D *= CG; M *= CG; Post *= CG;
! 275: break;
! 276: }
! 277: }
! 278: if ( U ) {
! 279: G = U;
! 280: } else if ( Post1 ) {
! 281: Post += Post1; G = dp_rest(G);
! 282: } else {
! 283: D += dp_hm(G); G = dp_rest(G);
! 284: }
! 285: }
! 286: return [D,Post,M];
! 287: }
! 288:
! 289: def nf_marked_rec(B,G,PS,HPS)
! 290: {
! 291: D = 0; M = 1;
! 292: while ( 1 ) {
! 293: L = nf_marked2(B,G,PS,HPS);
! 294: D1 = L[0]; P1 = L[1]; M1 = L[2];
! 295: /* M1*G = D1+P1 */
! 296: D = M1*D+D1;
! 297: M *= M1;
! 298: if ( !P1 ) return remove_cont([D,M]);
! 299: else G = P1;
! 300: }
! 301: }
! 302:
! 303: #if 0
! 304: /* increasing order */
! 305:
! 306: def comp_second(A,B)
! 307: {
! 308: if ( A[1] > B[1] ) return 1;
! 309: else if ( A[1] < B[1] ) return -1;
! 310: else return 0;
! 311: }
! 312: #else
! 313: /* decreasing order */
! 314:
! 315: def comp_second(A,B)
! 316: {
! 317: if ( A[1] > B[1] ) return -1;
! 318: else if ( A[1] < B[1] ) return 1;
! 319: else return 0;
! 320: }
! 321: #endif
! 322:
! 323: def sort_by_order(G,GH)
! 324: {
! 325: N = length(G);
! 326: V = vector(N);
! 327: for ( I = 0; I < N; I++ )
! 328: V[I] = [G[I],GH[I]];
! 329: V1 = qsort(V,comp_second);
! 330: for ( I = 0; I < N; I++ ) {
! 331: G[I] = V1[I][0]; GH[I] = V1[I][1];
! 332: }
! 333:
! 334: }
! 335:
! 336: def generic_walk(G1,V,M1,M2)
! 337: {
! 338: /* G is always a set of DP w.r.t. order W1 */
! 339: NV = length(V);
! 340: dp_ord(M1);
! 341: W1 = mat_to_vlist(M1);
! 342: W2 = mat_to_vlist(M2);
! 343: G = ltov(map(dp_ptod,G1,V));
! 344: GH = map(dp_hm,G);
! 345: dp_ord(M2);
! 346: G = ltov(map(dp_ptod,G1,V));
! 347:
! 348: Tw = Tgb = Tcw = Tlift = Tred = 0;
! 349: Tnf = Tcont = 0;
! 350: W = 0;
! 351: Step = 0;
! 352:
! 353: S2 = size(W2); R1 = S2[0]+1;
! 354: CW2 = matrix(R1,NV);
! 355: for ( I = 1; I < R1; I++ )
! 356: for ( J = 0; J < NV; J++ )
! 357: CW2[I][J] = W2[I-1][J];
! 358: CW = compute_compat_weight(G,GH);
! 359: for ( I = 0; I < NV; I++ )
! 360: CW2[0][I] = CW[I];
! 361:
! 362: while ( 1 ) {
! 363: print("step ",0); print(Step++);
! 364: T0 = time()[0];
! 365: L = dp_compute_last_w(vtol(G),vtol(GH),W,M1,M2);
! 366: if ( !L ) {
! 367: print(["w",Tw,"gb",Tgb,"lift",Tlift,"cw",Tcw,"red",Tred]);
! 368: return vtol(map(dp_dtop,G,V));
! 369: }
! 370: W = L[0]; In = map(dp_dtop,L[1],V);
! 371: dp_ord(CW2); InVec = ltov(map(dp_ptod,In,V));
! 372: Tw += time()[0]-T0;
! 373: T0 = time()[0];
! 374: H = nd_gr(In,V,0,M2);
! 375: N = length(H);
! 376: H1 = vector(N);
! 377: HH1 = vector(N);
! 378:
! 379: dp_ord(M2);
! 380: for ( I = 0; I < N; I++ )
! 381: HH1[I] = dp_hm(dp_ptod(H[I],V));
! 382: NG = length(G);
! 383: for ( Ind = [], I = 0; I < NG; I++ )
! 384: Ind = cons(I,Ind);
! 385: Ind = reverse(Ind);
! 386: Tgb += time()[0]-T0;
! 387:
! 388: dp_ord(CW2);
! 389: T0 = time()[0];
! 390: G = map(dp_ptod,map(dp_dtop,G,V),V);
! 391: for ( I = 0; I < N; I++ ) {
! 392: F = dp_ptod(H[I],V);
! 393: T = dp_true_nf_and_quotient_marked(Ind,F,InVec,GH);
! 394: NF = T[0]; Den = T[1]; Q = T[2];
! 395: for ( J = 0, U = 0; J < NG; J++ )
! 396: U += Q[J]*G[J];
! 397: H1[I] = U;
! 398: HH1[I] = Den*HH1[I];
! 399: }
! 400: T1 = time()[0]; Tlift += T1-T0;
! 401:
! 402: T0 = time()[0];
! 403: CW = compute_compat_weight(H1,HH1);
! 404: for ( I = 0; I < NV; I++ )
! 405: CW2[0][I] = CW[I];
! 406: dp_ord(CW2);
! 407: H1 = map(dp_ptod,map(dp_dtop,H1,V),V);
! 408:
! 409: sort_by_order(H1,HH1);
! 410: Tcw += time()[0]-T0;
! 411:
! 412: T0 = time()[0];
! 413: for ( I = 0; I < N; I++ ) {
! 414: for ( Ind = [], J = 0; J < N; J++ )
! 415: if ( J != I ) Ind = cons(J,Ind);
! 416: Ind = reverse(Ind);
! 417: T = dp_true_nf_marked(Ind,H1[I],H1,HH1);
! 418: HT = dp_ht(HH1[I]);
! 419: for ( S = S0 = T[0]; S; S = dp_rest(S) )
! 420: if ( dp_ht(S) == HT )
! 421: break;
! 422: if ( !S )
! 423: error("cannot happen");
! 424: H1[I] = S0;
! 425: HH1[I] = dp_hm(S);
! 426: }
! 427: T1 = time()[0]; Tred += T1-T0;
! 428: G = H1;
! 429: GH = HH1;
! 430: }
! 431: }
! 432:
! 433: def generic_walk_mod(G1,V,M1,M2,Mod)
! 434: {
! 435: /* G is always a set of DP w.r.t. order W1 */
! 436: setmod(Mod);
! 437: NV = length(V);
! 438: dp_ord(M1);
! 439: W1 = mat_to_vlist(M1);
! 440: W2 = mat_to_vlist(M2);
! 441: G = ltov(map(dp_mod,map(dp_ptod,G1,V),Mod,[]));
! 442: GH = map(dp_hm,G);
! 443: dp_ord(M2);
! 444: G = ltov(map(dp_mod,map(dp_ptod,G1,V),Mod,[]));
! 445:
! 446: Tw = Tgb = Tlift = Tred = 0;
! 447: Tnf = Tcont = 0;
! 448: W = 0;
! 449: Step = 0;
! 450:
! 451: S2 = size(W2); R1 = S2[0]+1;
! 452: CW2 = matrix(R1,NV);
! 453: for ( I = 1; I < R1; I++ )
! 454: for ( J = 0; J < NV; J++ )
! 455: CW2[I][J] = W2[I-1][J];
! 456: CW = compute_compat_weight(G,GH);
! 457: for ( I = 0; I < NV; I++ )
! 458: CW2[0][I] = CW[I];
! 459:
! 460: while ( 1 ) {
! 461: print("step ",0); print(Step++);
! 462: T0 = time()[0];
! 463: L = dp_compute_last_w(vtol(G),vtol(GH),W,M1,M2);
! 464: if ( !L ) {
! 465: print(["w",Tw,"gb",Tgb,"lift",Tlift,"cw",Tcw,"red",Tred]);
! 466: return vtol(map(dp_dtop,G,V));
! 467: }
! 468: W = L[0]; In = map(dp_dtop,L[1],V);
! 469: dp_ord(CW2); InVec = ltov(map(dp_mod,map(dp_ptod,In,V),Mod,[]));
! 470: Tw += time()[0]-T0;
! 471: T0 = time()[0];
! 472: H = nd_gr(In,V,Mod,M2);
! 473: N = length(H);
! 474: H1 = vector(N);
! 475: HH1 = vector(N);
! 476:
! 477: dp_ord(M2);
! 478: for ( I = 0; I < N; I++ )
! 479: HH1[I] = dp_hm(dp_mod(dp_ptod(H[I],V),Mod,[]));
! 480: NG = length(G);
! 481: for ( Ind = [], I = 0; I < NG; I++ )
! 482: Ind = cons(I,Ind);
! 483: Ind = reverse(Ind);
! 484: Tgb += time()[0]-T0;
! 485:
! 486: dp_ord(CW2);
! 487: T0 = time()[0];
! 488: G = map(dp_mod,map(dp_ptod,map(dp_dtop,G,V),V),Mod,[]);
! 489: for ( I = 0; I < N; I++ ) {
! 490: F = dp_mod(dp_ptod(H[I],V),Mod,[]);
! 491: T = dp_true_nf_and_quotient_marked_mod(Ind,F,InVec,GH,Mod);
! 492: NF = T[0]; Den = T[1]; Q = T[2];
! 493: for ( J = 0, U = 0; J < NG; J++ )
! 494: U += Q[J]*G[J];
! 495: H1[I] = U;
! 496: HH1[I] = Den*HH1[I];
! 497: }
! 498: T1 = time()[0]; Tlift += T1-T0;
! 499:
! 500: T0 = time()[0];
! 501: CW = compute_compat_weight(H1,HH1);
! 502: for ( I = 0; I < NV; I++ )
! 503: CW2[0][I] = CW[I];
! 504: dp_ord(CW2);
! 505:
! 506: H1 = map(dp_mod,map(dp_ptod,map(dp_dtop,H1,V),V),Mod,[]);
! 507: sort_by_order(H1,HH1);
! 508: Tcw += time()[0]-T0;
! 509:
! 510: T0 = time()[0];
! 511: for ( I = 0; I < N; I++ ) {
! 512: for ( Ind = [], J = 0; J < N; J++ )
! 513: if ( J != I ) Ind = cons(J,Ind);
! 514: Ind = reverse(Ind);
! 515: T = dp_true_nf_marked_mod(Ind,H1[I],H1,HH1,Mod);
! 516: HT = dp_ht(HH1[I]);
! 517: for ( S = S0 = T[0]; S; S = dp_rest(S) )
! 518: if ( dp_ht(S) == HT )
! 519: break;
! 520: if ( !S )
! 521: error("cannot happen");
! 522: H1[I] = S0;
! 523: HH1[I] = dp_hm(S);
! 524: }
! 525: T1 = time()[0]; Tred += T1-T0;
! 526: G = H1;
! 527: GH = HH1;
! 528: }
! 529: }
! 530:
! 531: def gw(G,V,O1,O2)
! 532: {
! 533: N = length(V);
! 534: W1 = create_order_mat(N,O1);
! 535: W2 = create_order_mat(N,O2);
! 536: R = generic_walk(G,V,W1,W2);
! 537: return R;
! 538: }
! 539:
! 540: def gw_mod(G,V,O1,O2,Mod)
! 541: {
! 542: N = length(V);
! 543: W1 = create_order_mat(N,O1);
! 544: W2 = create_order_mat(N,O2);
! 545: R = generic_walk_mod(G,V,W1,W2,Mod);
! 546: return R;
! 547: }
! 548:
! 549: def tolex_gw(G,V,O1)
! 550: {
! 551:
! 552: dp_ord(O1);
! 553: D1 = ltov(map(dp_ptod,G,V));
! 554: H1 = map(dp_ht,D1);
! 555: W1 = compute_compat_weight_gmp(D1,H1);
! 556:
! 557: if ( zero_dim(G,V,O1) )
! 558: G2M = tolexm(G,V,O1,V,31991);
! 559: else {
! 560: GH = map(homogenize,G,V,hhh);
! 561: VH = append(V,[hhh]);
! 562: G2M = nd_gr(GH,VH,31991,1);
! 563: G2M = map(subst,G2M,hhh,1);
! 564: G2M = nd_gr_postproc(G2M,V,31991,2,0);
! 565: }
! 566: dp_ord(2);
! 567: D2 = ltov(map(dp_ptod,G2M,V));
! 568: H2 = map(dp_ht,D2);
! 569: W2 = compute_compat_weight_gmp(D2,H2);
! 570:
! 571: N = length(V);
! 572: M1 = create_order_mat(N,0);
! 573: for ( I = 0; I < N; I++ ) M1[0][I] = W1[I];
! 574: M2 = create_order_mat(N,0);
! 575: for ( I = 0; I < N; I++ ) M2[0][I] = W2[I];
! 576: dp_set_top_weight(0);
! 577: R = generic_walk(G,V,M1,M2);
! 578: dp_set_top_weight(0);
! 579: return R;
! 580: }
! 581: endmodule$
! 582:
! 583: end$
! 584: