![[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]];
1.3 ! ohara 329: V1 = qsort(V,noro_gw.comp_second);
1.1 noro 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: endmodule$
550:
551: end$
552: