Annotation of OpenXM/src/asir-contrib/testing/noro/wishartd.rr, Revision 1.4
1.4 ! noro 1: /* $OpenXM: OpenXM/src/asir-contrib/testing/noro/wishartd.rr,v 1.3 2016/04/26 00:55:02 noro Exp $ */
1.1 noro 2: /* A package for computing PDEs satisfied by a matrix 1F1 on diagonal regions */
1.3 noro 3: if (version()<20160401) {error("The version of Risa/Asir must be after 20160401 to run this package.");} else{}
1.1 noro 4: module n_wishartd$
5: localf compf1$
6: localf my_comp_f$
7: localf compc1, compc, compd, compt, addpf, subpf, addc, ftorat, ctorat, mulf$
8: localf lcmf, mulf, divf, mulc, mulcpf, diffc, muldpf, normalizef1, normalizec1, normalizepf$
9: localf wsetup, mtot, wishartpf, degf, removef, subst_f, lopitalpf, simpop, simppf$
10: localf eliminatetop, diagpf, diag0pf, adjop1, adjop, reducepf, reduceotherspf$
11: localf print_f, printc, printt, printpf, ftop, ctord, comppfrd, multpf, spolypf, nfpf$
1.3 noro 12: localf nfpf0,nfpf1,nfpf_zero$
1.1 noro 13: localf substc$
14: localf dpf,dpf2,dpf3,abs$
15: localf pftord$
16: localf lopital_others_pf$
17: localf pftop$
18: localf dnc,pftozpf,addzpf,zpftopf$
19: localf mul_lopitalpf$
20: localf indep_eqs$
21: localf mulpf$
22: localf pftoeuler$
23: localf gammam,prc,prc2$
24: localf ps,psvalue,act_op,decomp_op,create_homo,bpsvalue$
25: localf pfaffian,gen_basis$
26: localf evalpf,eval_pfaffian,genrk45$
27: localf solve_leq,ptom$
28: localf msubst$
29: localf vton,encodef1,encodef,encodec1,encodec$
30: localf vton,ftotex,ctotex,pftotex,optotex,eqtotex$
31: localf genprc,prob_by_hgm,prob_by_ps$
32: localf partition, all_partition, coef_partition, count_dup,simp_power1$
33: localf lop1,lopn,loph,mul_lopitalpf2$
34: localf muldmpf$
35: localf diagcheck,eliminatecheck$
36: localf getchi$
37: localf message$
38: localf pfaffian2, eval_pfaffian2$
1.3 noro 39: localf mapleout$
40: localf gbcheck$
41: localf partition,all_partition,partition_to_pattern$
42: localf degpf,pftodpf,all_diag,y1y2$
1.4 ! noro 43: localf help$
1.1 noro 44:
45: /*
46: * pfpoly = [[C,<<...>>],...]
47: * C = [[N,L],...] L = [[F,M],...]
48: * lc(F) = 1
49: */
50:
51: static Print$
52: static PF_N,PF_V,PF_DV$
53: static Tlopital,Tred,Tlineq,Tactmul$
54: static Tp,Tps,Trk$
55:
56: /* XXX execute ord([y1,y2,...,dy1,dy2,...]) */
57:
58: /* F1>F2 <=> var(F1)>var(F2) || var(F1)=var(F2) && rest(F1)>rest(F2) */
59: /* F1^M1 > F2^M2 <=> F1>F2 || F1=F2 && M1>M2 */
60:
61: def abs(A) { return A>0 ? A : -A; }
62:
63: def compf1(E1,E2)
64: {
65: F1 = E1[0]; F2 = E2[0];
66: if ( F1>F2 ) return 1;
67: if ( F1<F2 ) return -1;
68: M1 = E1[1]; M2 = E2[1];
69: if ( M1 > M2 ) return 1;
70: if ( M1 < M2 ) return -1;
71: return 0;
72: }
73:
74: def compc1(ND1,ND2)
75: {
76: if ( R = comp_f(ND1[1],ND2[1]) ) return R;
77: N1 = ND1[0]; N2 = ND2[0];
78: if ( N1 > N2 ) return 1;
79: if ( N1 < N2 ) return -1;
80: return 0;
81: }
82:
83: /* compare ND lists */
84: def compc(L1,L2)
85: {
86: for ( ; L1 != [] && L2 != []; L1 = cdr(L1), L2 = cdr(L2) ) {
87: E1 = car(L1); E2 = car(L2);
88: if ( R = compc1(E1,E2) ) return R;
89: }
90: if ( L1 != [] ) return 1;
91: if ( L2 != [] ) return -1;
92: return 0;
93: }
94:
95: def compd(D1,D2)
96: {
97: if ( D1 > D2 ) return 1;
98: if ( D1 < D2 ) return -1;
99: return 0;
100: }
101:
102: /* compare terms */
103: def compt(T1,T2)
104: {
105: if ( R = compd(T1[1],T2[1]) ) return R;
106: if ( R = compc(T1[0],T2[0]) ) return R;
107: return 0;
108: }
109:
110: def addpf(F1,F2)
111: {
112: R = [];
113: while ( F1 != [] && F2 != [] ) {
114: T1 = car(F1); T2 = car(F2);
115: C = compd(T1[1],T2[1]);
116: if ( C > 0 ) {
117: R = cons(T1,R); F1 = cdr(F1);
118: } else if ( C < 0 ) {
119: R = cons(T2,R); F2 = cdr(F2);
120: } else {
121: S = addc(T1[0],T2[0]);
122: if ( S != [] )
123: R = cons([S,T1[1]],R);
124: F1 = cdr(F1); F2 = cdr(F2);
125: }
126: }
127: R = reverse(R);
128: if ( F1 != [] ) R = append(R,F1);
129: else if ( F2 != [] ) R = append(R,F2);
130: return R;
131: }
132:
133: def subpf(F1,F2)
134: {
135: T = mulcpf([[-1,[]]],F2);
136: T = addpf(F1,T);
137: return T;
138: }
139:
140: def addc(F1,F2)
141: {
142: R = [];
143: while ( F1 != [] && F2 != [] ) {
144: T1 = car(F1); T2 = car(F2);
145: C = comp_f(T1[1],T2[1]);
146: if ( !T1[0] || !T2[0] )
147: error("afo");
148: if ( C > 0 ) {
149: R = cons(T1,R); F1 = cdr(F1);
150: } else if ( C < 0 ) {
151: R = cons(T2,R); F2 = cdr(F2);
152: } else {
153: S = T1[0]+T2[0];
154: if ( S )
155: R = cons([S,T1[1]],R);
156: F1 = cdr(F1); F2 = cdr(F2);
157: }
158: }
159: R = reverse(R);
160: if ( F1 != [] ) R = append(R,F1);
161: else if ( F2 != [] ) R = append(R,F2);
162: return R;
163: }
164:
165: def ftorat(F)
166: {
167: R = 1;
168: for ( ; F != []; F = cdr(F) ) {
169: F0 = car(F);
170: R *= F0[0]^F0[1];
171: }
172: return R;
173: }
174:
175: def ctorat(C)
176: {
177: R = 0;
178: for ( ; C != []; C = cdr(C) ) {
179: C0 = car(C);
180: R += C0[0]/ftorat(C0[1]);
181: R = red(R);
182: }
183: return R;
184: }
185:
186: def mulf(L1,L2)
187: {
188: R = [];
189: while ( L1 != [] && L2 != [] ) {
190: A1 = car(L1); A2 = car(L2);
191: if ( A1[0] > A2[0] ) {
192: R = cons(A1,R); L1 = cdr(L1);
193: } else if ( A1[0] < A2[0] ) {
194: R = cons(A2,R); L2 = cdr(L2);
195: } else {
196: R = cons([A1[0],A1[1]+A2[1]],R);
197: L1 = cdr(L1); L2 = cdr(L2);
198: }
199: }
200: R = reverse(R);
201: if ( L2 == [] ) return append(R,L1);
202: else if ( L1 == [] ) return append(R,L2);
203: else return R;
204: }
205:
206: def lcmf(L1,L2)
207: {
208: R = [];
209: while ( L1 != [] && L2 != [] ) {
210: A1 = car(L1); A2 = car(L2);
211: if ( A1[0] > A2[0] ) {
212: R = cons(A1,R); L1 = cdr(L1);
213: } else if ( A1[0] < A2[0] ) {
214: R = cons(A2,R); L2 = cdr(L2);
215: } else {
216: M = A1[1]>A2[1]?A1[1]:A2[1];
217: R = cons([A1[0],M],R);
218: L1 = cdr(L1); L2 = cdr(L2);
219: }
220: }
221: R = reverse(R);
222: if ( L2 == [] ) return append(R,L1);
223: else if ( L1 == [] ) return append(R,L2);
224: else return R;
225: }
226:
227: def mulf(L1,L2)
228: {
229: R = [];
230: while ( L1 != [] && L2 != [] ) {
231: A1 = car(L1); A2 = car(L2);
232: if ( A1[0] > A2[0] ) {
233: R = cons(A1,R); L1 = cdr(L1);
234: } else if ( A1[0] < A2[0] ) {
235: R = cons(A2,R); L2 = cdr(L2);
236: } else {
237: R = cons([A1[0],A1[1]+A2[1]],R);
238: L1 = cdr(L1); L2 = cdr(L2);
239: }
240: }
241: R = reverse(R);
242: if ( L2 == [] ) return append(R,L1);
243: else if ( L1 == [] ) return append(R,L2);
244: else return R;
245: }
246:
247: def divf(L1,L2)
248: {
249: R = [];
250: while ( L1 != [] && L2 != [] ) {
251: A1 = car(L1); A2 = car(L2);
252: if ( A1[0] > A2[0] ) {
253: R = cons(A1,R); L1 = cdr(L1);
254: } else if ( A1[0] < A2[0] ) {
255: error("divf : cannot happen");
256: } else {
257: M = A1[1]-A2[1];
258: if ( M < 0 )
259: error("divf : cannot happen");
260: if ( M > 0 )
261: R = cons([A1[0],M],R);
262: L1 = cdr(L1); L2 = cdr(L2);
263: }
264: }
265: R = reverse(R);
266: if ( L2 == [] ) return append(R,L1);
267: else if ( L1 == [] ) return append(R,L2);
268: else return R;
269: }
270:
271: def mulc(C1,C2)
272: {
273: R = [];
274: C1 = reverse(C1);
275: for ( ; C1 != []; C1 = cdr(C1) ) {
276: S = [];
277: A1 = car(C1);
278: for ( T = C2 ; T != []; T = cdr(T) ) {
279: A2 = car(T);
280: S = cons([A1[0]*A2[0],mulf(A1[1],A2[1])],S);
281: }
282: S = reverse(S);
283: R = addc(R,S);
284: }
285: return R;
286: }
287:
288: def vton(V)
289: {
290: for ( I = 1; I <= PF_N; I++ )
291: if ( V == PF_V[I] ) return I;
292: error("vton : no such variable");
293: }
294:
295: def encodef1(F)
296: {
297: A = F[0];
298: V1 = var(A);
299: N1 = vton(V1);
300: R = A-V1;
301: if ( !R )
302: return [N1,F[1]];
303: else
304: return [N1*65536+vton(var(R)),F[1]];
305: }
306:
307: def encodef(F)
308: {
309: return map(encodef1,F);
310: }
311:
312: def encodec1(C)
313: {
314: return cons(C[0],encodef(C[1]));
315: }
316:
317: def encodec(C)
318: {
319: R = map(encodec1,C);
320: }
321:
322: def mulcpf(C,F)
323: {
324: R = [];
325: for ( ; F != []; F = cdr(F) ) {
326: F0 = car(F);
327: C1 = mulc(C,F0[0]);
328: R = cons([C1,F0[1]],R);
329: }
330: return reverse(R);
331: }
332:
333: def mulpf(F1,F2)
334: {
335: R = [];
336: N = length(PF_V);
337: for ( T = reverse(F1); T != []; T = cdr(T) ) {
338: T0 = car(T); C = T0[0]; Op = T0[1];
339: E = dp_etov(Op);
340: S = F2;
341: for ( I = 0; I < N; I++ )
342: for ( J = 0; J < E[I]; J++ ) S = muldpf(PF_V[I],S);
343: S = mulcpf(C,S);
344: R = addpf(R,S);
345: }
346: return S;
347: }
348:
349: def diffc(X,C)
350: {
351: R = [];
352: for ( T = C; T != []; T = cdr(T) ) {
353: T0 = car(T);
354: N = T0[0];
355: S = [];
356: for ( L = T0[1]; L != []; S = cons(car(L),S), L = cdr(L) ) {
357: L0 = car(L);
358: F = L0[0]; M = L0[1];
359: DF = diff(F,X);
360: if ( DF ) {
361: V = cons([F,M+1],cdr(L));
362: for ( U = S; U != []; U = cdr(U) ) V = cons(car(U),V);
363: R = addc([[-M*N*DF,V]],R);
364: }
365: }
366: }
367: return R;
368: }
369:
370: def muldpf(X,F)
371: {
372: R = [];
373: DX = dp_ptod(strtov("d"+rtostr(X)),PF_DV);
374: for ( T = F; T != []; T = cdr(T) ) {
375: T0 = car(T);
376: C = diffc(X,T0[0]);
377: if ( C != [] )
378: R = addpf(R,[[T0[0],T0[1]*DX],[C,T0[1]]]);
379: else
380: R = addpf(R,[[T0[0],T0[1]*DX]]);
381: }
382: return R;
383: }
384:
385: /* d^m*c = sum_{i=0}^m m!/i!/(m-i)!*c^(i)*d^(m-i) */
386: def muldmpf(X,M,F)
387: {
388: R = [];
389: DX = strtov("d"+rtostr(X));
390: FM = fac(M);
391: for ( T = F; T != []; T = cdr(T) ) {
392: T0 = car(T);
393: C = T0[0]; Op = T0[1];
394: U = [];
395: for ( I = 0; I <= M; I++ ) {
396: if ( C == [] ) break;
397: C0 = [[FM/fac(I)/fac(M-I),[]]];
398: U = cons([mulc(C0,C),dp_ptod(DX^(M-I),PF_DV)*Op],U);
399: C = diffc(X,C);
400: }
401: U = reverse(U);
402: R = addpf(U,R);
403: }
404: return R;
405: }
406:
407: def normalizef1(F)
408: {
409: if ( coef(F,1,var(F)) < 0 ) F = -F;
410: return F;
411: }
412:
413: def normalizec1(C)
414: {
415: N = C[0]; L = C[1];
416: S = 1;
417: R = [];
418: for ( ; L != []; L = cdr(L) ) {
419: L0 = car(L);
420: F = L0[0]; M = L0[1];
421: if ( coef(F,1,var(F)) < 0 ) {
422: F = -F;
423: if ( M%2 ) S = -S;
424: }
425: R = cons([F,M],R);
426: }
427: return [S*N,reverse(qsort(R,n_wishartd.compf1))];
428: }
429:
430: def normalizepf(F)
431: {
432: R = [];
433: for ( ; F != []; F = cdr(F) ) {
434: F0 = car(F);
435: C = map(normalizec1,F[0]);
436: R = cons([C,F[1]],R);
437: }
438: return reverse(R);
439: }
440:
441: def substc(C,Y2,Y1)
442: {
443: R = [];
444: for ( T = C; T != []; T = cdr(T) ) {
445: T0 = car(T); N = T0[0]; D = T0[1];
446: Z = subst_f(D,Y2,Y1);
447: R = addc([[Z[0]*N,Z[1]]],R);
448: }
449: return R;
450: }
451:
452: /* dy0 : dummy variable for dy */
453:
454: def wsetup(N)
455: {
456: V = [];
457: DV = [];
458: for ( I = N; I>= 0; I-- ) {
459: YI = strtov("y"+rtostr(I));
460: DYI = strtov("dy"+rtostr(I));
461: V = cons(YI,V);
462: DV = cons(DYI,DV);
463: }
464: PF_N = N;
465: PF_V = V;
466: PF_DV = DV;
467: ord(append(V,DV));
468: }
469:
470: def mtot(M,Dn)
471: {
472: D = dp_ptod(M,PF_DV);
473: F = fctr(Dn); C = car(F)[0];
474: Dn = reverse(qsort(cdr(F),n_wishartd.compf1));
475: return [[[dp_hc(D)/C,Dn]],dp_ht(D)];
476: }
477:
478: def wishartpf(N,I)
479: {
480: YI = PF_V[I]; DYI = PF_DV[I];
481: R = [mtot(DYI^2,1)];
482: R = addpf([mtot(c*DYI,YI)],R);
483: R = addpf([mtot(-DYI,1)],R);
484: for ( J = 1; J <= N; J++ ) {
485: if ( J == I ) continue;
486: YJ = PF_V[J]; DYJ = PF_DV[J];
487: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
488: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
489: R = addpf([mtot(-1/2*DYI,YI)],R);
490: R = addpf([mtot(1/2*DYJ,YI)],R);
491: }
492: R = addpf([mtot(-a,YI)],R);
493: return R;
494: }
495:
496: def degf(F,D)
497: {
498: for ( ; F != []; F = cdr(F) ) {
499: F0 = car(F);
500: if ( F0[0] == D ) return F0[1];
501: }
502: return 0;
503: }
504:
505: def removef(F,D)
506: {
507: R = [];
508: for ( ; F != []; F = cdr(F) ) {
509: F0 = car(F);
510: if ( F0[0] != D ) R = cons(F0,R);
511: }
512: return reverse(R);
513: }
514:
515: def subst_f(F,Y2,Y1)
516: {
517: R = [];
518: Sgn = 1;
519: for ( ; F != []; F = cdr(F) ) {
520: F0 = car(F);
521: T = subst(F0[0],Y2,Y1);
522: if ( coef(T,1,var(T)) < 0 ) {
523: T = -T;
524: if ( F0[1]%2 ) Sgn = -Sgn;
525: }
526: R = cons([T,F0[1]],R);
527: }
528: if ( R == [] ) return [Sgn,R];
529: R = qsort(R,n_wishartd.compf1);
530: S = [car(R)]; R = cdr(R);
531: while ( R != [] ) {
532: R0 = car(R);
533: S0 = car(S);
534: if ( R0[0] == S0[0] )
535: S = cons([S0[0],S0[1]+R0[1]],cdr(S));
536: else
537: S = cons(R0,S);
538: R = cdr(R);
539: }
540: return [Sgn,S];
541: }
542:
543: /* Y2 -> Y1 */
544: def lopitalpf(P,Y1,Y2)
545: {
546: // if ( !member(Y2,vars(P)) ) return P;
547: D = normalizef1(Y2-Y1);
548: if ( D == Y2-Y1 ) Sgn = 1;
549: else Sgn = -1;
550: DY2 = strtov("d"+rtostr(Y2));
551: R = [];
552: for ( T = reverse(P); T != []; T = cdr(T) ) {
553: T0 = car(T);
554: C = T0[0]; Op = T0[1];
555: for ( U = reverse(C); U != []; U = cdr(U) ) {
556: U0 = car(U);
557: Nm = U0[0]; Dn = U0[1];
558: K = degf(Dn,D);
559: if ( !K ) {
560: Z = subst_f(Dn,Y2,Y1);
561: Sgn1 = Z[0]; Dn1 = Z[1];
562: R = addpf([[[[Sgn1*Nm,Dn1]],Op]],R);
563: } else {
564: Dn1 = removef(Dn,D);
565: C1 = [[1,Dn1]];
566: Z = [];
567: for ( J = 0; J <= K; J++ ) {
568: B = fac(K)/fac(J)/fac(K-J);
569: C1s = substc(C1,Y2,Y1);
570: if ( C1s != [] ) {
571: W = [[C1s,dp_ptod(DY2^(K-J),PF_DV)*Op]];
572: W = mulcpf([[B,[]]],W);
573: Z = addpf(W,Z);
574: }
575: C1 = diffc(Y2,C1);
576: if ( C1 == [] ) break;
577: }
578: Z = mulcpf([[Sgn^K*Nm/fac(K),[]]],Z);
579: R = addpf(Z,R);
580: }
581: }
582: }
583: return R;
584: }
585:
586: def lopital_others_pf(P,L) {
587: Q = P;
588: for ( T = L; T != []; T = cdr(T) ) {
589: T0 = car(T);
590: I0 = T0[0]; I1 = T0[1]; I = I1-I0+1;
591: for ( M = I0; M <= I1; M++ ) {
592: Q = lopitalpf(Q,PF_V[I0],PF_V[M]);
593: }
594: Q = simppf(Q,I0,I);
595: Q = adjop(Q,I0,I);
596: }
597: return Q;
598: }
599:
600: def simpop(Op,I0,NV)
601: {
602: E = dp_etov(Op);
603: R = [];
604: for ( J = 0, I = I0; J < NV; I++, J++ ) R = cons(E[I],R);
605: R = reverse(qsort(R));
606: for ( J = 0, I = I0; J < NV; I++, J++ ) E[I] = R[J];
607: return dp_vtoe(E);
608: }
609:
610: def simppf(P,I0,NV)
611: {
612: R = [];
613: for ( T = P; T != []; T = cdr(T) ) {
614: T0 = car(T);
615: R = addpf([[T0[0],simpop(T0[1],I0,NV)]],R);
616: }
617: return R;
618: }
619:
620: /* simplify (dy1+...+dyn)^k */
621:
622: def simp_power1(K,I0,NV)
623: {
624: P = all_partition(K,NV);
625: M = map(coef_partition,P,K,NV,I0);
626: for ( R = 0, T = M; T != []; T = cdr(T) )
627: R += car(T);
628: return R;
629: }
630:
631: def indep_eqs(Eq)
632: {
633: M = length(Eq);
634: D = 0;
635: for ( I = 0; I < M; I++ )
636: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
637: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
638: Terms = vector(N);
639: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
640: A = matrix(M,N);
641: for ( I = 0; I < M; I++ ) {
642: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
643: T = car(H)[1];
644: for ( K = 0; K < N; K++ )
645: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
646: }
647: }
648: for ( I = 0; I < M; I++ ) {
649: Dn = 1;
650: for ( J = 0; J < N; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
651: for ( J = 0; J < N; J++ ) A[I][J] *= Dn;
652: }
653: R = indep_rows_mod(A,lprime(0));
654: if ( length(R) == N ) {
655: L = [];
656: for ( I = N-1; I >= 0; I-- )
657: L = cons(Eq[R[I]],L);
658: return L;
659: } else
660: return 0;
661: }
662:
663: def eliminatetop(Eq)
664: {
665: Eq = indep_eqs(Eq);
666: if ( !Eq ) return 0;
667:
668: M = length(Eq);
669: R = vector(M);
670: D = 0;
671: for ( I = 0; I < M; I++ )
672: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
673: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
674: if ( N != M )
675: return 0;
676: N2 = 2*N;
677: Terms = vector(N);
678: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
679: A = matrix(N,N2);
680: for ( I = 0; I < N; I++ ) {
681: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
682: T = car(H)[1];
683: for ( K = 0; K < N; K++ )
684: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
685: }
686: A[I][N+I] = 1;
687: }
688: for ( I = 0; I < N; I++ ) {
689: Dn = 1;
690: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
691: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
692: }
693: Ret = generic_gauss_elim(A);
694: Num = Ret[0]; Den = Ret[1]; Ind0 = Ret[2]; Ind1 = Ret[3];
695: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
696: return 0;
697: R = [];
698: if ( Print ) print(["N=",N]);
699: for ( I = 0; I < N; I++ ) {
700: if ( Print ) print(["I=",I],2);
701: T = [];
702: for ( L = 0; L < N; L++ )
703: if ( Num[I][L] ) T = addpf(T,mulcpf([[Num[I][L]/Den,[]]],Eq[L][1]));
704: R = cons([Terms[I],T],R);
705: }
706: if ( Print ) print("");
707: R = qsort(R,n_wishartd.compt);
708: return reverse(R);
709: }
710:
711: def eliminatecheck(Eq,Top)
712: {
713: Eq = indep_eqs(Eq);
714: if ( !Eq ) return 0;
715:
716: M = length(Eq);
717: R = vector(M);
718: D = 0;
719: for ( I = 0; I < M; I++ )
720: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
721: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
722: if ( N != M )
723: return 0;
724: N2 = 2*N;
725: Terms = vector(N);
726: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
727: A = matrix(N,N2);
728: for ( I = 0; I < N; I++ ) {
729: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
730: T = car(H)[1];
731: for ( K = 0; K < N; K++ )
732: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
733: }
734: A[I][N+I] = 1;
735: }
736: for ( I = 0; I < N; I++ ) {
737: Dn = 1;
738: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
739: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
740: }
741: if ( Top ) {
742: for ( I = 0; I < N; I++ ) {
743: for ( T = J = 0; J < N; J++ ) if ( J != I ) T += abs(A[I][J]);
744: if ( abs(A[I][I]) < T )
745: if ( Print ) print(["ng",I]);
746: }
747: } else {
748: for ( I = 1; I < N; I++ ) {
749: for ( T = J = 0; J < N-2; J++ ) if ( J != I-1 ) T += abs(A[I][J]);
750: if ( abs(A[I][I-1]) < T )
751: if ( Print ) print(["ng",I]);
752: }
753: }
754:
755: Ret = generic_gauss_elim_mod(A,lprime(0));
756: Num = Ret[0]; Ind0 = Ret[1]; Ind1 = Ret[2];
757: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
758: return 0;
759: else
760: return 1;
761: }
762:
763: def mul_lopitalpf(Z,N,K,I0,I1,E)
764: {
765: I = I1-I0+1;
766: for ( J = I0; J <= N; J++ )
767: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
768: for ( J = I0+1; J <= I1; J++ ) {
769: Z = lopitalpf(Z,PF_V[I0],PF_V[J]);
770: }
771: S = simppf(Z,I0,I);
772: H = []; L = [];
773: for ( ; S != []; S = cdr(S) ) {
774: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
775: else L = cons(car(S),L);
776: }
777: return [reverse(H),reverse(L)];
778: }
779:
780: /* Blocks = [B1,B2,...] */
781: /* Bi=[s,e] ys=...=ye f */
782:
783: def diag0pf(N,Blocks)
784: {
785: Tlopital = Tred = Tlineq = 0;
786: wsetup(N);
787: P = vector(N+1);
788: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
789: Simp = [];
790: Done = [];
791: Len = length(Blocks);
792: for ( A = 0; A < Len; A++ ) {
793: B = Blocks[A];
794: Blocks0 = setminus(Blocks,[B]);
795: I0 = B[0];
796: I1 = B[1];
797: I = I1-I0+1;
798: for ( K = I0; K <= I1; K++ )
799: P[K] = lopital_others_pf(P[K],Blocks0);
800: Rall = [];
801: for ( K = I+1; K >= 2; K-- ) {
802: if ( Print ) print(["K=",K],2);
803: DK = simp_power1(K,I0,I);
804: Eq = [];
805: TlopitalK = 0;
806: for ( T = DK; T; T = dp_rest(T) ) {
807: E = dp_etov(T);
808: for ( U = I0; U <= I1; U++ )
809: if ( E[U] >= 2 )
810: break;
811: if ( U > I1 ) continue;
812: E[U] -= 2;
813: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
814: Eq = cons(Ret,Eq);
815: }
816: Eq = reverse(Eq);
817:
818: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
819: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
820:
821: DY = mtot(-dy0^K,1);
822: if ( K == I+1 ) {
823: EqTop = addpf([DY],Eq0);
824: } else {
825: H = []; L = [];
826: for ( S = Eq0 ; S != []; S = cdr(S) ) {
827: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
828: else L = cons(car(S),L);
829: }
830: L = reverse(L); H = reverse(H);
831: Eq = cons([H,addpf([DY],L)],Eq);
832: }
833: T0 = time();
834: R = eliminatetop(Eq);
835: if ( R )
836: Rall = cons(R,Rall);
837: else
838: return [];
839: T1 = time(); Tlineq += T1[0]-T0[0];
840: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
841: }
842: Rall = reverse(Rall);
843:
844: /* EqTop : dyi0 -> dy -> dyi0 */
845: for ( T = Rall; T != []; T = cdr(T) ) {
846: if ( Print ) print(["red",length(T)+1],2);
847: T0 = time();
848: EqTop = reducepf(EqTop,car(T));
849: T1 = time(); Tred += T1[0]-T0[0];
850: if ( Print ) print([T1[0]-T0[0],"sec"]);
851: }
852: EqTop = adjop(EqTop,I0,I);
853: Done = cons(EqTop,Done);
854:
855: }
856: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
857: Done = ltov(Done);
858: Len = length(Done);
859: for ( I = 0; I < Len; I++ ) {
860: for ( G = [], J = 0; J < Len; J++ )
861: if ( J != I ) G = cons(Done[J],G);
862: Done[I] = nfpf(Done[I],G);
863: }
864: return vtol(Done);
865: }
866:
867: def diagpf(N,Blocks)
868: {
869: Tlopital = Tred = Tlineq = 0;
870: wsetup(N);
871: P = vector(N+1);
872: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
873: Simp = [];
874: Done = [];
875: Len = length(Blocks);
876: for ( A = 0; A < Len; A++ ) {
877: B = Blocks[A];
878: Blocks0 = setminus(Blocks,[B]);
879: I0 = B[0];
880: I1 = B[1];
881: I = I1-I0+1;
882: for ( K = I0; K <= I1; K++ )
883: P[K] = lopital_others_pf(P[K],Blocks0);
884: Rall = [];
885: for ( K = I+1; K >= 2; K-- ) {
886: if ( Print ) print(["K=",K],2);
887: DK = simp_power1(K,I0,I);
888: Eq = [];
889: TlopitalK = 0;
890: for ( T = DK; T; T = dp_rest(T) ) {
891: E = dp_etov(T);
892: for ( U = I1; U >= I0; U-- )
893: if ( E[U] >= 2 )
894: break;
895: if ( U < I0 ) continue;
896: E[U] -= 2;
897: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
898: Eq = cons(Ret,Eq);
899: }
900: Eq = reverse(Eq);
901:
902: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
903: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
904:
905: DY = mtot(-dy0^K,1);
906: if ( K == I+1 ) {
907: EqTop = addpf([DY],Eq0);
908: } else {
909: H = []; L = [];
910: for ( S = Eq0 ; S != []; S = cdr(S) ) {
911: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
912: else L = cons(car(S),L);
913: }
914: L = reverse(L); H = reverse(H);
915: Eq = cons([H,addpf([DY],L)],Eq);
916: }
917: T0 = time();
918: R = eliminatetop(Eq);
919: if ( R )
920: Rall = cons(R,Rall);
921: else {
922: if ( Print ) print(["lineq retry.."],2);
923: Eq1 = [];
924: Terms = [];
925: for ( T = dp_rest(DK); T; T = dp_rest(T) )
926: Terms = cons(dp_ht(T),Terms);
927: while ( Terms != [] ) {
928: for ( Q = 0; Terms != [] && Q < 10; Q++, Terms = cdr(Terms) ) {
929: E = dp_etov(car(Terms));
930: if ( E[I0] >= 2 ) {
931: E[I0] -= 2;
932: Ret = mul_lopitalpf(P[I0],N,K,I0,I1,E);
933: Eq1 = cons(Ret,Eq1);
934: }
935: }
936: Eq = append(Eq,Eq1);
937: R = eliminatetop(Eq);
938: if ( R ) break;
939: }
940: if ( !R )
941: error("diagpf : failed to find a solvable linear system");
942: Rall = cons(R,Rall);
943: }
944: T1 = time(); Tlineq += T1[0]-T0[0];
945: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
946: }
947: Rall = reverse(Rall);
948:
949: /* EqTop : dyi0 -> dy -> dyi0 */
950: for ( T = Rall; T != []; T = cdr(T) ) {
951: if ( Print ) print(["red",length(T)+1],2);
952: T0 = time();
953: EqTop = reducepf(EqTop,car(T));
954: T1 = time(); Tred += T1[0]-T0[0];
955: if ( Print ) print([T1[0]-T0[0],"sec"]);
956: }
957: EqTop = adjop(EqTop,I0,I);
958: Done = cons(EqTop,Done);
959:
960: }
961: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
962: Done = ltov(Done);
963: Len = length(Done);
964: for ( I = 0; I < Len; I++ ) {
965: for ( G = [], J = 0; J < Len; J++ )
966: if ( J != I ) G = cons(Done[J],G);
1.3 noro 967: Done[I] = nfpf1(Done[I],G);
1.1 noro 968: }
969: return vtol(Done);
970: }
971:
972: def diagcheck(N)
973: {
974: Tmuld = Tlopital = 0;
975: MHist = [];
976: wsetup(N);
977: P = vector(N+1);
978: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
979: Simp = [];
980: Done = [];
981: B = [1,N];
982: I0 = B[0];
983: I1 = B[1];
984: I = I1-I0+1;
985: for ( K = 2; K <= I+1; K++ ) {
986: if ( Print ) print(["K=",K],2);
987: DK = simp_power1(K,I0,I);
988: Eq = [];
989: TlopitalK = 0;
990: for ( T = DK; T; T = dp_rest(T) ) {
991: E = dp_etov(T);
992: for ( U = I0; U <= I1; U++ )
993: if ( E[U] >= 2 )
994: break;
995: if ( U > I1 ) continue;
996: E[U] -= 2;
997: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E|high=1);
998: Eq = cons(Ret,Eq);
999: }
1000: Eq = reverse(Eq);
1001:
1002: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
1003: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
1004:
1005: DY = mtot(-dy0^K,1);
1006: if ( K == I+1 ) {
1007: EqTop = addpf([DY],Eq0);
1008: } else {
1009: H = []; L = [];
1010: for ( S = Eq0 ; S != []; S = cdr(S) ) {
1011: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
1012: else L = cons(car(S),L);
1013: }
1014: L = reverse(L); H = reverse(H);
1015: Eq = cons([H,addpf([DY],L)],Eq);
1016: }
1017: R = eliminatecheck(Eq,K==I+1);
1018: if ( !R )
1019: return 0;
1020: }
1021: if ( Print ) print(["muld",Tmuld,"lopital",Tlopital]);
1022: return 1;
1023: }
1024:
1025: /* dyi -> dyi/K, dy->dyi */
1026: def adjop1(T,I,K)
1027: {
1028: C = T[0];
1029: Op = dp_etov(T[1]);
1030: if ( Op[I] ) C = mulc([[1/K,[]]],C);
1031: Op[I] += Op[0];
1032: Op[0] = 0;
1033: return [C,dp_vtoe(Op)];
1034: }
1035:
1036: def adjop(P,I,K)
1037: {
1038: P = map(adjop1,P,I,K);
1039: P = qsort(P,n_wishartd.compt);
1040: return reverse(P);
1041: }
1042:
1043: def dnc(C)
1044: {
1045: D = 1;
1046: for ( T = C; T != []; T = cdr(T) ) {
1047: T0 = car(T); N = T0[0];
1048: M = sdiv(ptozp(N),N);
1049: D = ilcm(D,M);
1050: }
1051: return D;
1052: }
1053:
1054: def pftozpf(F)
1055: {
1056: D = 1;
1057: for ( T = F; T != []; T = cdr(T) ) {
1058: T0 = car(T);
1059: D = ilcm(D,dnc(T0[0]));
1060: }
1061: return [D,mulcpf([[D,[]]],F)];
1062: }
1063:
1064: def zpftopf(F)
1065: {
1066: return mulcpf([[1/F[0],[]]],F[1]);
1067: }
1068:
1069: def addzpf(F1,F2)
1070: {
1071: D1 = F1[0]; D2 = F2[0];
1072: L = ilcm(D1,D2); C1 = L/D1; C2 = L/D2;
1073: N = addpf(mulcpf([[L/D1,[]]],F1[1]),mulcpf([[L/D2,[]]],F2[1]));
1074: return [L,N];
1075: }
1076:
1077: /* R : reducers */
1078: def reducepf(Eq,R)
1079: {
1080: Ret = pftozpf(Eq); Eq = Ret[1]; Dn = Ret[0];
1081: S = [1,[]];
1082: for ( T = R, U = []; T != []; T = cdr(T) )
1083: U = cons([car(T)[0],pftozpf(car(T)[1])],U);
1084: R = reverse(U);
1085: for ( T = reverse(Eq); T != []; T = cdr(T) ) {
1086: T0 = car(T); C = T0[0]; Op = T0[1];
1087: for ( U = (R); U != []; U = cdr(U) ) {
1088: U0 = car(U);
1089: if ( dp_redble(Op,U0[0]) ) {
1090: E = dp_etov(dp_subd(Op,U0[0]));
1091: Red = U0[1];
1092: N = length(E);
1093: for ( J = 1; J < N; J++ )
1094: for ( L = 0; L < E[J]; L++ ) Red = [Red[0],muldpf(PF_V[J],Red[1])];
1095: Red = [Red[0],mulcpf(C,Red[1])];
1096: Red = [Red[0],mulcpf([[-1,[]]],Red[1])];
1097: S = addzpf(S,Red);
1098: break;
1099: }
1100: }
1101: if ( U == [] ) S = addzpf([1,[T0]],S);
1102: }
1103: S = [S[0]*Dn,S[1]];
1104: return zpftopf(S);
1105: }
1106:
1107: def reduceotherspf(Eq,P,I1,N)
1108: {
1109: S = [];
1110: Z = Eq;
1111: while ( Z != [] ) {
1112: T0 = car(Z); C = T0[0]; Op = T0[1];
1113: for ( I = I1; I <= N; I++ ) {
1114: U0 = P[I];
1115: M = car(U0)[0][0][0];
1116: H = car(U0)[1];
1117: if ( dp_redble(Op,H) ) {
1118: E = dp_etov(dp_subd(Op,H));
1119: Red = U0;
1120: Len = length(E);
1121: for ( J = 0; J < Len; J++ )
1122: for ( L = 0; L < E[J]; L++ ) Red = muldpf(PF_V[J],Red);
1123: C1 = mulc([[1/M,[]]],C);
1124: Red = mulcpf(C1,Red);
1125: Z = subpf(Z,Red);
1126: break;
1127: }
1128: }
1129: if ( I > N ) {
1130: S = cons(T0,S);
1131: Z = cdr(Z);
1132: }
1133: }
1134: return reverse(S);
1135: }
1136:
1137: def print_f(F)
1138: {
1139: for ( ; F != []; F = cdr(F) ) {
1140: F0 = car(F);
1141: print("(",0); print(F0[0],0); print(")",0);
1142: if ( F0[1] > 1 ) {
1143: print("^",0); print(F0[1],0);
1144: }
1145: if ( cdr(F) != [] ) print("*",0);
1146: }
1147: }
1148:
1149: def printc(C)
1150: {
1151: print("(",0);
1152: for ( ; C != []; C = cdr(C) ) {
1153: C0 = car(C);
1154: print("(",0); print(C0[0],0); print(")",0);
1155: if ( C0[1] != [] ) {
1156: print("/(",0);
1157: print_f(C0[1]);
1158: print(")",0);
1159: }
1160: if ( cdr(C) != [] ) print("+",0);
1161: }
1162: print(")",0);
1163: }
1164:
1165: def printt(T)
1166: {
1167: printc(T[0]); print("*",0); print(dp_dtop(T[1],PF_DV),0);
1168: }
1169:
1170: def printpf(F)
1171: {
1172: for ( ; F != []; F = cdr(F) ) {
1173: printt(car(F));
1174: if ( cdr(F) != [] ) print("+",0);
1175: }
1176: print("");
1177: }
1178:
1179: def ftop(F)
1180: {
1181: P = 1;
1182: for ( ; F != []; F = cdr(F) ) {
1183: F0 = car(F);
1184: P *= F0[0]^F0[1];
1185: }
1186: return P;
1187: }
1188:
1189: def pftord(F)
1190: {
1191: R = [];
1192: for ( T = F; T != []; T = cdr(T) ) {
1193: T0 = car(T); C = T0[0]; Op = T0[1];
1194: R = cons([ctord(C),Op],R);
1195: }
1196: return reverse(R);
1197: }
1198:
1199: def pftop(F)
1200: {
1201: R = pftord(F);
1202: Op = F[0][1];
1203: N = length(dp_etov(Op));
1204: DY = [];
1205: for ( I = N; I >= 0; I-- ) DY = cons(strtov("dy"+rtostr(I)),DY);
1206: Lcm = [];
1207: for ( T = R; T != []; T = cdr(T) )
1208: Lcm = lcmf(Lcm,T[0][0][1]);
1209: S = 0;
1210: for ( T = R; T != []; T = cdr(T) ) {
1211: N = T[0][0][0]*ftop(divf(Lcm,T[0][0][1]));
1212: Op = dp_dtop(T[0][1],DY);
1213: S += N*Op;
1214: }
1215: return S;
1216: }
1217:
1218: def ctord(C)
1219: {
1220: N = 0; D = [];
1.3 noro 1221: Len = length(C);
1222: for ( I = 0, T = reverse(C); T != []; T = cdr(T),I++ ) {
1223: // print([Len-I],2);
1.1 noro 1224: T0 = car(T); N0 = T0[0]; D0 = T0[1];
1225: L = lcmf(D,D0);
1226: /* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
1227: N = N*ftop(divf(L,D))+N0*ftop(divf(L,D0));
1228: D = L;
1229: }
1.3 noro 1230:
1.1 noro 1231: R = [];
1.3 noro 1232: for ( S = D; S != []; S = cdr(S) ) {
1233: T0 = car(S); F = T0[0]; M = T0[1];
1.1 noro 1234: while ( M > 0 && (Q = tdiv(N,F)) ) {
1235: N = Q;
1236: M--;
1237: }
1238: if ( M ) R = cons([F,M],R);
1239: }
1240: D = reverse(R);
1241: return [N,D];
1242: }
1243:
1244:
1245: def comppfrd(P,R)
1246: {
1247: P = qsort(P,n_wishartd.compt); P=reverse(P);
1248: R = qsort(R,n_wishartd.compt); R=reverse(R);
1249: if ( length(P) != length(R) ) return 0;
1250: for ( ; P != []; P = cdr(P), R = cdr(R) ) {
1251: P0 = car(P); R0 = car(R);
1252: C0 = ctord(P0[0]); C1 = R0[0];
1253: D0 = ftop(C0[1]); D1 = ftop(C1[1]);
1254: if ( (D0 == D1) && C0[0] == -C1[0] ) continue;
1255: if ( (D0 == -D1) && C0[0] == C1[0] ) continue;
1256: return 0;
1257: }
1258: return 1;
1259: }
1260:
1261: def multpf(T,F)
1262: {
1263: E = dp_etov(T[1]);
1264: N = length(E);
1265: Z = F;
1266: for ( J = 1; J < N; J++ )
1267: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
1268: Z = mulcpf(T[0],Z);
1269: return Z;
1270: }
1271:
1272: def spolypf(F,G)
1273: {
1274: F0 = car(F); G0 = car(G);
1275: DF = F0[1]; DG = G0[1];
1276: L = dp_lcm(DF,DG);
1277: F1 = multpf([F0[0],dp_subd(L,DF)],F);
1278: G1 = multpf([G0[0],dp_subd(L,DG)],G);
1279: S = subpf(F1,G1);
1280: return S;
1281: }
1282:
1283: def nfpf(F,G)
1284: {
1285: Z = F;
1286: R = [];
1287: while ( Z != [] ) {
1288: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1289: for ( T = G; T != []; T = cdr(T) ) {
1290: T0 = car(T);
1291: if ( dp_redble(D,T0[0][1]) ) break;
1292: }
1293: if ( T != [] ) {
1294: CG = ctorat(T0[0][0]);
1295: C = mulc(C,[[-1/CG,[]]]);
1296: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1297: Z = addpf(Z,S);
1298: } else {
1299: R = cons(Z0,R);
1300: Z = cdr(Z);
1301: }
1302: }
1303: S = [];
1304: for ( T = R; T != []; T = cdr(T) ) {
1305: U = ctord(T[0][0]);
1306: if ( U[0] )
1307: S = cons(T[0],S);
1308: }
1309: return S;
1310: }
1311:
1312: def nfpf0(F,G)
1313: {
1314: Z = F;
1315: R = [];
1316: while ( Z != [] ) {
1317: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1318: for ( T = G; T != []; T = cdr(T) ) {
1319: T0 = car(T);
1320: if ( dp_redble(D,T0[0][1]) ) break;
1321: }
1322: if ( T != [] ) {
1323: CG = ctorat(T0[0][0]);
1324: C = mulc(C,[[-1/CG,[]]]);
1325: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1326: Z = addpf(Z,S);
1327: } else {
1328: R = cons(Z0,R);
1329: Z = cdr(Z);
1330: }
1331: }
1332: S = [];
1333: for ( T = R; T != []; T = cdr(T) ) {
1334: U = ctord(T[0][0]);
1335: if ( U[0] )
1336: S = cons(T[0],S);
1337: }
1338: return S;
1339: }
1340:
1.3 noro 1341: def nfpf1(F,G)
1342: {
1343: Z = F;
1344: R = [];
1345: while ( Z != [] ) {
1346: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1347: for ( T = G; T != []; T = cdr(T) ) {
1348: T0 = car(T);
1349: if ( dp_redble(D,T0[0][1]) ) break;
1350: }
1351: if ( T != [] ) {
1352: CG = ctorat(T0[0][0]);
1353: C = mulc(C,[[-1/CG,[]]]);
1354: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1355: Z = addpf(Z,S);
1356: } else {
1357: R = cons(Z0,R);
1358: Z = cdr(Z);
1359: }
1360: }
1361: return reverse(R);
1362: }
1363:
1364: def nfpf_zero(F,G)
1365: {
1366: Z = F;
1367: R = [];
1368: while ( Z != [] ) {
1369: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1370: U = ctord(C);
1371: if ( !U[0] ) { Z = cdr(Z); continue;}
1372: for ( T = G; T != []; T = cdr(T) ) {
1373: T0 = car(T);
1374: if ( dp_redble(D,T0[0][1]) ) break;
1375: }
1376: if ( T != [] ) {
1377: CG = ctorat(T0[0][0]);
1378: C = mulc(C,[[-1/CG,[]]]);
1379: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1380: Z = addpf(Z,S);
1381: } else {
1382: R = cons(Z0,R);
1383: Z = cdr(Z);
1384: }
1385: }
1386: S = [];
1387: for ( T = R; T != []; T = cdr(T) ) {
1388: U = ctord(T[0][0]);
1389: if ( U[0] )
1390: S = cons(T[0],S);
1391: }
1392: print("");
1393: return S;
1394: }
1395:
1.1 noro 1396: def pftoeuler(F)
1397: {
1398: VDV = [y1,dy1];
1399: P = pftop(F);
1400: Z = dp_ptod(P,VDV);
1401: E = dp_etov(dp_ht(Z));
1402: S = E[1]-E[0];
1403: if ( S < 0 )
1404: error("pftoeuler : invalid input");
1405: P *= y1^S;
1406: Z = dp_ptod(P,VDV);
1407: E = dp_etov(dp_ht(Z));
1408: D = E[0];
1409: C = vector(D+1);
1410: C[0] = 1;
1411: for ( I = 1; I <= D; I++ )
1412: C[I] = C[I-1]*(t-I+1);
1413: R = 0;
1414: for ( T = Z; T; T = dp_rest(T) ) {
1415: E = dp_etov(dp_ht(T));
1416: S = E[0]-E[1];
1417: if ( S < 0 )
1418: error("pftoeuler : invalid input");
1419: R += dp_hc(T)*y1^S*C[E[1]];
1420: }
1421: return R;
1422: }
1423:
1424: def gammam(M,A)
1425: {
1426: R = 1;
1427: for ( I = 1; I <= M; I++ )
1428: R *= mpfr_gamma(A-(I-1)/2);
1429: R *= eval(@pi^(1/4*M*(M-1)));
1430: return R;
1431: }
1432:
1.3 noro 1433: def genprc(M,N,PL,SL)
1.1 noro 1434: {
1435: A = (M+1)/2; C = (N+M+1)/2;
1436: DetS = 1;
1437: TrIS = 0;
1438: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1439: DetS *= car(S)^car(T);
1440: TrIS += car(T)/car(S);
1441: }
1442: C = 1/eval(DetS^(1/2*N)*2^(1/2*N*M));
1443: return C;
1444: }
1445:
1446: /*
1447: * PL=[P1,P2,...]: sizes of blocks
1448: * SL=[S1,S2,...]: [S1xP1,S2xP2,..]
1449: */
1450:
1451: def prob_by_hgm(M,N,PL,SL,X)
1452: {
1453: A = (M+1)/2; C = (N+M+1)/2;
1454:
1455: B = []; V = []; Beta = []; SBeta = 0;
1456: Prev = 1;
1457: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1458: V = cons(strtov("y"+rtostr(Prev)),V);
1459: B = cons([Prev,Prev+car(T)-1],B);
1460: Prev += car(T);
1461: Beta = cons(1/(2*car(S)),Beta);
1462: SBeta += car(T)/(2*car(S));
1463: }
1464: if ( Prev != M+1 )
1465: error("invalid block specification");
1466: B = reverse(B); V = reverse(V);
1467: Beta = reverse(Beta);
1468:
1469: T0 = time();
1470: if ( type(Z=getopt(eq)) == -1 )
1471: Z = diagpf(M,B);
1472: T1 = time(); Tdiag = (T1[0]-T0[0])+(T1[1]-T0[1]);
1473:
1474: if ( type(Step=getopt(step)) == -1 )
1475: Step = 10000;
1476:
1477: Z = subst(Z,a,A,c,C);
1.3 noro 1478: if ( type(X0=getopt(x0)) == -1 ) {
1479: if ( type(Init=getopt(init)) == -1 ) Init = 1;
1480: X0 = 0;
1481: Len = length(Beta);
1482: for ( I = 0; I < Len; I++ )
1483: if ( Beta[I] > X0 ) X0 = Beta[I];
1484: X0 = Init/X0;
1485: }
1486: if ( type(Rk=getopt(rk)) == -1 ) Rk = 5;
1487: if ( type(Eps=getopt(eps)) == -1 ) Eps = 10^(-4);
1488: if ( type(MapleOut=getopt(mapleout)) == -1 ) MapleOut = 0;
1489: F = genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,Rk,Step,X0,Eps|mapleout=MapleOut);
1.1 noro 1490: if ( Print ) print(["diag",Tdiag,"pfaffian",Tp,"ps",Tps,"rk",Trk]);
1491: return F;
1492: }
1493:
1494: def prob_by_ps(M,N,PL,SL,X)
1495: {
1496: A = (M+1)/2; C = (N+M+1)/2;
1497:
1498: B = []; V = []; Beta = []; SBeta = 0;
1499: Prev = 1;
1500: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1501: V = cons(strtov("y"+rtostr(Prev)),V);
1502: B = cons([Prev,Prev+car(T)-1],B);
1503: Prev += car(T);
1504: Beta = cons(1/(2*car(S)),Beta);
1505: SBeta += car(T)/(2*car(S));
1506: }
1507: if ( Prev != M+1 )
1508: error("invalid block specification");
1509: B = reverse(B); V = reverse(V);
1510: Beta = reverse(Beta);
1511:
1512: if ( type(Z=getopt(eq)) == -1 )
1513: Z = diagpf(M,B);
1514:
1515: Z = subst(Z,a,A,c,C);
1516: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
1.3 noro 1517: C = GM*genprc(M,N,PL,SL);
1.1 noro 1518:
1519: Beta = ltov(Beta)*eval(exp(0));
1520: X *= eval(exp(0));
1521: L = psvalue(Z,V,Beta*X); PS = L[0];
1522: MN2 = M*N/2;
1523: ExpF = eval(X^(MN2)*exp(-SBeta*X));
1524: return C*L[1]*ExpF;
1525: }
1526:
1527: def ps(Z,V,TD)
1528: {
1529: Tact = Tgr = Tactmul = 0;
1530: G = map(ptozp,map(pftop,Z));
1531: M = length(V);
1532: N = length(G);
1533: for ( I = 0, DV = []; I < M; I++ )
1534: DV = cons(strtov("d"+rtostr(V[I])),DV);
1535: DV = reverse(DV);
1536: VDV = append(V,DV);
1537: DG = map(decomp_op,G,V,DV);
1538: F = [1];
1539: R = 1;
1540: Den = 1;
1541: for ( I = 1 ; I <= TD; I++ ) {
1542: L = create_homo(V,I);
1543: FI = L[0]; CV = L[1];
1544: Eq = [];
1545: for ( K = 0; K < N; K++ ) {
1546: P = DG[K]; Len = length(P);
1547: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
1548: T0=time();
1549: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
1550: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
1551: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
1552: T1=time(); Tact += T1[0]-T0[0];
1553: if ( Lhs ) {
1554: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
1555: Eq = cons(dp_hc(T),Eq);
1556: }
1557: }
1558: }
1559: T0=time();
1560: #if 0
1561: Sol = solve_leq(Eq,CV);
1562: #else
1563: Sol = nd_f4(Eq,CV,0,0);
1564: #endif
1565: L = p_true_nf(FI,Sol,CV,0);
1566: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
1567: FI = L[0]*(Den1/L[1]);
1568: T1=time(); Tgr += T1[0]-T0[0];
1569:
1570: MJ = Den1/Den; Den = Den1;
1571: for ( S = [], T = F; T != []; T = cdr(T) )
1572: S = cons(MJ*car(T),S);
1573: F = cons(FI,reverse(S));
1574: R = R*MJ+car(F);
1575: }
1576: return R/Den;
1577: }
1578:
1579: def msubst(F,V,X)
1580: {
1581: M = length(V);
1582: for ( J = 0, F0 = F*eval(exp(0)); J < M; J++ )
1583: F0 = subst(F0,V[J],X[J]);
1584: return F0;
1585: }
1586:
1587: def psvalue(Z,V,X)
1588: {
1589: Tact = Tgr = Tactmul = 0;
1590: G = map(ptozp,map(pftop,Z));
1591: M = length(V);
1592: N = length(G);
1593: for ( I = 0, DV = []; I < M; I++ )
1594: DV = cons(strtov("d"+rtostr(V[I])),DV);
1595: DV = reverse(DV);
1596: VDV = append(V,DV);
1597: DG = map(decomp_op,G,V,DV);
1598: Prev = getopt(prev);
1599: if ( type(Prev) == -1 ) {
1600: F = [1];
1601: R = 1;
1602: Val = eval(exp(0));
1603: Den = 1;
1604: I = 1;
1605: } else {
1606: /* Prev = [R/Den,Val1,F,Den] */
1607: Den = Prev[3];
1608: F = Prev[2];
1609: R = Prev[0]*Den;
1610: I = length(F);
1611: Val = msubst(Prev[0],V,X);
1612: ValI = msubst(F[0],V,X)/Den;
1613: if ( Val-ValI == Val )
1614: return [Prev[0],Val,F,Den];
1615: }
1616: for ( ; ; I++ ) {
1617: L = create_homo(V,I);
1618: FI = L[0]; CV = L[1];
1619: Eq = [];
1620: for ( K = 0; K < N; K++ ) {
1621: P = DG[K]; Len = length(P);
1622: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
1623: T0=time();
1624: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
1625: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
1626: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
1627: T1=time(); Tact += T1[0]-T0[0];
1628: if ( Lhs ) {
1629: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
1630: Eq = cons(dp_hc(T),Eq);
1631: }
1632: }
1633: }
1634: T0=time();
1635: Sol = solve_leq(Eq,CV);
1636: L = p_true_nf(FI,Sol,CV,0);
1.3 noro 1637: delete_uc();
1.1 noro 1638: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
1639: FI = L[0]*(Den1/L[1]);
1640: T1=time(); Tgr += T1[0]-T0[0];
1641:
1642: MJ = Den1/Den; Den = Den1;
1643: for ( S = [], T = F; T != []; T = cdr(T) )
1644: S = cons(MJ*car(T),S);
1645: F = cons(FI,reverse(S));
1646: R = R*MJ+car(F);
1647: F0 = msubst(FI,V,X)/Den;
1648: Val1 = Val+F0;
1649: if ( Val1 == Val ) {
1650: if ( Print ) print(["I=",I,"act",Tact,"actmul",Tactmul,"gr",Tgr]);
1651: return [R/Den,Val1,F,Den];
1652: } else {
1653: if ( Print ) print([F0],2);
1654: Val = Val1;
1655: }
1656: }
1657: }
1658:
1659: /* P -> [P0,P1,...] Pi = y^a dy^b, |a|-|b|=i */
1660:
1661: def decomp_op(P,V,DV)
1662: {
1663: M = length(V);
1664: VDV = append(V,DV);
1665: D = dp_ptod(P,VDV);
1666: for ( I = 0, T = D; T; T = dp_rest(T), I++ );
1667: for ( T = D, NotSet = 1; T; T = dp_rest(T) ) {
1668: E = dp_etov(T);
1669: for ( K = 0, J = 0; J < M; J++ )
1670: K += E[J]-E[M+J];
1671: if ( NotSet ) {
1672: Max = Min = K; NotSet = 0;
1673: } else if ( K > Max ) Max = K;
1674: else if ( K < Min ) Min = K;
1675: }
1676: W = vector(Max-Min+1);
1677: for ( T = D; T; T = dp_rest(T) ) {
1678: E = dp_etov(T);
1679: for ( K = 0, J = 0; J < M; J++ )
1680: K += E[J]-E[M+J];
1681: W[K-Min] += dp_hm(T);
1682: }
1683: W = map(dp_ptod,map(dp_dtop,W,VDV),DV);
1684: return W;
1685: }
1686:
1687: def create_homo(V,TD)
1688: {
1689: M = length(V);
1690: for ( S = 0, I = 0; I < M; I++ ) S += V[I];
1691: Tmp = 0;
1692: Nc = 0;
1693: for ( RI = 0, D = dp_ptod(S^TD,V); D; D = dp_rest(D), Nc++ ) {
1694: E = dp_etov(D);
1695: T = uc();
1696: U = dp_dtop(dp_ht(D),V);
1697: RI += T*U;
1698: Tmp += T;
1699: }
1700: CV = vector(Nc);
1701: for ( I = 0; I < Nc; I++ ) {
1702: CV[I] = var(Tmp); Tmp -= CV[I];
1703: }
1704: return [RI,vtol(CV)];
1705: }
1706:
1707: def act_op(P,V,F)
1708: {
1709: N = length(V);
1710: for ( R = 0; P; P = dp_rest(P) ) {
1711: C = dp_hc(P); E = dp_etov(P);
1712: for ( I = 0, S = F; I < N; I++ )
1713: for ( J = 0; J < E[I]; J++ ) S = diff(S,V[I]);
1714: T0 = time();
1715: R += C*S;
1716: T1 = time(); Tactmul += T1[0]-T0[0];
1717: }
1718: return R;
1719: }
1720:
1721: def gen_basis(D,DV)
1722: {
1723: N = length(D);
1724: B = [];
1725: E = vector(N);
1726: while ( 1 ) {
1727: B = cons(dp_dtop(dp_vtoe(E),DV),B);
1728: for ( I = N-1; I >= 0; I-- )
1729: if ( E[I] < D[I]-1 ) break;
1730: if ( I < 0 ) return reverse(B);
1731: E[I]++;
1732: for ( J = I+1; J < N; J++ ) E[J] = 0;
1733: }
1734: }
1735:
1736: def pfaffian(Z)
1737: {
1738: N = length(Z);
1739: D = vector(N);
1740: Mat = vector(N);
1741: V = vars(Z);
1742: DV = vector(N);
1743: for ( I = 0; I < N; I++ )
1744: DV[I] = strtov("d"+rtostr(V[I]));
1745: DV = vtol(DV);
1746: for ( I = 0; I < N; I++ ) {
1747: ZI = Z[I];
1748: HI = ZI[0][1];
1749: DI = dp_dtop(HI,PF_DV);
1750: for ( J = 0; J < N; J++ )
1751: if ( var(DI) == DV[J] )
1752: break;
1753: D[J] = deg(DI,DV[J]);
1754: }
1755: B = gen_basis(D,DV);
1756: Dim = length(B);
1757: Hist = [];
1758: for ( I = 0; I < N; I++ ) {
1759: if ( Print ) print(["I=",I]);
1760: A = matrix(Dim,Dim);
1761: for ( K = 0; K < Dim; K++ )
1762: for ( L = 0; L < Dim; L++ )
1763: A[K][L] = [];
1764: for ( J = 0; J < Dim; J++ ) {
1765: if ( Print ) print(["J=",J],2);
1766: Mon = DV[I]*B[J];
1767: for ( T = Hist; T != []; T = cdr(T) )
1768: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
1769: if ( (T != []) ) {
1770: for ( L = 0; L < N; L++ )
1771: if ( Q == DV[L] ) break;
1772: Red1 = muldpf(V[L],car(T)[1]);
1773: Red = nfpf0(Red1,Z);
1774: } else
1775: Red = nfpf0([mtot(Mon,1)],Z);
1776: Hist = cons([Mon,Red],Hist);
1777: for ( T = Red; T != []; T = cdr(T) ) {
1778: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
1779: for ( K = 0; K < Dim; K++ )
1780: if ( B[K] == Op ) break;
1781: if ( K == Dim )
1782: error("afo");
1783: else
1784: A[J][K] = C;
1785: }
1786: }
1787: if ( Print ) print("");
1788: A = map(ctord,A);
1789: Lcm = [];
1790: for ( K = 0; K < Dim; K++ )
1791: for ( L = 0; L < Dim; L++ )
1792: Lcm = lcmf(Lcm,A[K][L][1]);
1793: for ( K = 0; K < Dim; K++ )
1794: for ( L = 0; L < Dim; L++ ) {
1795: Num = A[K][L][0]; Den = A[K][L][1];
1796: A[K][L] = Num*ftorat(divf(Lcm,Den));
1797: }
1798: Lcm = ftorat(Lcm);
1799: if ( !Lcm ) Lcm = 1;
1800: Mat[I] = [A,Lcm];
1801: }
1802: return [Mat,B];
1803: }
1804:
1.3 noro 1805: def pfaffian2(Z,V,Beta,SBeta,MN2,XV)
1.1 noro 1806: {
1807: N = length(Z);
1808: D = vector(N);
1809: Mat = vector(N);
1810: DV = vector(N);
1811: for ( I = 0; I < N; I++ )
1812: DV[I] = strtov("d"+rtostr(V[I]));
1813: DV = vtol(DV);
1814: for ( I = 0; I < N; I++ ) {
1815: ZI = Z[I];
1816: HI = ZI[0][1];
1817: DI = dp_dtop(HI,PF_DV);
1818: for ( J = 0; J < N; J++ )
1819: if ( var(DI) == DV[J] )
1820: break;
1821: D[J] = deg(DI,DV[J]);
1822: }
1823: B = gen_basis(D,DV);
1824: Dim = length(B);
1825: Hist = [];
1826: for ( I = 0; I < N; I++ ) {
1827: if ( Print ) print(["I=",I]);
1828: A = matrix(Dim,Dim);
1829: for ( K = 0; K < Dim; K++ )
1830: for ( L = 0; L < Dim; L++ )
1831: A[K][L] = [];
1832: for ( J = 0; J < Dim; J++ ) {
1833: if ( Print ) print(["J=",J],2);
1834: Mon = DV[I]*B[J];
1835: for ( T = Hist; T != []; T = cdr(T) )
1836: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
1837: if ( (T != []) ) {
1838: for ( L = 0; L < N; L++ )
1839: if ( Q == DV[L] ) break;
1840: Red1 = muldpf(V[L],car(T)[1]);
1.3 noro 1841: Red = nfpf1(Red1,Z);
1.1 noro 1842: } else
1.3 noro 1843: Red = nfpf1([mtot(Mon,1)],Z);
1.1 noro 1844: Hist = cons([Mon,Red],Hist);
1845: for ( T = Red; T != []; T = cdr(T) ) {
1846: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
1847: for ( K = 0; K < Dim; K++ )
1848: if ( B[K] == Op ) break;
1849: if ( K == Dim )
1850: error("afo");
1851: else
1852: A[J][K] = C;
1853: }
1854: }
1.3 noro 1855: for ( K = 0; K < N; K++ )
1856: A = subst(A,V[K],Beta[K]*XV);
1.1 noro 1857: A = map(ctord,A);
1858: Lcm = [];
1859: for ( K = 0; K < Dim; K++ )
1860: for ( L = 0; L < Dim; L++ )
1861: Lcm = lcmf(Lcm,A[K][L][1]);
1862: for ( K = 0; K < Dim; K++ )
1863: for ( L = 0; L < Dim; L++ ) {
1864: Num = A[K][L][0]; Den = A[K][L][1];
1865: A[K][L] = Num*ftorat(divf(Lcm,Den));
1866: }
1867: Lcm = ftorat(Lcm);
1868: if ( !Lcm ) Lcm = 1;
1869: A = map(red,A/Lcm);
1870: Lcm = 1;
1871: for ( K = 0; K < Dim; K++ )
1872: for ( L = 0; L < Dim; L++ )
1873: Lcm = lcm(dn(A[K][L]),Lcm);
1874: for ( K = 0; K < Dim; K++ )
1875: for ( L = 0; L < Dim; L++ )
1876: A[K][L] = nm(A[K][L])*sdiv(Lcm,dn(A[K][L]));
1877: Mat[I] = [A,Lcm];
1878: }
1.3 noro 1879: Lcm = XV;
1.1 noro 1880: for ( I = 0; I < N; I++ )
1881: Lcm = lcm(Mat[I][1],Lcm);
1.3 noro 1882: /* R = (MN2/XV-SBeta)*Id = (MN2-SBeta*XV)*(Lcm/XV)/Lcm*Id */
1.1 noro 1883: R = matrix(Dim,Dim);
1.3 noro 1884: D = (MN2-SBeta*XV)*sdiv(Lcm,XV);
1885: for ( I = 0; I < Dim; I++ ) R[I][I] = D;
1.1 noro 1886: for ( I = 0; I < N; I++ ) {
1887: A = Mat[I][0];
1888: for ( K = 0; K < Dim; K++ )
1.3 noro 1889: for ( L = 0; L < Dim; L++ ) {
1890: R[K][L] += Beta[I]*(A[K][L])*sdiv(Lcm,Mat[I][1]);
1891: }
1.1 noro 1892: }
1893: Deg = 0;
1894: for ( K = 0; K < Dim; K++ )
1895: for ( L = 0; L < Dim; L++ ) {
1896: DegT = deg(R[K][L],XV);
1897: if ( DegT > Deg ) Deg = DegT;
1898: }
1899: RA = vector(Deg+1);
1900: for ( I = 0; I <= Deg; I++ )
1901: RA[I] = map(coef,R,I);
1902: return [[RA,Lcm],B];
1903: }
1904:
1905: def evalpf(F,V,X)
1906: {
1907: if ( !F ) return 0;
1908: F0 = F;
1909: N = length(V);
1910: for ( I = 0; I < N; I++ )
1911: F0 = subst(F0,V[I],X[I]);
1912: return F0;
1913: }
1914:
1915: def eval_pfaffian(Mat,Beta,SBeta,MN2,V,XI,F)
1916: {
1917: N = length(V);
1918: Dim = size(Mat[0][0])[0];
1919: R = vector(Dim);
1920: Mul = vector(N);
1921: X = XI*Beta;
1922: X = vtol(X);
1923: for ( K = 0; K < N; K++ )
1924: Mul[K] = Beta[K]/evalpf(Mat[K][1],V,X);
1925: for ( K = 0; K < N; K++ ) {
1926: MatK = Mat[K][0];
1927: for ( I = 0; I < Dim; I++ ) {
1928: for ( U = J = 0; J < Dim; J++ )
1929: U += substr2np(MatK[I][J],V,X)*F[J];
1930: R[I] += Mul[K]*U;
1931: }
1932: }
1933: for ( I = 0; I < Dim; I++ ) R[I] -= (SBeta-MN2/XI)*F[I];
1934: return R;
1935: }
1936:
1.3 noro 1937: def eval_pfaffian2(Mat,XI,F)
1.1 noro 1938: {
1.3 noro 1939: MA = Mat[0];
1940: Den = Mat[1];
1941: if ( T = var(Den) ) Den = subst(Den,T,XI);
1.1 noro 1942: Len = length(MA);
1943: Dim = size(MA[0])[0];
1944: R = vector(Dim);
1.3 noro 1945: for ( I = Len-1; I >= 0; I-- ) {
1946: R = XI*R+MA[I]*F;
1.1 noro 1947: }
1.3 noro 1948: R = R/Den;
1.1 noro 1949: return R;
1950: }
1951:
1.3 noro 1952: def genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,RK,Step,X0,Eps)
1.1 noro 1953: {
1.3 noro 1954: if ( type(MapleOut=getopt(mapleout)) == -1 ) MapleOut = 0;
1955: ctrl("bigfloat",1);
1956: O = eval(exp(0));
1.1 noro 1957: Len = length(V);
1958: DV = vector(Len);
1959: for ( I = 0; I < Len; I++ )
1960: DV[I] = strtov("d"+rtostr(V[I]));
1961: DV = vtol(DV);
1.3 noro 1962: A = (1+M)/2; C = (1+M+N)/2; MN2 = M*N/2;
1.1 noro 1963: Z0 = subst(Z,a,A,c,C);
1964: T0 = time();
1.3 noro 1965: Q = pfaffian2(Z0,V,Beta,SBeta,MN2,t); Mat = Q[0]; Base = Q[1];
1.1 noro 1966: MatLen = length(Q[0][0]);
1967: for ( I = 0; I < MatLen; I++ )
1.3 noro 1968: Mat[0][I] *= O;
1.1 noro 1969: T1 = time(); Tp = (T1[0]-T0[0])+(T1[1]-T0[1]);
1970: T0 = time();
1.3 noro 1971: Beta = ltov(Beta)*O;
1972: X *= O;
1973: Beta0 = Beta*X0;
1.1 noro 1974: L = psvalue(Z0,V,Beta0); PS = L[0];
1.3 noro 1975: if ( Print ) print(["x0=",X0]);
1.1 noro 1976: T1 = time(); Tps = (T1[0]-T0[0])+(T1[1]-T0[1]);
1977: T0 = time();
1978: Dim = length(Base);
1979: ExpF = eval(X0^(MN2)*exp(-SBeta*X0));
1.3 noro 1980: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2)*genprc(M,N,PL,SL);
1981: ExpF *= GM;
1982: F0 = vector(Dim);
1.1 noro 1983: for ( I = 0; I < Dim; I++ ) {
1984: DPS = act_op(dp_ptod(Base[I],DV),V,PS);
1985: for ( J = 0; J < Len; J++ )
1986: DPS = subst(DPS,V[J],Beta0[J]);
1.3 noro 1987: F0[I] = DPS*ExpF;
1988: }
1989: First = 1;
1990: Val = [];
1991: LVal = [];
1992: DVal = [];
1993: RDVal = [];
1994: if ( MapleOut ) {
1995: mapleout(Mat[0],Mat[1],X0,F0,X,MapleOut);
1996: return;
1.1 noro 1997: }
1.3 noro 1998: for ( I = 0; ; Step *= 2, I++ ) {
1999: if ( Print ) print("Step=",0); if ( Print ) print(Step);
1.1 noro 2000: F = F0*1;
1.3 noro 2001: T = F[0]; F /= T;
2002: R = [[X0,T]];
2003: R1 = rk_ratmat(RK,Mat[0],Mat[1],X0,X,Step,F);
2004: R = append(R1,R);
2005:
2006: setround(1);
2007: Adj = eval(exp(0));
2008: for ( T = R; T != []; T = cdr(T) )
2009: Adj *= car(T)[1];
2010: LVal = cons(Adj,LVal);
2011: setround(2);
2012: Adj = eval(exp(0));
2013: for ( T = R; T != []; T = cdr(T) )
2014: Adj *= car(T)[1];
2015: Val = cons(Adj,Val);
2016: setround(0);
2017:
2018: if ( Print ) { print(""); print([LVal[0],Val[0]]); }
2019: if ( I >= 1 ) DVal = cons(Val[0]-Val[1],DVal);
2020: if ( I >= 2 ) {
2021: RDVal = cons(DVal[1]/DVal[0],RDVal);
2022: print(["log ratio=",eval(log(RDVal[0])/log(2))]);
2023: }
2024: if ( I >= 1 && abs(1-Val[0]/Val[1])<Eps ) {
2025: T1 = time(); Trk = (T1[0]-T0[0])+(T1[1]-T0[1]);
2026: return Val[0];
1.1 noro 2027: }
2028: if ( Print ) print("");
2029: }
2030: }
2031:
2032: def solve_leq(L,V)
2033: {
1.3 noro 2034: #if 0
1.1 noro 2035: N = length(V);
2036: B = [];
2037: for ( I = 0; I < N; I++, L = cdr(L) ) B = cons(car(L),B);
2038: while ( 1 ) {
2039: Sol = nd_f4(B,V,0,0);
2040: if ( length(Sol) != N ) {
2041: B = cons(car(L),Sol); L = cdr(L);
2042: } else {
2043: return Sol;
2044: }
2045: }
1.3 noro 2046: #else
2047: Sol = nd_f4_trace(L,V,0,-1,0);
2048: return Sol;
2049: #endif
1.1 noro 2050: }
2051:
2052:
2053: def ptom(P,V)
2054: {
2055: M = length(P); N = length(V);
2056: A = matrix(M,N+1);
2057: for ( I = 0; I < M; I++ ) {
2058: F = ptozp(P[I]);
2059: VT = V;
2060: J = 0;
2061: while ( F )
2062: if ( type(F) <= 1 ) {
2063: A[I][N] = F; break;
2064: } else {
2065: for ( FV = var(F); FV != car(VT); VT = cdr(VT), J++ );
2066: A[I][J] = coef(F,1);
2067: F -= A[I][J]*FV;
2068: }
2069: }
2070: return A;
2071: }
2072:
2073: def vton(V1,V)
2074: {
2075: N = length(V);
2076: for ( I = 0; I < N; I++ )
2077: if ( V1 == V[I] ) return I;
2078: error("vton : no such variable");
2079: }
2080:
2081: def ftotex(F,V)
2082: {
2083: if ( F == [] ) return "1";
2084: R = "";
2085: for ( T = F; T != []; T = cdr(T) ) {
2086: Fac = car(T)[0]; M = car(T)[1];
2087: V1 = var(Fac); NV1 = vton(V1,V);
2088: if ( Fac == V1 ) SFac = "y_{"+rtostr(NV1)+"}";
2089: else {
2090: V2 = var(Fac-V1);
2091: NV2 = vton(V2,V);
2092: SFac = "(y_{"+rtostr(NV1)+"}-y_{"+rtostr(NV2)+"})";
2093: }
2094: if ( M == 1 ) R += SFac;
2095: else R += SFac+"^{"+rtostr(M)+"}";
2096: }
2097: return R;
2098: }
2099:
2100: def ctotex(C,V)
2101: {
2102: R = "";
2103: for ( T = C; T != []; T = cdr(T) ) {
2104: if ( T != C ) R += "+";
2105: N = quotetotex(objtoquote(car(T)[0]));
2106: if ( car(T)[1] == [] ) {
2107: R += "("+N+")";
2108: } else {
2109: D = ftotex(car(T)[1],V);
2110: R += "\\frac{"+N+"}{"+D+"}";
2111: }
2112: }
2113: return R;
2114: }
2115:
2116: def optotex(Op,DV)
2117: {
2118: E = dp_etov(Op);
2119: N = length(E);
2120: R = "";
2121: for ( I = 0; I < N; I++ )
2122: if ( E[I] )
2123: if ( E[I] == 1 )
2124: R += "\\partial_{"+rtostr(I)+"}";
2125: else
2126: R += "\\partial_{"+rtostr(I)+"}^{"+rtostr(E[I])+"}";
2127: return R;
2128: }
2129:
2130: def pftotex(P,V,DV)
2131: {
2132: R = "";
2133: for ( T = P; T != []; T = cdr(T) ) {
2134: if ( T != P ) R += "+";
2135: C = ctotex(car(T)[0],V); Op = optotex(car(T)[1],DV);
2136: R += "("+C+")"+Op+"\\\\";
2137: }
2138: return R;
2139: }
2140:
2141: def eqtotex(Z)
2142: {
2143: shell("rm -f __tmp__.tex");
2144: output("__tmp__.tex");
2145: N = length(Z);
2146: print("\\documentclass[12pt]{jsarticle}");
2147: print("\\begin{document}");
2148: print("\\large\\noindent");
2149: for ( I = 0; I < N; I++ ) {
2150: print("$h_"+rtostr(I)+"=",0);
2151: T = mulcpf([[-1,[]]],Z[I]);
2152: print(pftotex(T,PF_V,PF_DV),0);
2153: print("$\\\\");
2154: }
2155: print("\\end{document}");
2156: output();
2157: shell("latex __tmp__");
2158: shell("xdvi __tmp__");
2159: }
2160:
2161: /*
2162: * for a partition L=[N1,...,Nl]
2163: * compute the coefficient of x1^N1*...*xl^Nl in phi((x1+...+xM)^N)
2164: * where phi(x(s(1))^n1*...*x(s(k))^nk)=x1^n1*...*xk^nk (n1>=...>=nk)
2165: * if count_dup(L)=[I1,...,Is], then the coefficient is
2166: * C=N!/(N1!*...*Nl!)*M!/((M-l)!*I1!*...*Is!)
2167: * return C*<<0,...0,N1,...,Nl,0,...,0>>
2168: position 0 Base
2169: */
2170:
2171: def coef_partition(L,N,M,Base)
2172: {
2173: R = fac(N);
2174: for ( T = L; T != []; T = cdr(T) )
2175: R = sdiv(R,fac(car(T)));
2176: S = count_dup(L);
2177: R *= fac(M)/fac(M-length(L));
2178: for ( T = S; T != []; T = cdr(T) )
2179: R = sdiv(R,fac(car(T)));
2180: E = vector(length(PF_DV));
2181: for ( I = Base, T = L; T != []; T = cdr(T), I++ )
2182: E[I] = car(T);
2183: return R*dp_vtoe(E);
2184: }
2185:
2186: /*
2187: * for a partition L=[N1,...,Nk], return [I1,...,Is]
2188: * s.t. N1=...=N(I1)>N(I1+1)=...=N(I1+I2)>...
2189: */
2190:
2191: def count_dup(L)
2192: {
2193: D = [];
2194: T = L;
2195: while ( T != [] ) {
2196: T0 = car(T); T = cdr(T);
2197: for ( I = 1; T != [] && car(T) == T0; T = cdr(T) ) I++;
2198: D = cons(I,D);
2199: }
2200: return D;
2201: }
2202:
2203: /* generate (N_1,...,N_L) s.t. N_1+...+N_L=N, N_1>=...>=N_L>=1, L<= K */
2204:
2205: def all_partition(N,K)
2206: {
2207: R = [];
2208: for ( I = 1; I <= K; I++ )
2209: R = append(partition(N,I),R);
2210: return map(reverse,R);
2211: }
2212:
2213: /* generate (N_1,...,N_K) s.t. N_1+...+N_K=N, 1<=N_1<=...<=N_K */
2214:
2215: def partition(N,K)
2216: {
2217: if ( N < K ) return [];
2218: else if ( N == K ) {
2219: S = [];
2220: for ( I = 0; I < K; I++ )
2221: S = cons(1,S);
2222: return [S];
2223: } else if ( K == 0 )
2224: return [];
2225: else {
2226: R = [];
2227: L = partition(N-K,K);
2228: for ( T = L; T != []; T = cdr(T) ) {
2229: S = [];
2230: for ( U = car(T); U != []; U = cdr(U) )
2231: S = cons(car(U)+1,S);
2232: R = cons(reverse(S),R);
2233: }
2234: L = partition(N-1,K-1);
2235: for ( T = L; T != []; T = cdr(T) )
2236: R = cons(cons(1,car(T)),R);
2237: return R;
2238: }
2239: }
2240:
2241: def mul_lopitalpf2(P,I,N,K,I0,I1,E)
2242: {
2243: YI = PF_V[I]; DYI = PF_DV[I];
2244: R = [mtot(DYI^2,1)];
2245: for ( J = I0; J <= I1; J++ ) {
2246: if ( J == I ) continue;
2247: YJ = PF_V[J]; DYJ = PF_DV[J];
2248: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
2249: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
2250: }
2251: S = subpf(P[I],R);
2252: R = loph(E,I,I0,I1);
2253: if ( type(getopt(high)) == -1 )
2254: S = mul_lopitalpf(S,N,K,I0,I1,E)[1];
2255: else
2256: S = 0;
2257: return [R,S];
2258: }
2259:
2260: /* the highest part of lim_{Y(I+1),...,Y(I+N-1)->YI} E(PI) */
2261: /* where E1 = E+(0,...,2,...,0) */
2262:
2263: def loph(E,I,I0,I1)
2264: {
2265: E1 = E*1;
2266: E1[I] += 2;
2267: R = [[[[1,[]]],dp_vtoe(E1)]];
2268: S = lopn(E,I,I0,I1);
2269: S = mulcpf([[1/2,[]]],S);
2270: R = addpf(R,S);
2271: R = simppf(R,I0,I1-I0+1);
2272: return R;
2273: }
2274:
2275: /* lim_{Y(I+1),...,Y(I+N-1)->YI} dy^E(1/(YI-Y(I+1))+...+1/(YI-Y(I+N-1))) */
2276:
2277: def lopn(E,I,I0,I1)
2278: {
2279: R = [];
2280: for ( K = I0; K <= I1; K++ ) {
2281: if ( K != I ) {
2282: L = lop1(E,I,K);
2283: R = addpf(R,L);
2284: }
2285: }
2286: return R;
2287: }
2288:
2289: /* lim_{YJ->YI} dy^E(1/(YI-YJ) */
2290: def lop1(E,I,J)
2291: {
2292: YI = PF_V[I]; YJ = PF_V[J];
2293: DYI = PF_DV[I]; DYJ = PF_DV[J];
2294: R = [];
2295: R = addpf([mtot(DYI,YI-YJ)],R);
2296: R = addpf([mtot(-DYJ,YI-YJ)],R);
2297: N = length(PF_V);
2298: BI = E[I]; BJ = E[J];
2299: for ( K = 1, D = 1; K < N; K++ )
2300: if ( K != I && K != J )
2301: D *= PF_DV[K]^E[K];
2302: S = [mtot(-1/(BJ+1)*DYI^(BI+1)*DYJ^(BJ+1)*D,1)];
2303: for ( K = 0; K <= BI; K++ )
2304: if ( -BI+BJ+2*K+2 )
2305: S = addpf([mtot(fac(BI)*fac(BJ)*(-BI+BJ+2*K+2)/fac(BI-K)/fac(BJ+K+2)*DYI^(BI-K)*DYJ^(BJ+K+2)*D,1)],S);
2306: return S;
2307: }
2308:
2309: def getchi(N)
2310: {
2311: CHI=[
2312: [5 ,11.07049769,1.145476226],
2313: [10 ,18.30703805,3.940299136],
2314: [20 ,31.41043284,10.85081139],
2315: [50 ,67.50480655,34.76425168],
2316: [100,124.3421134,77.92946517],
2317: [500,553.1268089,449.1467564],
2318: [1000, 1074.679449,927.594363]];
2319: for ( T = CHI; T != []; T = cdr(T) )
2320: if ( car(T)[0] == N ) return car(T);
2321: return [];
2322: }
2323:
2324: def message(S) {
2325: dp_gr_print(S);
2326: Print = S;
2327: }
1.3 noro 2328:
2329: def mapleout(Nm,Dn,X0,F0,X,OutF)
2330: {
2331: N = length(F0);
2332: V = vector(N);
2333: VT = vector(N);
2334: M = length(Nm);
2335: R = matrix(N,N);
2336: for ( I = 0; I < M; I++ )
2337: R += Nm[I]*t^I;
2338: for ( I = 0; I < N; I++ ) {
2339: V[I] = strtov("x"+rtostr(I));
2340: VT[I] = strtov(rtostr(V[I])+"(t)");
2341: }
2342: R *= VT;
2343: output(OutF);
2344: print("dsys:={");
2345: for ( I = 0; I < N; I++ ) {
2346: print("diff("+rtostr(V[I])+"(t),t)=("+rtostr(R[I])+")/("+rtostr(Dn)+"),");
2347: }
2348: for ( I = 0; I < N; I++ ) {
2349: print(rtostr(V[I])+"("+rtostr(X0)+")="+rtostr(F0[I]),0);
2350: if ( I < N-1 ) print(",");
2351: }
2352: print("};");
2353: print("dsol:=dsolve(dsys,numeric,output=Array(["+rtostr(X)+"]));");
2354: output();
2355: }
2356:
2357: def gbcheck(Z)
2358: {
2359: N = length(Z);
2360: for ( I = 0; I < N; I++ )
2361: for ( J = I+1; J < N; J++ ) {
2362: S = spolypf(Z[I],Z[J]);
2363: R = nfpf_zero(S,Z);
2364: if ( R != [] ) return 0;
2365: else print([I,J],2);
2366: }
2367: print("");
2368: return 1;
2369: }
2370:
2371: def all_partition(N,K)
2372: {
2373: R = [];
2374: for ( I = 1; I <= K; I++ )
2375: R = append(partition(N,I),R);
2376: return map(reverse,R);
2377: }
2378:
2379: def partition(N,K)
2380: {
2381: if ( N < K ) return [];
2382: else if ( N == K ) {
2383: S = [];
2384: for ( I = 0; I < K; I++ )
2385: S = cons(1,S);
2386: return [S];
2387: } else if ( K == 0 )
2388: return [];
2389: else {
2390: R = [];
2391: L = partition(N-K,K);
2392: for ( T = L; T != []; T = cdr(T) ) {
2393: S = [];
2394: for ( U = car(T); U != []; U = cdr(U) )
2395: S = cons(car(U)+1,S);
2396: R = cons(reverse(S),R);
2397: }
2398: L = partition(N-1,K-1);
2399: for ( T = L; T != []; T = cdr(T) )
2400: R = cons(cons(1,car(T)),R);
2401: return R;
2402: }
2403: }
2404:
2405: def partition_to_pattern(L)
2406: {
2407: Cur = 1;
2408: R = [];
2409: for ( T = L; T != []; T = cdr(T) ) {
2410: R = cons([Cur,Cur+car(T)-1],R);
2411: Cur += car(T);
2412: }
2413: return reverse(R);
2414: }
2415:
2416: localf ctond$
2417:
2418: def ctond(C)
2419: {
2420: D = [];
2421: for ( T = C; T != []; T = cdr(T) )
2422: D = lcmf(D,car(T)[1]);
2423: N = [];
2424: for ( T = C; T != []; T = cdr(T) ) {
2425: D0 = divf(D,car(T)[1]);
2426: N = cons([car(T)[0],D0],N);
2427: }
2428: return [reverse(N),D];
2429: }
2430:
2431: localf tdegf$
2432:
2433: def tdegf(F)
2434: {
2435: D = 0;
2436: for ( L = F[1]; L != []; L = cdr(L) )
2437: D += car(L)[1];
2438: return D;
2439: }
2440:
2441: localf ntop$
2442:
2443: def ntop(N)
2444: {
2445: R = 0;
2446: for ( ; N != []; N = cdr(N) )
2447: R += car(N)[0]*ftop(car(N)[1]);
2448: return R;
2449: }
2450:
2451: #if 0
2452: def pwrtab(D)
2453: {
2454: Len = length(D);
2455: F = vector(Len);
2456: M = vector(Len);
2457: for ( I = 0; I < Len; I++ ) {
2458: F[I] = D[I][0]; M[I] = D[I][1];
2459: }
2460: P = vector(Len);
2461: for ( I = 0; I < Len; I++ ) {
2462: MI = M[I]; FI = F[I];
2463: PI = P[I] = vector(MI+1);
2464: for ( PI[0] = 1, J = 1; J <= MI; J++ ) PI[J] = PI[J-1]*FI;
2465: }
2466: }
2467:
2468: def ctord2(C)
2469: {
2470: D = [];
2471: for ( T = C; T != []; T = cdr(T) )
2472: D = lcmf(D,car(T)[1]);
2473: H = pwrtab(D);
2474: N = 0;
2475: for ( I = 0, T = reverse(C); T != []; T = cdr(T),I++ ) {
2476: print([Len-I],2);
2477: T0 = car(T); N0 = T0[0]; D0 = T0[1];
2478: L = lcmf(D,D0);
2479: /* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
2480: N = N*getpwr(divf(L,D))+N0*getpwr(divf(L,D0));
2481: D = L;
2482: }
2483: R = [];
2484: for ( T = D; T != []; T = cdr(T) ) {
2485: T0 = car(T); F = T0[0]; M = T0[1];
2486: while ( M > 0 && (Q = tdiv(N,F)) ) {
2487: N = Q;
2488: M--;
2489: }
2490: if ( M ) R = cons([F,M],R);
2491: }
2492: D = reverse(R);
2493: return [N,D];
2494: }
2495: #endif
2496:
2497: def degpf(P)
2498: {
2499: R = [];
2500: for ( T = P; T != []; T = cdr(T) ) {
2501: Term = car(T);
2502: C = Term[0]; D = Term[1];
2503: C0 = C[0];
2504: CN = C0[0]; CD = C0[1];
2505: R = cons([CD!=[]?CD[0][1]:0,D],R);
2506: }
2507: return reverse(R);
2508: }
2509:
2510: def pftodpf(F)
2511: {
2512: D0 = car(F)[1];
2513: N = length(dp_etov(D0));
2514: DY = [];
2515: for ( I = 0; I < N; I++ )
2516: DY = cons(strtov("dy"+rtostr(I)),DY);
2517: DY = reverse(DY);
2518: R = [];
2519: for ( T = mulcpf([[-1,[]]],F); T != []; T = cdr(T) ) {
2520: C = car(T)[0]; D = dp_dtop(car(T)[1],DY);
2521: R = cons([C,D],R);
2522: }
2523: return reverse(R);
2524: }
2525:
2526: def all_diag(N,Dir)
2527: {
2528: P = all_partition(N,N);
2529: U = map(partition_to_pattern,P);
2530: Len = length(P);
2531: for ( I = 0; I < Len; I++ ) {
2532: Z = map(pftodpf,diagpf(N,U[I]));
2533: for ( Name = rtostr(N)+"-", T = P[I]; T != []; T = cdr(T) ) {
2534: Name += rtostr(car(T));
2535: if ( length(T) > 1 ) Name += "_";
2536: }
2537: bsave(Z,Dir+"/"+Name);
2538: }
2539: }
2540:
2541: def y1y2(M)
2542: {
2543: wsetup(M);
2544: U = [[2*c-M,[[y1,1]]],[-2,[]]];
2545: for ( J = 3; J <= M; J++ )
2546: U = addc(U,[[1,[[PF_V[1]-PF_V[J],1]]]]);
2547: P = [];
2548: for ( J = 3; J <= M; J++ ) {
2549: C = addc([[1,[[PF_V[1]-PF_V[J],1]]]],[[-1,[[PF_V[1],1]]]]);
2550: P = addpf(P,[[C,dp_ptod(PF_DV[J],PF_DV)]]);
2551: }
2552: }
2553:
1.4 ! noro 2554: def help()
! 2555: {
! 2556: print("n_wishartd.diagpf(M,[[S1,E1],[S2,E2],...]) : computation of a system of PDEs satisfied by a diagonalized 1F1");
! 2557: print(" Arguments : M is the number of variables, [Si,Ei]'s are diagonal blocks s.t. S(i+1)=Ei+1.");
! 2558: print(" An example : n_wishartd.diagpf(10,[[1,9],[10,10]]) returns a system of PDEs satisfied by 1F1(y1,...,y1,y10).");
! 2559: print("");
! 2560: print("n_wishartd.prob_by_hgm(M,N,[P1,P2,...],[S1,S2,...],T|options) : computation of the probability Pr[l1<T|Sigma].");
! 2561: print(" Arguments : M is the number of variables, N is the degrees of freedom,");
! 2562: print(" Pi is the size of i-th block and Si is the value of i-th (repeated) eigenvalue of Sigma.");
! 2563: print(" Options : step=k => set the number of initial steps=k (default : 10^4)");
! 2564: print(" t0=val => start HGM from t=val (default : max(t0/(2Si),i=1,...)=1)");
! 2565: print(" eps=e => set the relative error bound=e (default : 10^(-4))");
! 2566: print("");
! 2567: print("n_wishartd.message(onoff) : if onoff=1 then various diagnostic messages are shown.");
! 2568: }
! 2569: print("n_wishartd.rr : a package for diagonalization of the matrix 1F1.")$
! 2570: print("n_wishartd.help() shows brief description of some important functions.")$
1.1 noro 2571: endmodule$
2572: end$
1.3 noro 2573:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to wishartd.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