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