/* $OpenXM: OpenXM/src/asir-contrib/testing/noro/wishartd.rr,v 1.3 2016/04/26 00:55:02 noro Exp $ */
/* A package for computing PDEs satisfied by a matrix 1F1 on diagonal regions */
if (version()<20160401) {error("The version of Risa/Asir must be after 20160401 to run this package.");} else{}
module n_wishartd$
localf compf1$
localf my_comp_f$
localf compc1, compc, compd, compt, addpf, subpf, addc, ftorat, ctorat, mulf$
localf lcmf, mulf, divf, mulc, mulcpf, diffc, muldpf, normalizef1, normalizec1, normalizepf$
localf wsetup, mtot, wishartpf, degf, removef, subst_f, lopitalpf, simpop, simppf$
localf eliminatetop, diagpf, diag0pf, adjop1, adjop, reducepf, reduceotherspf$
localf print_f, printc, printt, printpf, ftop, ctord, comppfrd, multpf, spolypf, nfpf$
localf nfpf0,nfpf1,nfpf_zero$
localf substc$
localf dpf,dpf2,dpf3,abs$
localf pftord$
localf lopital_others_pf$
localf pftop$
localf dnc,pftozpf,addzpf,zpftopf$
localf mul_lopitalpf$
localf indep_eqs$
localf mulpf$
localf pftoeuler$
localf gammam,prc,prc2$
localf ps,psvalue,act_op,decomp_op,create_homo,bpsvalue$
localf pfaffian,gen_basis$
localf evalpf,eval_pfaffian,genrk45$
localf solve_leq,ptom$
localf msubst$
localf vton,encodef1,encodef,encodec1,encodec$
localf vton,ftotex,ctotex,pftotex,optotex,eqtotex$
localf genprc,prob_by_hgm,prob_by_ps$
localf partition, all_partition, coef_partition, count_dup,simp_power1$
localf lop1,lopn,loph,mul_lopitalpf2$
localf muldmpf$
localf diagcheck,eliminatecheck$
localf getchi$
localf message$
localf pfaffian2, eval_pfaffian2$
localf mapleout$
localf gbcheck$
localf partition,all_partition,partition_to_pattern$
localf degpf,pftodpf,all_diag,y1y2$
/*
* pfpoly = [[C,<<...>>],...]
* C = [[N,L],...] L = [[F,M],...]
* lc(F) = 1
*/
static Print$
static PF_N,PF_V,PF_DV$
static Tlopital,Tred,Tlineq,Tactmul$
static Tp,Tps,Trk$
/* XXX execute ord([y1,y2,...,dy1,dy2,...]) */
/* F1>F2 <=> var(F1)>var(F2) || var(F1)=var(F2) && rest(F1)>rest(F2) */
/* F1^M1 > F2^M2 <=> F1>F2 || F1=F2 && M1>M2 */
def abs(A) { return A>0 ? A : -A; }
def compf1(E1,E2)
{
F1 = E1[0]; F2 = E2[0];
if ( F1>F2 ) return 1;
if ( F1<F2 ) return -1;
M1 = E1[1]; M2 = E2[1];
if ( M1 > M2 ) return 1;
if ( M1 < M2 ) return -1;
return 0;
}
def compc1(ND1,ND2)
{
if ( R = comp_f(ND1[1],ND2[1]) ) return R;
N1 = ND1[0]; N2 = ND2[0];
if ( N1 > N2 ) return 1;
if ( N1 < N2 ) return -1;
return 0;
}
/* compare ND lists */
def compc(L1,L2)
{
for ( ; L1 != [] && L2 != []; L1 = cdr(L1), L2 = cdr(L2) ) {
E1 = car(L1); E2 = car(L2);
if ( R = compc1(E1,E2) ) return R;
}
if ( L1 != [] ) return 1;
if ( L2 != [] ) return -1;
return 0;
}
def compd(D1,D2)
{
if ( D1 > D2 ) return 1;
if ( D1 < D2 ) return -1;
return 0;
}
/* compare terms */
def compt(T1,T2)
{
if ( R = compd(T1[1],T2[1]) ) return R;
if ( R = compc(T1[0],T2[0]) ) return R;
return 0;
}
def addpf(F1,F2)
{
R = [];
while ( F1 != [] && F2 != [] ) {
T1 = car(F1); T2 = car(F2);
C = compd(T1[1],T2[1]);
if ( C > 0 ) {
R = cons(T1,R); F1 = cdr(F1);
} else if ( C < 0 ) {
R = cons(T2,R); F2 = cdr(F2);
} else {
S = addc(T1[0],T2[0]);
if ( S != [] )
R = cons([S,T1[1]],R);
F1 = cdr(F1); F2 = cdr(F2);
}
}
R = reverse(R);
if ( F1 != [] ) R = append(R,F1);
else if ( F2 != [] ) R = append(R,F2);
return R;
}
def subpf(F1,F2)
{
T = mulcpf([[-1,[]]],F2);
T = addpf(F1,T);
return T;
}
def addc(F1,F2)
{
R = [];
while ( F1 != [] && F2 != [] ) {
T1 = car(F1); T2 = car(F2);
C = comp_f(T1[1],T2[1]);
if ( !T1[0] || !T2[0] )
error("afo");
if ( C > 0 ) {
R = cons(T1,R); F1 = cdr(F1);
} else if ( C < 0 ) {
R = cons(T2,R); F2 = cdr(F2);
} else {
S = T1[0]+T2[0];
if ( S )
R = cons([S,T1[1]],R);
F1 = cdr(F1); F2 = cdr(F2);
}
}
R = reverse(R);
if ( F1 != [] ) R = append(R,F1);
else if ( F2 != [] ) R = append(R,F2);
return R;
}
def ftorat(F)
{
R = 1;
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
R *= F0[0]^F0[1];
}
return R;
}
def ctorat(C)
{
R = 0;
for ( ; C != []; C = cdr(C) ) {
C0 = car(C);
R += C0[0]/ftorat(C0[1]);
R = red(R);
}
return R;
}
def mulf(L1,L2)
{
R = [];
while ( L1 != [] && L2 != [] ) {
A1 = car(L1); A2 = car(L2);
if ( A1[0] > A2[0] ) {
R = cons(A1,R); L1 = cdr(L1);
} else if ( A1[0] < A2[0] ) {
R = cons(A2,R); L2 = cdr(L2);
} else {
R = cons([A1[0],A1[1]+A2[1]],R);
L1 = cdr(L1); L2 = cdr(L2);
}
}
R = reverse(R);
if ( L2 == [] ) return append(R,L1);
else if ( L1 == [] ) return append(R,L2);
else return R;
}
def lcmf(L1,L2)
{
R = [];
while ( L1 != [] && L2 != [] ) {
A1 = car(L1); A2 = car(L2);
if ( A1[0] > A2[0] ) {
R = cons(A1,R); L1 = cdr(L1);
} else if ( A1[0] < A2[0] ) {
R = cons(A2,R); L2 = cdr(L2);
} else {
M = A1[1]>A2[1]?A1[1]:A2[1];
R = cons([A1[0],M],R);
L1 = cdr(L1); L2 = cdr(L2);
}
}
R = reverse(R);
if ( L2 == [] ) return append(R,L1);
else if ( L1 == [] ) return append(R,L2);
else return R;
}
def mulf(L1,L2)
{
R = [];
while ( L1 != [] && L2 != [] ) {
A1 = car(L1); A2 = car(L2);
if ( A1[0] > A2[0] ) {
R = cons(A1,R); L1 = cdr(L1);
} else if ( A1[0] < A2[0] ) {
R = cons(A2,R); L2 = cdr(L2);
} else {
R = cons([A1[0],A1[1]+A2[1]],R);
L1 = cdr(L1); L2 = cdr(L2);
}
}
R = reverse(R);
if ( L2 == [] ) return append(R,L1);
else if ( L1 == [] ) return append(R,L2);
else return R;
}
def divf(L1,L2)
{
R = [];
while ( L1 != [] && L2 != [] ) {
A1 = car(L1); A2 = car(L2);
if ( A1[0] > A2[0] ) {
R = cons(A1,R); L1 = cdr(L1);
} else if ( A1[0] < A2[0] ) {
error("divf : cannot happen");
} else {
M = A1[1]-A2[1];
if ( M < 0 )
error("divf : cannot happen");
if ( M > 0 )
R = cons([A1[0],M],R);
L1 = cdr(L1); L2 = cdr(L2);
}
}
R = reverse(R);
if ( L2 == [] ) return append(R,L1);
else if ( L1 == [] ) return append(R,L2);
else return R;
}
def mulc(C1,C2)
{
R = [];
C1 = reverse(C1);
for ( ; C1 != []; C1 = cdr(C1) ) {
S = [];
A1 = car(C1);
for ( T = C2 ; T != []; T = cdr(T) ) {
A2 = car(T);
S = cons([A1[0]*A2[0],mulf(A1[1],A2[1])],S);
}
S = reverse(S);
R = addc(R,S);
}
return R;
}
def vton(V)
{
for ( I = 1; I <= PF_N; I++ )
if ( V == PF_V[I] ) return I;
error("vton : no such variable");
}
def encodef1(F)
{
A = F[0];
V1 = var(A);
N1 = vton(V1);
R = A-V1;
if ( !R )
return [N1,F[1]];
else
return [N1*65536+vton(var(R)),F[1]];
}
def encodef(F)
{
return map(encodef1,F);
}
def encodec1(C)
{
return cons(C[0],encodef(C[1]));
}
def encodec(C)
{
R = map(encodec1,C);
}
def mulcpf(C,F)
{
R = [];
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
C1 = mulc(C,F0[0]);
R = cons([C1,F0[1]],R);
}
return reverse(R);
}
def mulpf(F1,F2)
{
R = [];
N = length(PF_V);
for ( T = reverse(F1); T != []; T = cdr(T) ) {
T0 = car(T); C = T0[0]; Op = T0[1];
E = dp_etov(Op);
S = F2;
for ( I = 0; I < N; I++ )
for ( J = 0; J < E[I]; J++ ) S = muldpf(PF_V[I],S);
S = mulcpf(C,S);
R = addpf(R,S);
}
return S;
}
def diffc(X,C)
{
R = [];
for ( T = C; T != []; T = cdr(T) ) {
T0 = car(T);
N = T0[0];
S = [];
for ( L = T0[1]; L != []; S = cons(car(L),S), L = cdr(L) ) {
L0 = car(L);
F = L0[0]; M = L0[1];
DF = diff(F,X);
if ( DF ) {
V = cons([F,M+1],cdr(L));
for ( U = S; U != []; U = cdr(U) ) V = cons(car(U),V);
R = addc([[-M*N*DF,V]],R);
}
}
}
return R;
}
def muldpf(X,F)
{
R = [];
DX = dp_ptod(strtov("d"+rtostr(X)),PF_DV);
for ( T = F; T != []; T = cdr(T) ) {
T0 = car(T);
C = diffc(X,T0[0]);
if ( C != [] )
R = addpf(R,[[T0[0],T0[1]*DX],[C,T0[1]]]);
else
R = addpf(R,[[T0[0],T0[1]*DX]]);
}
return R;
}
/* d^m*c = sum_{i=0}^m m!/i!/(m-i)!*c^(i)*d^(m-i) */
def muldmpf(X,M,F)
{
R = [];
DX = strtov("d"+rtostr(X));
FM = fac(M);
for ( T = F; T != []; T = cdr(T) ) {
T0 = car(T);
C = T0[0]; Op = T0[1];
U = [];
for ( I = 0; I <= M; I++ ) {
if ( C == [] ) break;
C0 = [[FM/fac(I)/fac(M-I),[]]];
U = cons([mulc(C0,C),dp_ptod(DX^(M-I),PF_DV)*Op],U);
C = diffc(X,C);
}
U = reverse(U);
R = addpf(U,R);
}
return R;
}
def normalizef1(F)
{
if ( coef(F,1,var(F)) < 0 ) F = -F;
return F;
}
def normalizec1(C)
{
N = C[0]; L = C[1];
S = 1;
R = [];
for ( ; L != []; L = cdr(L) ) {
L0 = car(L);
F = L0[0]; M = L0[1];
if ( coef(F,1,var(F)) < 0 ) {
F = -F;
if ( M%2 ) S = -S;
}
R = cons([F,M],R);
}
return [S*N,reverse(qsort(R,n_wishartd.compf1))];
}
def normalizepf(F)
{
R = [];
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
C = map(normalizec1,F[0]);
R = cons([C,F[1]],R);
}
return reverse(R);
}
def substc(C,Y2,Y1)
{
R = [];
for ( T = C; T != []; T = cdr(T) ) {
T0 = car(T); N = T0[0]; D = T0[1];
Z = subst_f(D,Y2,Y1);
R = addc([[Z[0]*N,Z[1]]],R);
}
return R;
}
/* dy0 : dummy variable for dy */
def wsetup(N)
{
V = [];
DV = [];
for ( I = N; I>= 0; I-- ) {
YI = strtov("y"+rtostr(I));
DYI = strtov("dy"+rtostr(I));
V = cons(YI,V);
DV = cons(DYI,DV);
}
PF_N = N;
PF_V = V;
PF_DV = DV;
ord(append(V,DV));
}
def mtot(M,Dn)
{
D = dp_ptod(M,PF_DV);
F = fctr(Dn); C = car(F)[0];
Dn = reverse(qsort(cdr(F),n_wishartd.compf1));
return [[[dp_hc(D)/C,Dn]],dp_ht(D)];
}
def wishartpf(N,I)
{
YI = PF_V[I]; DYI = PF_DV[I];
R = [mtot(DYI^2,1)];
R = addpf([mtot(c*DYI,YI)],R);
R = addpf([mtot(-DYI,1)],R);
for ( J = 1; J <= N; J++ ) {
if ( J == I ) continue;
YJ = PF_V[J]; DYJ = PF_DV[J];
R = addpf([mtot(1/2*DYI,YI-YJ)],R);
R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
R = addpf([mtot(-1/2*DYI,YI)],R);
R = addpf([mtot(1/2*DYJ,YI)],R);
}
R = addpf([mtot(-a,YI)],R);
return R;
}
def degf(F,D)
{
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
if ( F0[0] == D ) return F0[1];
}
return 0;
}
def removef(F,D)
{
R = [];
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
if ( F0[0] != D ) R = cons(F0,R);
}
return reverse(R);
}
def subst_f(F,Y2,Y1)
{
R = [];
Sgn = 1;
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
T = subst(F0[0],Y2,Y1);
if ( coef(T,1,var(T)) < 0 ) {
T = -T;
if ( F0[1]%2 ) Sgn = -Sgn;
}
R = cons([T,F0[1]],R);
}
if ( R == [] ) return [Sgn,R];
R = qsort(R,n_wishartd.compf1);
S = [car(R)]; R = cdr(R);
while ( R != [] ) {
R0 = car(R);
S0 = car(S);
if ( R0[0] == S0[0] )
S = cons([S0[0],S0[1]+R0[1]],cdr(S));
else
S = cons(R0,S);
R = cdr(R);
}
return [Sgn,S];
}
/* Y2 -> Y1 */
def lopitalpf(P,Y1,Y2)
{
// if ( !member(Y2,vars(P)) ) return P;
D = normalizef1(Y2-Y1);
if ( D == Y2-Y1 ) Sgn = 1;
else Sgn = -1;
DY2 = strtov("d"+rtostr(Y2));
R = [];
for ( T = reverse(P); T != []; T = cdr(T) ) {
T0 = car(T);
C = T0[0]; Op = T0[1];
for ( U = reverse(C); U != []; U = cdr(U) ) {
U0 = car(U);
Nm = U0[0]; Dn = U0[1];
K = degf(Dn,D);
if ( !K ) {
Z = subst_f(Dn,Y2,Y1);
Sgn1 = Z[0]; Dn1 = Z[1];
R = addpf([[[[Sgn1*Nm,Dn1]],Op]],R);
} else {
Dn1 = removef(Dn,D);
C1 = [[1,Dn1]];
Z = [];
for ( J = 0; J <= K; J++ ) {
B = fac(K)/fac(J)/fac(K-J);
C1s = substc(C1,Y2,Y1);
if ( C1s != [] ) {
W = [[C1s,dp_ptod(DY2^(K-J),PF_DV)*Op]];
W = mulcpf([[B,[]]],W);
Z = addpf(W,Z);
}
C1 = diffc(Y2,C1);
if ( C1 == [] ) break;
}
Z = mulcpf([[Sgn^K*Nm/fac(K),[]]],Z);
R = addpf(Z,R);
}
}
}
return R;
}
def lopital_others_pf(P,L) {
Q = P;
for ( T = L; T != []; T = cdr(T) ) {
T0 = car(T);
I0 = T0[0]; I1 = T0[1]; I = I1-I0+1;
for ( M = I0; M <= I1; M++ ) {
Q = lopitalpf(Q,PF_V[I0],PF_V[M]);
}
Q = simppf(Q,I0,I);
Q = adjop(Q,I0,I);
}
return Q;
}
def simpop(Op,I0,NV)
{
E = dp_etov(Op);
R = [];
for ( J = 0, I = I0; J < NV; I++, J++ ) R = cons(E[I],R);
R = reverse(qsort(R));
for ( J = 0, I = I0; J < NV; I++, J++ ) E[I] = R[J];
return dp_vtoe(E);
}
def simppf(P,I0,NV)
{
R = [];
for ( T = P; T != []; T = cdr(T) ) {
T0 = car(T);
R = addpf([[T0[0],simpop(T0[1],I0,NV)]],R);
}
return R;
}
/* simplify (dy1+...+dyn)^k */
def simp_power1(K,I0,NV)
{
P = all_partition(K,NV);
M = map(coef_partition,P,K,NV,I0);
for ( R = 0, T = M; T != []; T = cdr(T) )
R += car(T);
return R;
}
def indep_eqs(Eq)
{
M = length(Eq);
D = 0;
for ( I = 0; I < M; I++ )
for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
for ( N = 0, T = D; T; T = dp_rest(T), N++ );
Terms = vector(N);
for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
A = matrix(M,N);
for ( I = 0; I < M; I++ ) {
for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
T = car(H)[1];
for ( K = 0; K < N; K++ )
if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
}
}
for ( I = 0; I < M; I++ ) {
Dn = 1;
for ( J = 0; J < N; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
for ( J = 0; J < N; J++ ) A[I][J] *= Dn;
}
R = indep_rows_mod(A,lprime(0));
if ( length(R) == N ) {
L = [];
for ( I = N-1; I >= 0; I-- )
L = cons(Eq[R[I]],L);
return L;
} else
return 0;
}
def eliminatetop(Eq)
{
Eq = indep_eqs(Eq);
if ( !Eq ) return 0;
M = length(Eq);
R = vector(M);
D = 0;
for ( I = 0; I < M; I++ )
for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
for ( N = 0, T = D; T; T = dp_rest(T), N++ );
if ( N != M )
return 0;
N2 = 2*N;
Terms = vector(N);
for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
A = matrix(N,N2);
for ( I = 0; I < N; I++ ) {
for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
T = car(H)[1];
for ( K = 0; K < N; K++ )
if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
}
A[I][N+I] = 1;
}
for ( I = 0; I < N; I++ ) {
Dn = 1;
for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
}
Ret = generic_gauss_elim(A);
Num = Ret[0]; Den = Ret[1]; Ind0 = Ret[2]; Ind1 = Ret[3];
if ( length(Ind0) != N || Ind0[N-1] != N-1 )
return 0;
R = [];
if ( Print ) print(["N=",N]);
for ( I = 0; I < N; I++ ) {
if ( Print ) print(["I=",I],2);
T = [];
for ( L = 0; L < N; L++ )
if ( Num[I][L] ) T = addpf(T,mulcpf([[Num[I][L]/Den,[]]],Eq[L][1]));
R = cons([Terms[I],T],R);
}
if ( Print ) print("");
R = qsort(R,n_wishartd.compt);
return reverse(R);
}
def eliminatecheck(Eq,Top)
{
Eq = indep_eqs(Eq);
if ( !Eq ) return 0;
M = length(Eq);
R = vector(M);
D = 0;
for ( I = 0; I < M; I++ )
for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
for ( N = 0, T = D; T; T = dp_rest(T), N++ );
if ( N != M )
return 0;
N2 = 2*N;
Terms = vector(N);
for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
A = matrix(N,N2);
for ( I = 0; I < N; I++ ) {
for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
T = car(H)[1];
for ( K = 0; K < N; K++ )
if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
}
A[I][N+I] = 1;
}
for ( I = 0; I < N; I++ ) {
Dn = 1;
for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
}
if ( Top ) {
for ( I = 0; I < N; I++ ) {
for ( T = J = 0; J < N; J++ ) if ( J != I ) T += abs(A[I][J]);
if ( abs(A[I][I]) < T )
if ( Print ) print(["ng",I]);
}
} else {
for ( I = 1; I < N; I++ ) {
for ( T = J = 0; J < N-2; J++ ) if ( J != I-1 ) T += abs(A[I][J]);
if ( abs(A[I][I-1]) < T )
if ( Print ) print(["ng",I]);
}
}
Ret = generic_gauss_elim_mod(A,lprime(0));
Num = Ret[0]; Ind0 = Ret[1]; Ind1 = Ret[2];
if ( length(Ind0) != N || Ind0[N-1] != N-1 )
return 0;
else
return 1;
}
def mul_lopitalpf(Z,N,K,I0,I1,E)
{
I = I1-I0+1;
for ( J = I0; J <= N; J++ )
for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
for ( J = I0+1; J <= I1; J++ ) {
Z = lopitalpf(Z,PF_V[I0],PF_V[J]);
}
S = simppf(Z,I0,I);
H = []; L = [];
for ( ; S != []; S = cdr(S) ) {
if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
else L = cons(car(S),L);
}
return [reverse(H),reverse(L)];
}
/* Blocks = [B1,B2,...] */
/* Bi=[s,e] ys=...=ye f */
def diag0pf(N,Blocks)
{
Tlopital = Tred = Tlineq = 0;
wsetup(N);
P = vector(N+1);
for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
Simp = [];
Done = [];
Len = length(Blocks);
for ( A = 0; A < Len; A++ ) {
B = Blocks[A];
Blocks0 = setminus(Blocks,[B]);
I0 = B[0];
I1 = B[1];
I = I1-I0+1;
for ( K = I0; K <= I1; K++ )
P[K] = lopital_others_pf(P[K],Blocks0);
Rall = [];
for ( K = I+1; K >= 2; K-- ) {
if ( Print ) print(["K=",K],2);
DK = simp_power1(K,I0,I);
Eq = [];
TlopitalK = 0;
for ( T = DK; T; T = dp_rest(T) ) {
E = dp_etov(T);
for ( U = I0; U <= I1; U++ )
if ( E[U] >= 2 )
break;
if ( U > I1 ) continue;
E[U] -= 2;
Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
Eq = cons(Ret,Eq);
}
Eq = reverse(Eq);
for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
DY = mtot(-dy0^K,1);
if ( K == I+1 ) {
EqTop = addpf([DY],Eq0);
} else {
H = []; L = [];
for ( S = Eq0 ; S != []; S = cdr(S) ) {
if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
else L = cons(car(S),L);
}
L = reverse(L); H = reverse(H);
Eq = cons([H,addpf([DY],L)],Eq);
}
T0 = time();
R = eliminatetop(Eq);
if ( R )
Rall = cons(R,Rall);
else
return [];
T1 = time(); Tlineq += T1[0]-T0[0];
if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
}
Rall = reverse(Rall);
/* EqTop : dyi0 -> dy -> dyi0 */
for ( T = Rall; T != []; T = cdr(T) ) {
if ( Print ) print(["red",length(T)+1],2);
T0 = time();
EqTop = reducepf(EqTop,car(T));
T1 = time(); Tred += T1[0]-T0[0];
if ( Print ) print([T1[0]-T0[0],"sec"]);
}
EqTop = adjop(EqTop,I0,I);
Done = cons(EqTop,Done);
}
if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
Done = ltov(Done);
Len = length(Done);
for ( I = 0; I < Len; I++ ) {
for ( G = [], J = 0; J < Len; J++ )
if ( J != I ) G = cons(Done[J],G);
Done[I] = nfpf(Done[I],G);
}
return vtol(Done);
}
def diagpf(N,Blocks)
{
Tlopital = Tred = Tlineq = 0;
wsetup(N);
P = vector(N+1);
for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
Simp = [];
Done = [];
Len = length(Blocks);
for ( A = 0; A < Len; A++ ) {
B = Blocks[A];
Blocks0 = setminus(Blocks,[B]);
I0 = B[0];
I1 = B[1];
I = I1-I0+1;
for ( K = I0; K <= I1; K++ )
P[K] = lopital_others_pf(P[K],Blocks0);
Rall = [];
for ( K = I+1; K >= 2; K-- ) {
if ( Print ) print(["K=",K],2);
DK = simp_power1(K,I0,I);
Eq = [];
TlopitalK = 0;
for ( T = DK; T; T = dp_rest(T) ) {
E = dp_etov(T);
for ( U = I1; U >= I0; U-- )
if ( E[U] >= 2 )
break;
if ( U < I0 ) continue;
E[U] -= 2;
Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
Eq = cons(Ret,Eq);
}
Eq = reverse(Eq);
for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
DY = mtot(-dy0^K,1);
if ( K == I+1 ) {
EqTop = addpf([DY],Eq0);
} else {
H = []; L = [];
for ( S = Eq0 ; S != []; S = cdr(S) ) {
if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
else L = cons(car(S),L);
}
L = reverse(L); H = reverse(H);
Eq = cons([H,addpf([DY],L)],Eq);
}
T0 = time();
R = eliminatetop(Eq);
if ( R )
Rall = cons(R,Rall);
else {
if ( Print ) print(["lineq retry.."],2);
Eq1 = [];
Terms = [];
for ( T = dp_rest(DK); T; T = dp_rest(T) )
Terms = cons(dp_ht(T),Terms);
while ( Terms != [] ) {
for ( Q = 0; Terms != [] && Q < 10; Q++, Terms = cdr(Terms) ) {
E = dp_etov(car(Terms));
if ( E[I0] >= 2 ) {
E[I0] -= 2;
Ret = mul_lopitalpf(P[I0],N,K,I0,I1,E);
Eq1 = cons(Ret,Eq1);
}
}
Eq = append(Eq,Eq1);
R = eliminatetop(Eq);
if ( R ) break;
}
if ( !R )
error("diagpf : failed to find a solvable linear system");
Rall = cons(R,Rall);
}
T1 = time(); Tlineq += T1[0]-T0[0];
if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
}
Rall = reverse(Rall);
/* EqTop : dyi0 -> dy -> dyi0 */
for ( T = Rall; T != []; T = cdr(T) ) {
if ( Print ) print(["red",length(T)+1],2);
T0 = time();
EqTop = reducepf(EqTop,car(T));
T1 = time(); Tred += T1[0]-T0[0];
if ( Print ) print([T1[0]-T0[0],"sec"]);
}
EqTop = adjop(EqTop,I0,I);
Done = cons(EqTop,Done);
}
if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
Done = ltov(Done);
Len = length(Done);
for ( I = 0; I < Len; I++ ) {
for ( G = [], J = 0; J < Len; J++ )
if ( J != I ) G = cons(Done[J],G);
Done[I] = nfpf1(Done[I],G);
}
return vtol(Done);
}
def diagcheck(N)
{
Tmuld = Tlopital = 0;
MHist = [];
wsetup(N);
P = vector(N+1);
for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
Simp = [];
Done = [];
B = [1,N];
I0 = B[0];
I1 = B[1];
I = I1-I0+1;
for ( K = 2; K <= I+1; K++ ) {
if ( Print ) print(["K=",K],2);
DK = simp_power1(K,I0,I);
Eq = [];
TlopitalK = 0;
for ( T = DK; T; T = dp_rest(T) ) {
E = dp_etov(T);
for ( U = I0; U <= I1; U++ )
if ( E[U] >= 2 )
break;
if ( U > I1 ) continue;
E[U] -= 2;
Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E|high=1);
Eq = cons(Ret,Eq);
}
Eq = reverse(Eq);
for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
DY = mtot(-dy0^K,1);
if ( K == I+1 ) {
EqTop = addpf([DY],Eq0);
} else {
H = []; L = [];
for ( S = Eq0 ; S != []; S = cdr(S) ) {
if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
else L = cons(car(S),L);
}
L = reverse(L); H = reverse(H);
Eq = cons([H,addpf([DY],L)],Eq);
}
R = eliminatecheck(Eq,K==I+1);
if ( !R )
return 0;
}
if ( Print ) print(["muld",Tmuld,"lopital",Tlopital]);
return 1;
}
/* dyi -> dyi/K, dy->dyi */
def adjop1(T,I,K)
{
C = T[0];
Op = dp_etov(T[1]);
if ( Op[I] ) C = mulc([[1/K,[]]],C);
Op[I] += Op[0];
Op[0] = 0;
return [C,dp_vtoe(Op)];
}
def adjop(P,I,K)
{
P = map(adjop1,P,I,K);
P = qsort(P,n_wishartd.compt);
return reverse(P);
}
def dnc(C)
{
D = 1;
for ( T = C; T != []; T = cdr(T) ) {
T0 = car(T); N = T0[0];
M = sdiv(ptozp(N),N);
D = ilcm(D,M);
}
return D;
}
def pftozpf(F)
{
D = 1;
for ( T = F; T != []; T = cdr(T) ) {
T0 = car(T);
D = ilcm(D,dnc(T0[0]));
}
return [D,mulcpf([[D,[]]],F)];
}
def zpftopf(F)
{
return mulcpf([[1/F[0],[]]],F[1]);
}
def addzpf(F1,F2)
{
D1 = F1[0]; D2 = F2[0];
L = ilcm(D1,D2); C1 = L/D1; C2 = L/D2;
N = addpf(mulcpf([[L/D1,[]]],F1[1]),mulcpf([[L/D2,[]]],F2[1]));
return [L,N];
}
/* R : reducers */
def reducepf(Eq,R)
{
Ret = pftozpf(Eq); Eq = Ret[1]; Dn = Ret[0];
S = [1,[]];
for ( T = R, U = []; T != []; T = cdr(T) )
U = cons([car(T)[0],pftozpf(car(T)[1])],U);
R = reverse(U);
for ( T = reverse(Eq); T != []; T = cdr(T) ) {
T0 = car(T); C = T0[0]; Op = T0[1];
for ( U = (R); U != []; U = cdr(U) ) {
U0 = car(U);
if ( dp_redble(Op,U0[0]) ) {
E = dp_etov(dp_subd(Op,U0[0]));
Red = U0[1];
N = length(E);
for ( J = 1; J < N; J++ )
for ( L = 0; L < E[J]; L++ ) Red = [Red[0],muldpf(PF_V[J],Red[1])];
Red = [Red[0],mulcpf(C,Red[1])];
Red = [Red[0],mulcpf([[-1,[]]],Red[1])];
S = addzpf(S,Red);
break;
}
}
if ( U == [] ) S = addzpf([1,[T0]],S);
}
S = [S[0]*Dn,S[1]];
return zpftopf(S);
}
def reduceotherspf(Eq,P,I1,N)
{
S = [];
Z = Eq;
while ( Z != [] ) {
T0 = car(Z); C = T0[0]; Op = T0[1];
for ( I = I1; I <= N; I++ ) {
U0 = P[I];
M = car(U0)[0][0][0];
H = car(U0)[1];
if ( dp_redble(Op,H) ) {
E = dp_etov(dp_subd(Op,H));
Red = U0;
Len = length(E);
for ( J = 0; J < Len; J++ )
for ( L = 0; L < E[J]; L++ ) Red = muldpf(PF_V[J],Red);
C1 = mulc([[1/M,[]]],C);
Red = mulcpf(C1,Red);
Z = subpf(Z,Red);
break;
}
}
if ( I > N ) {
S = cons(T0,S);
Z = cdr(Z);
}
}
return reverse(S);
}
def print_f(F)
{
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
print("(",0); print(F0[0],0); print(")",0);
if ( F0[1] > 1 ) {
print("^",0); print(F0[1],0);
}
if ( cdr(F) != [] ) print("*",0);
}
}
def printc(C)
{
print("(",0);
for ( ; C != []; C = cdr(C) ) {
C0 = car(C);
print("(",0); print(C0[0],0); print(")",0);
if ( C0[1] != [] ) {
print("/(",0);
print_f(C0[1]);
print(")",0);
}
if ( cdr(C) != [] ) print("+",0);
}
print(")",0);
}
def printt(T)
{
printc(T[0]); print("*",0); print(dp_dtop(T[1],PF_DV),0);
}
def printpf(F)
{
for ( ; F != []; F = cdr(F) ) {
printt(car(F));
if ( cdr(F) != [] ) print("+",0);
}
print("");
}
def ftop(F)
{
P = 1;
for ( ; F != []; F = cdr(F) ) {
F0 = car(F);
P *= F0[0]^F0[1];
}
return P;
}
def pftord(F)
{
R = [];
for ( T = F; T != []; T = cdr(T) ) {
T0 = car(T); C = T0[0]; Op = T0[1];
R = cons([ctord(C),Op],R);
}
return reverse(R);
}
def pftop(F)
{
R = pftord(F);
Op = F[0][1];
N = length(dp_etov(Op));
DY = [];
for ( I = N; I >= 0; I-- ) DY = cons(strtov("dy"+rtostr(I)),DY);
Lcm = [];
for ( T = R; T != []; T = cdr(T) )
Lcm = lcmf(Lcm,T[0][0][1]);
S = 0;
for ( T = R; T != []; T = cdr(T) ) {
N = T[0][0][0]*ftop(divf(Lcm,T[0][0][1]));
Op = dp_dtop(T[0][1],DY);
S += N*Op;
}
return S;
}
def ctord(C)
{
N = 0; D = [];
Len = length(C);
for ( I = 0, T = reverse(C); T != []; T = cdr(T),I++ ) {
// print([Len-I],2);
T0 = car(T); N0 = T0[0]; D0 = T0[1];
L = lcmf(D,D0);
/* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
N = N*ftop(divf(L,D))+N0*ftop(divf(L,D0));
D = L;
}
R = [];
for ( S = D; S != []; S = cdr(S) ) {
T0 = car(S); F = T0[0]; M = T0[1];
while ( M > 0 && (Q = tdiv(N,F)) ) {
N = Q;
M--;
}
if ( M ) R = cons([F,M],R);
}
D = reverse(R);
return [N,D];
}
def comppfrd(P,R)
{
P = qsort(P,n_wishartd.compt); P=reverse(P);
R = qsort(R,n_wishartd.compt); R=reverse(R);
if ( length(P) != length(R) ) return 0;
for ( ; P != []; P = cdr(P), R = cdr(R) ) {
P0 = car(P); R0 = car(R);
C0 = ctord(P0[0]); C1 = R0[0];
D0 = ftop(C0[1]); D1 = ftop(C1[1]);
if ( (D0 == D1) && C0[0] == -C1[0] ) continue;
if ( (D0 == -D1) && C0[0] == C1[0] ) continue;
return 0;
}
return 1;
}
def multpf(T,F)
{
E = dp_etov(T[1]);
N = length(E);
Z = F;
for ( J = 1; J < N; J++ )
for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
Z = mulcpf(T[0],Z);
return Z;
}
def spolypf(F,G)
{
F0 = car(F); G0 = car(G);
DF = F0[1]; DG = G0[1];
L = dp_lcm(DF,DG);
F1 = multpf([F0[0],dp_subd(L,DF)],F);
G1 = multpf([G0[0],dp_subd(L,DG)],G);
S = subpf(F1,G1);
return S;
}
def nfpf(F,G)
{
Z = F;
R = [];
while ( Z != [] ) {
Z0 = car(Z); C = Z0[0]; D = Z0[1];
for ( T = G; T != []; T = cdr(T) ) {
T0 = car(T);
if ( dp_redble(D,T0[0][1]) ) break;
}
if ( T != [] ) {
CG = ctorat(T0[0][0]);
C = mulc(C,[[-1/CG,[]]]);
S = multpf([C,dp_subd(D,T0[0][1])],T0);
Z = addpf(Z,S);
} else {
R = cons(Z0,R);
Z = cdr(Z);
}
}
S = [];
for ( T = R; T != []; T = cdr(T) ) {
U = ctord(T[0][0]);
if ( U[0] )
S = cons(T[0],S);
}
return S;
}
def nfpf0(F,G)
{
Z = F;
R = [];
while ( Z != [] ) {
Z0 = car(Z); C = Z0[0]; D = Z0[1];
for ( T = G; T != []; T = cdr(T) ) {
T0 = car(T);
if ( dp_redble(D,T0[0][1]) ) break;
}
if ( T != [] ) {
CG = ctorat(T0[0][0]);
C = mulc(C,[[-1/CG,[]]]);
S = multpf([C,dp_subd(D,T0[0][1])],T0);
Z = addpf(Z,S);
} else {
R = cons(Z0,R);
Z = cdr(Z);
}
}
S = [];
for ( T = R; T != []; T = cdr(T) ) {
U = ctord(T[0][0]);
if ( U[0] )
S = cons(T[0],S);
}
return S;
}
def nfpf1(F,G)
{
Z = F;
R = [];
while ( Z != [] ) {
Z0 = car(Z); C = Z0[0]; D = Z0[1];
for ( T = G; T != []; T = cdr(T) ) {
T0 = car(T);
if ( dp_redble(D,T0[0][1]) ) break;
}
if ( T != [] ) {
CG = ctorat(T0[0][0]);
C = mulc(C,[[-1/CG,[]]]);
S = multpf([C,dp_subd(D,T0[0][1])],T0);
Z = addpf(Z,S);
} else {
R = cons(Z0,R);
Z = cdr(Z);
}
}
return reverse(R);
}
def nfpf_zero(F,G)
{
Z = F;
R = [];
while ( Z != [] ) {
Z0 = car(Z); C = Z0[0]; D = Z0[1];
U = ctord(C);
if ( !U[0] ) { Z = cdr(Z); continue;}
for ( T = G; T != []; T = cdr(T) ) {
T0 = car(T);
if ( dp_redble(D,T0[0][1]) ) break;
}
if ( T != [] ) {
CG = ctorat(T0[0][0]);
C = mulc(C,[[-1/CG,[]]]);
S = multpf([C,dp_subd(D,T0[0][1])],T0);
Z = addpf(Z,S);
} else {
R = cons(Z0,R);
Z = cdr(Z);
}
}
S = [];
for ( T = R; T != []; T = cdr(T) ) {
U = ctord(T[0][0]);
if ( U[0] )
S = cons(T[0],S);
}
print("");
return S;
}
def pftoeuler(F)
{
VDV = [y1,dy1];
P = pftop(F);
Z = dp_ptod(P,VDV);
E = dp_etov(dp_ht(Z));
S = E[1]-E[0];
if ( S < 0 )
error("pftoeuler : invalid input");
P *= y1^S;
Z = dp_ptod(P,VDV);
E = dp_etov(dp_ht(Z));
D = E[0];
C = vector(D+1);
C[0] = 1;
for ( I = 1; I <= D; I++ )
C[I] = C[I-1]*(t-I+1);
R = 0;
for ( T = Z; T; T = dp_rest(T) ) {
E = dp_etov(dp_ht(T));
S = E[0]-E[1];
if ( S < 0 )
error("pftoeuler : invalid input");
R += dp_hc(T)*y1^S*C[E[1]];
}
return R;
}
def gammam(M,A)
{
R = 1;
for ( I = 1; I <= M; I++ )
R *= mpfr_gamma(A-(I-1)/2);
R *= eval(@pi^(1/4*M*(M-1)));
return R;
}
def genprc(M,N,PL,SL)
{
A = (M+1)/2; C = (N+M+1)/2;
DetS = 1;
TrIS = 0;
for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
DetS *= car(S)^car(T);
TrIS += car(T)/car(S);
}
C = 1/eval(DetS^(1/2*N)*2^(1/2*N*M));
return C;
}
/*
* PL=[P1,P2,...]: sizes of blocks
* SL=[S1,S2,...]: [S1xP1,S2xP2,..]
*/
def prob_by_hgm(M,N,PL,SL,X)
{
A = (M+1)/2; C = (N+M+1)/2;
B = []; V = []; Beta = []; SBeta = 0;
Prev = 1;
for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
V = cons(strtov("y"+rtostr(Prev)),V);
B = cons([Prev,Prev+car(T)-1],B);
Prev += car(T);
Beta = cons(1/(2*car(S)),Beta);
SBeta += car(T)/(2*car(S));
}
if ( Prev != M+1 )
error("invalid block specification");
B = reverse(B); V = reverse(V);
Beta = reverse(Beta);
T0 = time();
if ( type(Z=getopt(eq)) == -1 )
Z = diagpf(M,B);
T1 = time(); Tdiag = (T1[0]-T0[0])+(T1[1]-T0[1]);
if ( type(Step=getopt(step)) == -1 )
Step = 10000;
Z = subst(Z,a,A,c,C);
if ( type(X0=getopt(x0)) == -1 ) {
if ( type(Init=getopt(init)) == -1 ) Init = 1;
X0 = 0;
Len = length(Beta);
for ( I = 0; I < Len; I++ )
if ( Beta[I] > X0 ) X0 = Beta[I];
X0 = Init/X0;
}
if ( type(Rk=getopt(rk)) == -1 ) Rk = 5;
if ( type(Eps=getopt(eps)) == -1 ) Eps = 10^(-4);
if ( type(MapleOut=getopt(mapleout)) == -1 ) MapleOut = 0;
F = genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,Rk,Step,X0,Eps|mapleout=MapleOut);
if ( Print ) print(["diag",Tdiag,"pfaffian",Tp,"ps",Tps,"rk",Trk]);
return F;
}
def prob_by_ps(M,N,PL,SL,X)
{
A = (M+1)/2; C = (N+M+1)/2;
B = []; V = []; Beta = []; SBeta = 0;
Prev = 1;
for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
V = cons(strtov("y"+rtostr(Prev)),V);
B = cons([Prev,Prev+car(T)-1],B);
Prev += car(T);
Beta = cons(1/(2*car(S)),Beta);
SBeta += car(T)/(2*car(S));
}
if ( Prev != M+1 )
error("invalid block specification");
B = reverse(B); V = reverse(V);
Beta = reverse(Beta);
if ( type(Z=getopt(eq)) == -1 )
Z = diagpf(M,B);
Z = subst(Z,a,A,c,C);
GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
C = GM*genprc(M,N,PL,SL);
Beta = ltov(Beta)*eval(exp(0));
X *= eval(exp(0));
L = psvalue(Z,V,Beta*X); PS = L[0];
MN2 = M*N/2;
ExpF = eval(X^(MN2)*exp(-SBeta*X));
return C*L[1]*ExpF;
}
def ps(Z,V,TD)
{
Tact = Tgr = Tactmul = 0;
G = map(ptozp,map(pftop,Z));
M = length(V);
N = length(G);
for ( I = 0, DV = []; I < M; I++ )
DV = cons(strtov("d"+rtostr(V[I])),DV);
DV = reverse(DV);
VDV = append(V,DV);
DG = map(decomp_op,G,V,DV);
F = [1];
R = 1;
Den = 1;
for ( I = 1 ; I <= TD; I++ ) {
L = create_homo(V,I);
FI = L[0]; CV = L[1];
Eq = [];
for ( K = 0; K < N; K++ ) {
P = DG[K]; Len = length(P);
/* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
T0=time();
Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
T1=time(); Tact += T1[0]-T0[0];
if ( Lhs ) {
for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
Eq = cons(dp_hc(T),Eq);
}
}
}
T0=time();
#if 0
Sol = solve_leq(Eq,CV);
#else
Sol = nd_f4(Eq,CV,0,0);
#endif
L = p_true_nf(FI,Sol,CV,0);
Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
FI = L[0]*(Den1/L[1]);
T1=time(); Tgr += T1[0]-T0[0];
MJ = Den1/Den; Den = Den1;
for ( S = [], T = F; T != []; T = cdr(T) )
S = cons(MJ*car(T),S);
F = cons(FI,reverse(S));
R = R*MJ+car(F);
}
return R/Den;
}
def msubst(F,V,X)
{
M = length(V);
for ( J = 0, F0 = F*eval(exp(0)); J < M; J++ )
F0 = subst(F0,V[J],X[J]);
return F0;
}
def psvalue(Z,V,X)
{
Tact = Tgr = Tactmul = 0;
G = map(ptozp,map(pftop,Z));
M = length(V);
N = length(G);
for ( I = 0, DV = []; I < M; I++ )
DV = cons(strtov("d"+rtostr(V[I])),DV);
DV = reverse(DV);
VDV = append(V,DV);
DG = map(decomp_op,G,V,DV);
Prev = getopt(prev);
if ( type(Prev) == -1 ) {
F = [1];
R = 1;
Val = eval(exp(0));
Den = 1;
I = 1;
} else {
/* Prev = [R/Den,Val1,F,Den] */
Den = Prev[3];
F = Prev[2];
R = Prev[0]*Den;
I = length(F);
Val = msubst(Prev[0],V,X);
ValI = msubst(F[0],V,X)/Den;
if ( Val-ValI == Val )
return [Prev[0],Val,F,Den];
}
for ( ; ; I++ ) {
L = create_homo(V,I);
FI = L[0]; CV = L[1];
Eq = [];
for ( K = 0; K < N; K++ ) {
P = DG[K]; Len = length(P);
/* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
T0=time();
Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
T1=time(); Tact += T1[0]-T0[0];
if ( Lhs ) {
for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
Eq = cons(dp_hc(T),Eq);
}
}
}
T0=time();
Sol = solve_leq(Eq,CV);
L = p_true_nf(FI,Sol,CV,0);
delete_uc();
Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
FI = L[0]*(Den1/L[1]);
T1=time(); Tgr += T1[0]-T0[0];
MJ = Den1/Den; Den = Den1;
for ( S = [], T = F; T != []; T = cdr(T) )
S = cons(MJ*car(T),S);
F = cons(FI,reverse(S));
R = R*MJ+car(F);
F0 = msubst(FI,V,X)/Den;
Val1 = Val+F0;
if ( Val1 == Val ) {
if ( Print ) print(["I=",I,"act",Tact,"actmul",Tactmul,"gr",Tgr]);
return [R/Den,Val1,F,Den];
} else {
if ( Print ) print([F0],2);
Val = Val1;
}
}
}
/* P -> [P0,P1,...] Pi = y^a dy^b, |a|-|b|=i */
def decomp_op(P,V,DV)
{
M = length(V);
VDV = append(V,DV);
D = dp_ptod(P,VDV);
for ( I = 0, T = D; T; T = dp_rest(T), I++ );
for ( T = D, NotSet = 1; T; T = dp_rest(T) ) {
E = dp_etov(T);
for ( K = 0, J = 0; J < M; J++ )
K += E[J]-E[M+J];
if ( NotSet ) {
Max = Min = K; NotSet = 0;
} else if ( K > Max ) Max = K;
else if ( K < Min ) Min = K;
}
W = vector(Max-Min+1);
for ( T = D; T; T = dp_rest(T) ) {
E = dp_etov(T);
for ( K = 0, J = 0; J < M; J++ )
K += E[J]-E[M+J];
W[K-Min] += dp_hm(T);
}
W = map(dp_ptod,map(dp_dtop,W,VDV),DV);
return W;
}
def create_homo(V,TD)
{
M = length(V);
for ( S = 0, I = 0; I < M; I++ ) S += V[I];
Tmp = 0;
Nc = 0;
for ( RI = 0, D = dp_ptod(S^TD,V); D; D = dp_rest(D), Nc++ ) {
E = dp_etov(D);
T = uc();
U = dp_dtop(dp_ht(D),V);
RI += T*U;
Tmp += T;
}
CV = vector(Nc);
for ( I = 0; I < Nc; I++ ) {
CV[I] = var(Tmp); Tmp -= CV[I];
}
return [RI,vtol(CV)];
}
def act_op(P,V,F)
{
N = length(V);
for ( R = 0; P; P = dp_rest(P) ) {
C = dp_hc(P); E = dp_etov(P);
for ( I = 0, S = F; I < N; I++ )
for ( J = 0; J < E[I]; J++ ) S = diff(S,V[I]);
T0 = time();
R += C*S;
T1 = time(); Tactmul += T1[0]-T0[0];
}
return R;
}
def gen_basis(D,DV)
{
N = length(D);
B = [];
E = vector(N);
while ( 1 ) {
B = cons(dp_dtop(dp_vtoe(E),DV),B);
for ( I = N-1; I >= 0; I-- )
if ( E[I] < D[I]-1 ) break;
if ( I < 0 ) return reverse(B);
E[I]++;
for ( J = I+1; J < N; J++ ) E[J] = 0;
}
}
def pfaffian(Z)
{
N = length(Z);
D = vector(N);
Mat = vector(N);
V = vars(Z);
DV = vector(N);
for ( I = 0; I < N; I++ )
DV[I] = strtov("d"+rtostr(V[I]));
DV = vtol(DV);
for ( I = 0; I < N; I++ ) {
ZI = Z[I];
HI = ZI[0][1];
DI = dp_dtop(HI,PF_DV);
for ( J = 0; J < N; J++ )
if ( var(DI) == DV[J] )
break;
D[J] = deg(DI,DV[J]);
}
B = gen_basis(D,DV);
Dim = length(B);
Hist = [];
for ( I = 0; I < N; I++ ) {
if ( Print ) print(["I=",I]);
A = matrix(Dim,Dim);
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
A[K][L] = [];
for ( J = 0; J < Dim; J++ ) {
if ( Print ) print(["J=",J],2);
Mon = DV[I]*B[J];
for ( T = Hist; T != []; T = cdr(T) )
if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
if ( (T != []) ) {
for ( L = 0; L < N; L++ )
if ( Q == DV[L] ) break;
Red1 = muldpf(V[L],car(T)[1]);
Red = nfpf0(Red1,Z);
} else
Red = nfpf0([mtot(Mon,1)],Z);
Hist = cons([Mon,Red],Hist);
for ( T = Red; T != []; T = cdr(T) ) {
T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
for ( K = 0; K < Dim; K++ )
if ( B[K] == Op ) break;
if ( K == Dim )
error("afo");
else
A[J][K] = C;
}
}
if ( Print ) print("");
A = map(ctord,A);
Lcm = [];
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
Lcm = lcmf(Lcm,A[K][L][1]);
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ ) {
Num = A[K][L][0]; Den = A[K][L][1];
A[K][L] = Num*ftorat(divf(Lcm,Den));
}
Lcm = ftorat(Lcm);
if ( !Lcm ) Lcm = 1;
Mat[I] = [A,Lcm];
}
return [Mat,B];
}
def pfaffian2(Z,V,Beta,SBeta,MN2,XV)
{
N = length(Z);
D = vector(N);
Mat = vector(N);
DV = vector(N);
for ( I = 0; I < N; I++ )
DV[I] = strtov("d"+rtostr(V[I]));
DV = vtol(DV);
for ( I = 0; I < N; I++ ) {
ZI = Z[I];
HI = ZI[0][1];
DI = dp_dtop(HI,PF_DV);
for ( J = 0; J < N; J++ )
if ( var(DI) == DV[J] )
break;
D[J] = deg(DI,DV[J]);
}
B = gen_basis(D,DV);
Dim = length(B);
Hist = [];
for ( I = 0; I < N; I++ ) {
if ( Print ) print(["I=",I]);
A = matrix(Dim,Dim);
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
A[K][L] = [];
for ( J = 0; J < Dim; J++ ) {
if ( Print ) print(["J=",J],2);
Mon = DV[I]*B[J];
for ( T = Hist; T != []; T = cdr(T) )
if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
if ( (T != []) ) {
for ( L = 0; L < N; L++ )
if ( Q == DV[L] ) break;
Red1 = muldpf(V[L],car(T)[1]);
Red = nfpf1(Red1,Z);
} else
Red = nfpf1([mtot(Mon,1)],Z);
Hist = cons([Mon,Red],Hist);
for ( T = Red; T != []; T = cdr(T) ) {
T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
for ( K = 0; K < Dim; K++ )
if ( B[K] == Op ) break;
if ( K == Dim )
error("afo");
else
A[J][K] = C;
}
}
for ( K = 0; K < N; K++ )
A = subst(A,V[K],Beta[K]*XV);
A = map(ctord,A);
Lcm = [];
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
Lcm = lcmf(Lcm,A[K][L][1]);
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ ) {
Num = A[K][L][0]; Den = A[K][L][1];
A[K][L] = Num*ftorat(divf(Lcm,Den));
}
Lcm = ftorat(Lcm);
if ( !Lcm ) Lcm = 1;
A = map(red,A/Lcm);
Lcm = 1;
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
Lcm = lcm(dn(A[K][L]),Lcm);
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ )
A[K][L] = nm(A[K][L])*sdiv(Lcm,dn(A[K][L]));
Mat[I] = [A,Lcm];
}
Lcm = XV;
for ( I = 0; I < N; I++ )
Lcm = lcm(Mat[I][1],Lcm);
/* R = (MN2/XV-SBeta)*Id = (MN2-SBeta*XV)*(Lcm/XV)/Lcm*Id */
R = matrix(Dim,Dim);
D = (MN2-SBeta*XV)*sdiv(Lcm,XV);
for ( I = 0; I < Dim; I++ ) R[I][I] = D;
for ( I = 0; I < N; I++ ) {
A = Mat[I][0];
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ ) {
R[K][L] += Beta[I]*(A[K][L])*sdiv(Lcm,Mat[I][1]);
}
}
Deg = 0;
for ( K = 0; K < Dim; K++ )
for ( L = 0; L < Dim; L++ ) {
DegT = deg(R[K][L],XV);
if ( DegT > Deg ) Deg = DegT;
}
RA = vector(Deg+1);
for ( I = 0; I <= Deg; I++ )
RA[I] = map(coef,R,I);
return [[RA,Lcm],B];
}
def evalpf(F,V,X)
{
if ( !F ) return 0;
F0 = F;
N = length(V);
for ( I = 0; I < N; I++ )
F0 = subst(F0,V[I],X[I]);
return F0;
}
def eval_pfaffian(Mat,Beta,SBeta,MN2,V,XI,F)
{
N = length(V);
Dim = size(Mat[0][0])[0];
R = vector(Dim);
Mul = vector(N);
X = XI*Beta;
X = vtol(X);
for ( K = 0; K < N; K++ )
Mul[K] = Beta[K]/evalpf(Mat[K][1],V,X);
for ( K = 0; K < N; K++ ) {
MatK = Mat[K][0];
for ( I = 0; I < Dim; I++ ) {
for ( U = J = 0; J < Dim; J++ )
U += substr2np(MatK[I][J],V,X)*F[J];
R[I] += Mul[K]*U;
}
}
for ( I = 0; I < Dim; I++ ) R[I] -= (SBeta-MN2/XI)*F[I];
return R;
}
def eval_pfaffian2(Mat,XI,F)
{
MA = Mat[0];
Den = Mat[1];
if ( T = var(Den) ) Den = subst(Den,T,XI);
Len = length(MA);
Dim = size(MA[0])[0];
R = vector(Dim);
for ( I = Len-1; I >= 0; I-- ) {
R = XI*R+MA[I]*F;
}
R = R/Den;
return R;
}
def genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,RK,Step,X0,Eps)
{
if ( type(MapleOut=getopt(mapleout)) == -1 ) MapleOut = 0;
ctrl("bigfloat",1);
O = eval(exp(0));
Len = length(V);
DV = vector(Len);
for ( I = 0; I < Len; I++ )
DV[I] = strtov("d"+rtostr(V[I]));
DV = vtol(DV);
A = (1+M)/2; C = (1+M+N)/2; MN2 = M*N/2;
Z0 = subst(Z,a,A,c,C);
T0 = time();
Q = pfaffian2(Z0,V,Beta,SBeta,MN2,t); Mat = Q[0]; Base = Q[1];
MatLen = length(Q[0][0]);
for ( I = 0; I < MatLen; I++ )
Mat[0][I] *= O;
T1 = time(); Tp = (T1[0]-T0[0])+(T1[1]-T0[1]);
T0 = time();
Beta = ltov(Beta)*O;
X *= O;
Beta0 = Beta*X0;
L = psvalue(Z0,V,Beta0); PS = L[0];
if ( Print ) print(["x0=",X0]);
T1 = time(); Tps = (T1[0]-T0[0])+(T1[1]-T0[1]);
T0 = time();
Dim = length(Base);
ExpF = eval(X0^(MN2)*exp(-SBeta*X0));
GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2)*genprc(M,N,PL,SL);
ExpF *= GM;
F0 = vector(Dim);
for ( I = 0; I < Dim; I++ ) {
DPS = act_op(dp_ptod(Base[I],DV),V,PS);
for ( J = 0; J < Len; J++ )
DPS = subst(DPS,V[J],Beta0[J]);
F0[I] = DPS*ExpF;
}
First = 1;
Val = [];
LVal = [];
DVal = [];
RDVal = [];
if ( MapleOut ) {
mapleout(Mat[0],Mat[1],X0,F0,X,MapleOut);
return;
}
for ( I = 0; ; Step *= 2, I++ ) {
if ( Print ) print("Step=",0); if ( Print ) print(Step);
F = F0*1;
T = F[0]; F /= T;
R = [[X0,T]];
R1 = rk_ratmat(RK,Mat[0],Mat[1],X0,X,Step,F);
R = append(R1,R);
setround(1);
Adj = eval(exp(0));
for ( T = R; T != []; T = cdr(T) )
Adj *= car(T)[1];
LVal = cons(Adj,LVal);
setround(2);
Adj = eval(exp(0));
for ( T = R; T != []; T = cdr(T) )
Adj *= car(T)[1];
Val = cons(Adj,Val);
setround(0);
if ( Print ) { print(""); print([LVal[0],Val[0]]); }
if ( I >= 1 ) DVal = cons(Val[0]-Val[1],DVal);
if ( I >= 2 ) {
RDVal = cons(DVal[1]/DVal[0],RDVal);
print(["log ratio=",eval(log(RDVal[0])/log(2))]);
}
if ( I >= 1 && abs(1-Val[0]/Val[1])<Eps ) {
T1 = time(); Trk = (T1[0]-T0[0])+(T1[1]-T0[1]);
return Val[0];
}
if ( Print ) print("");
}
}
def solve_leq(L,V)
{
#if 0
N = length(V);
B = [];
for ( I = 0; I < N; I++, L = cdr(L) ) B = cons(car(L),B);
while ( 1 ) {
Sol = nd_f4(B,V,0,0);
if ( length(Sol) != N ) {
B = cons(car(L),Sol); L = cdr(L);
} else {
return Sol;
}
}
#else
Sol = nd_f4_trace(L,V,0,-1,0);
return Sol;
#endif
}
def ptom(P,V)
{
M = length(P); N = length(V);
A = matrix(M,N+1);
for ( I = 0; I < M; I++ ) {
F = ptozp(P[I]);
VT = V;
J = 0;
while ( F )
if ( type(F) <= 1 ) {
A[I][N] = F; break;
} else {
for ( FV = var(F); FV != car(VT); VT = cdr(VT), J++ );
A[I][J] = coef(F,1);
F -= A[I][J]*FV;
}
}
return A;
}
def vton(V1,V)
{
N = length(V);
for ( I = 0; I < N; I++ )
if ( V1 == V[I] ) return I;
error("vton : no such variable");
}
def ftotex(F,V)
{
if ( F == [] ) return "1";
R = "";
for ( T = F; T != []; T = cdr(T) ) {
Fac = car(T)[0]; M = car(T)[1];
V1 = var(Fac); NV1 = vton(V1,V);
if ( Fac == V1 ) SFac = "y_{"+rtostr(NV1)+"}";
else {
V2 = var(Fac-V1);
NV2 = vton(V2,V);
SFac = "(y_{"+rtostr(NV1)+"}-y_{"+rtostr(NV2)+"})";
}
if ( M == 1 ) R += SFac;
else R += SFac+"^{"+rtostr(M)+"}";
}
return R;
}
def ctotex(C,V)
{
R = "";
for ( T = C; T != []; T = cdr(T) ) {
if ( T != C ) R += "+";
N = quotetotex(objtoquote(car(T)[0]));
if ( car(T)[1] == [] ) {
R += "("+N+")";
} else {
D = ftotex(car(T)[1],V);
R += "\\frac{"+N+"}{"+D+"}";
}
}
return R;
}
def optotex(Op,DV)
{
E = dp_etov(Op);
N = length(E);
R = "";
for ( I = 0; I < N; I++ )
if ( E[I] )
if ( E[I] == 1 )
R += "\\partial_{"+rtostr(I)+"}";
else
R += "\\partial_{"+rtostr(I)+"}^{"+rtostr(E[I])+"}";
return R;
}
def pftotex(P,V,DV)
{
R = "";
for ( T = P; T != []; T = cdr(T) ) {
if ( T != P ) R += "+";
C = ctotex(car(T)[0],V); Op = optotex(car(T)[1],DV);
R += "("+C+")"+Op+"\\\\";
}
return R;
}
def eqtotex(Z)
{
shell("rm -f __tmp__.tex");
output("__tmp__.tex");
N = length(Z);
print("\\documentclass[12pt]{jsarticle}");
print("\\begin{document}");
print("\\large\\noindent");
for ( I = 0; I < N; I++ ) {
print("$h_"+rtostr(I)+"=",0);
T = mulcpf([[-1,[]]],Z[I]);
print(pftotex(T,PF_V,PF_DV),0);
print("$\\\\");
}
print("\\end{document}");
output();
shell("latex __tmp__");
shell("xdvi __tmp__");
}
/*
* for a partition L=[N1,...,Nl]
* compute the coefficient of x1^N1*...*xl^Nl in phi((x1+...+xM)^N)
* where phi(x(s(1))^n1*...*x(s(k))^nk)=x1^n1*...*xk^nk (n1>=...>=nk)
* if count_dup(L)=[I1,...,Is], then the coefficient is
* C=N!/(N1!*...*Nl!)*M!/((M-l)!*I1!*...*Is!)
* return C*<<0,...0,N1,...,Nl,0,...,0>>
position 0 Base
*/
def coef_partition(L,N,M,Base)
{
R = fac(N);
for ( T = L; T != []; T = cdr(T) )
R = sdiv(R,fac(car(T)));
S = count_dup(L);
R *= fac(M)/fac(M-length(L));
for ( T = S; T != []; T = cdr(T) )
R = sdiv(R,fac(car(T)));
E = vector(length(PF_DV));
for ( I = Base, T = L; T != []; T = cdr(T), I++ )
E[I] = car(T);
return R*dp_vtoe(E);
}
/*
* for a partition L=[N1,...,Nk], return [I1,...,Is]
* s.t. N1=...=N(I1)>N(I1+1)=...=N(I1+I2)>...
*/
def count_dup(L)
{
D = [];
T = L;
while ( T != [] ) {
T0 = car(T); T = cdr(T);
for ( I = 1; T != [] && car(T) == T0; T = cdr(T) ) I++;
D = cons(I,D);
}
return D;
}
/* generate (N_1,...,N_L) s.t. N_1+...+N_L=N, N_1>=...>=N_L>=1, L<= K */
def all_partition(N,K)
{
R = [];
for ( I = 1; I <= K; I++ )
R = append(partition(N,I),R);
return map(reverse,R);
}
/* generate (N_1,...,N_K) s.t. N_1+...+N_K=N, 1<=N_1<=...<=N_K */
def partition(N,K)
{
if ( N < K ) return [];
else if ( N == K ) {
S = [];
for ( I = 0; I < K; I++ )
S = cons(1,S);
return [S];
} else if ( K == 0 )
return [];
else {
R = [];
L = partition(N-K,K);
for ( T = L; T != []; T = cdr(T) ) {
S = [];
for ( U = car(T); U != []; U = cdr(U) )
S = cons(car(U)+1,S);
R = cons(reverse(S),R);
}
L = partition(N-1,K-1);
for ( T = L; T != []; T = cdr(T) )
R = cons(cons(1,car(T)),R);
return R;
}
}
def mul_lopitalpf2(P,I,N,K,I0,I1,E)
{
YI = PF_V[I]; DYI = PF_DV[I];
R = [mtot(DYI^2,1)];
for ( J = I0; J <= I1; J++ ) {
if ( J == I ) continue;
YJ = PF_V[J]; DYJ = PF_DV[J];
R = addpf([mtot(1/2*DYI,YI-YJ)],R);
R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
}
S = subpf(P[I],R);
R = loph(E,I,I0,I1);
if ( type(getopt(high)) == -1 )
S = mul_lopitalpf(S,N,K,I0,I1,E)[1];
else
S = 0;
return [R,S];
}
/* the highest part of lim_{Y(I+1),...,Y(I+N-1)->YI} E(PI) */
/* where E1 = E+(0,...,2,...,0) */
def loph(E,I,I0,I1)
{
E1 = E*1;
E1[I] += 2;
R = [[[[1,[]]],dp_vtoe(E1)]];
S = lopn(E,I,I0,I1);
S = mulcpf([[1/2,[]]],S);
R = addpf(R,S);
R = simppf(R,I0,I1-I0+1);
return R;
}
/* lim_{Y(I+1),...,Y(I+N-1)->YI} dy^E(1/(YI-Y(I+1))+...+1/(YI-Y(I+N-1))) */
def lopn(E,I,I0,I1)
{
R = [];
for ( K = I0; K <= I1; K++ ) {
if ( K != I ) {
L = lop1(E,I,K);
R = addpf(R,L);
}
}
return R;
}
/* lim_{YJ->YI} dy^E(1/(YI-YJ) */
def lop1(E,I,J)
{
YI = PF_V[I]; YJ = PF_V[J];
DYI = PF_DV[I]; DYJ = PF_DV[J];
R = [];
R = addpf([mtot(DYI,YI-YJ)],R);
R = addpf([mtot(-DYJ,YI-YJ)],R);
N = length(PF_V);
BI = E[I]; BJ = E[J];
for ( K = 1, D = 1; K < N; K++ )
if ( K != I && K != J )
D *= PF_DV[K]^E[K];
S = [mtot(-1/(BJ+1)*DYI^(BI+1)*DYJ^(BJ+1)*D,1)];
for ( K = 0; K <= BI; K++ )
if ( -BI+BJ+2*K+2 )
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);
return S;
}
def getchi(N)
{
CHI=[
[5 ,11.07049769,1.145476226],
[10 ,18.30703805,3.940299136],
[20 ,31.41043284,10.85081139],
[50 ,67.50480655,34.76425168],
[100,124.3421134,77.92946517],
[500,553.1268089,449.1467564],
[1000, 1074.679449,927.594363]];
for ( T = CHI; T != []; T = cdr(T) )
if ( car(T)[0] == N ) return car(T);
return [];
}
def message(S) {
dp_gr_print(S);
Print = S;
}
def mapleout(Nm,Dn,X0,F0,X,OutF)
{
N = length(F0);
V = vector(N);
VT = vector(N);
M = length(Nm);
R = matrix(N,N);
for ( I = 0; I < M; I++ )
R += Nm[I]*t^I;
for ( I = 0; I < N; I++ ) {
V[I] = strtov("x"+rtostr(I));
VT[I] = strtov(rtostr(V[I])+"(t)");
}
R *= VT;
output(OutF);
print("dsys:={");
for ( I = 0; I < N; I++ ) {
print("diff("+rtostr(V[I])+"(t),t)=("+rtostr(R[I])+")/("+rtostr(Dn)+"),");
}
for ( I = 0; I < N; I++ ) {
print(rtostr(V[I])+"("+rtostr(X0)+")="+rtostr(F0[I]),0);
if ( I < N-1 ) print(",");
}
print("};");
print("dsol:=dsolve(dsys,numeric,output=Array(["+rtostr(X)+"]));");
output();
}
def gbcheck(Z)
{
N = length(Z);
for ( I = 0; I < N; I++ )
for ( J = I+1; J < N; J++ ) {
S = spolypf(Z[I],Z[J]);
R = nfpf_zero(S,Z);
if ( R != [] ) return 0;
else print([I,J],2);
}
print("");
return 1;
}
def all_partition(N,K)
{
R = [];
for ( I = 1; I <= K; I++ )
R = append(partition(N,I),R);
return map(reverse,R);
}
def partition(N,K)
{
if ( N < K ) return [];
else if ( N == K ) {
S = [];
for ( I = 0; I < K; I++ )
S = cons(1,S);
return [S];
} else if ( K == 0 )
return [];
else {
R = [];
L = partition(N-K,K);
for ( T = L; T != []; T = cdr(T) ) {
S = [];
for ( U = car(T); U != []; U = cdr(U) )
S = cons(car(U)+1,S);
R = cons(reverse(S),R);
}
L = partition(N-1,K-1);
for ( T = L; T != []; T = cdr(T) )
R = cons(cons(1,car(T)),R);
return R;
}
}
def partition_to_pattern(L)
{
Cur = 1;
R = [];
for ( T = L; T != []; T = cdr(T) ) {
R = cons([Cur,Cur+car(T)-1],R);
Cur += car(T);
}
return reverse(R);
}
localf ctond$
def ctond(C)
{
D = [];
for ( T = C; T != []; T = cdr(T) )
D = lcmf(D,car(T)[1]);
N = [];
for ( T = C; T != []; T = cdr(T) ) {
D0 = divf(D,car(T)[1]);
N = cons([car(T)[0],D0],N);
}
return [reverse(N),D];
}
localf tdegf$
def tdegf(F)
{
D = 0;
for ( L = F[1]; L != []; L = cdr(L) )
D += car(L)[1];
return D;
}
localf ntop$
def ntop(N)
{
R = 0;
for ( ; N != []; N = cdr(N) )
R += car(N)[0]*ftop(car(N)[1]);
return R;
}
#if 0
def pwrtab(D)
{
Len = length(D);
F = vector(Len);
M = vector(Len);
for ( I = 0; I < Len; I++ ) {
F[I] = D[I][0]; M[I] = D[I][1];
}
P = vector(Len);
for ( I = 0; I < Len; I++ ) {
MI = M[I]; FI = F[I];
PI = P[I] = vector(MI+1);
for ( PI[0] = 1, J = 1; J <= MI; J++ ) PI[J] = PI[J-1]*FI;
}
}
def ctord2(C)
{
D = [];
for ( T = C; T != []; T = cdr(T) )
D = lcmf(D,car(T)[1]);
H = pwrtab(D);
N = 0;
for ( I = 0, T = reverse(C); T != []; T = cdr(T),I++ ) {
print([Len-I],2);
T0 = car(T); N0 = T0[0]; D0 = T0[1];
L = lcmf(D,D0);
/* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
N = N*getpwr(divf(L,D))+N0*getpwr(divf(L,D0));
D = L;
}
R = [];
for ( T = D; T != []; T = cdr(T) ) {
T0 = car(T); F = T0[0]; M = T0[1];
while ( M > 0 && (Q = tdiv(N,F)) ) {
N = Q;
M--;
}
if ( M ) R = cons([F,M],R);
}
D = reverse(R);
return [N,D];
}
#endif
def degpf(P)
{
R = [];
for ( T = P; T != []; T = cdr(T) ) {
Term = car(T);
C = Term[0]; D = Term[1];
C0 = C[0];
CN = C0[0]; CD = C0[1];
R = cons([CD!=[]?CD[0][1]:0,D],R);
}
return reverse(R);
}
def pftodpf(F)
{
D0 = car(F)[1];
N = length(dp_etov(D0));
DY = [];
for ( I = 0; I < N; I++ )
DY = cons(strtov("dy"+rtostr(I)),DY);
DY = reverse(DY);
R = [];
for ( T = mulcpf([[-1,[]]],F); T != []; T = cdr(T) ) {
C = car(T)[0]; D = dp_dtop(car(T)[1],DY);
R = cons([C,D],R);
}
return reverse(R);
}
def all_diag(N,Dir)
{
P = all_partition(N,N);
U = map(partition_to_pattern,P);
Len = length(P);
for ( I = 0; I < Len; I++ ) {
Z = map(pftodpf,diagpf(N,U[I]));
for ( Name = rtostr(N)+"-", T = P[I]; T != []; T = cdr(T) ) {
Name += rtostr(car(T));
if ( length(T) > 1 ) Name += "_";
}
bsave(Z,Dir+"/"+Name);
}
}
def y1y2(M)
{
wsetup(M);
U = [[2*c-M,[[y1,1]]],[-2,[]]];
for ( J = 3; J <= M; J++ )
U = addc(U,[[1,[[PF_V[1]-PF_V[J],1]]]]);
P = [];
for ( J = 3; J <= M; J++ ) {
C = addc([[1,[[PF_V[1]-PF_V[J],1]]]],[[-1,[[PF_V[1],1]]]]);
P = addpf(P,[[C,dp_ptod(PF_DV[J],PF_DV)]]);
}
}
endmodule$
end$
"usage : n_wishartd.diagpf(M,[[S1,E1],[S2,E2],...]),\n where M is the number of variables and [Si,Ei] is a diagonal block and S(i+1)=Ei + 1. For example, n_wishartd.diagpf(10,[[1,9],[10,10]) returns a system of PDEs satisfied by 1F1(y1,...,y1,y10)."
"usage : n_wishartd.prob_by_hgm(M,N,[P1,P2,...],[S1,S2,...]|options) where M is the number of variables, N is the degrees of freedom, Pi is the size of i-th block and Si is the value of i-th (repeated) eigenvalue of Sigma. options : rk=4|5 => use 4|5-th order Runge-Kutta method (default : rk=5) step=k => set the number of steps=k (default : 10^4) init=x => set the maximal coordinate of the initial point=x (default : 1), eps=e => set the approximate relative error bound=e (default : 10^(-4)). For example, n_wishartd.prob_by_hgm(3,10,[2,1],[1/10,1],10|rk=5) computes Pr[l1<10|diag(1/10,1/10,1)] by RK5."