File: [local] / OpenXM / src / hgm / so3 / src / sord.rr (download)
Revision 1.1, Fri Feb 8 02:59:42 2013 UTC (11 years, 4 months ago) by takayama
Branch: MAIN
CVS Tags: RELEASE_1_3_1_13b, HEAD
Added notes to develope packages.
Global variables are initialized in so3_main.
Added license files.
sord.rr is used to generate a table in so3_nc.c
|
/* sord.rd (for SO, Ohara project).
This is based on wishart/Prog/rd.rr.
cf. h-mle/SO/Notes/note-so3.tex
*/
import("ok_diff.rr")$
/*
p(yi) = yi
q(yi,yj) = 1/(yi^2-yj^2)
*/
dp_ord(0)$
function p(x)$
function q(x,y)$
def pp(Yi) { return 1/Yi; }
def qq(Yi,Yj) { return 1/(Yi^2-Yj^2); }
/* (modified) Muirhead operators */
def mh(I,M) {
E=newvect(M); E[I] = 2;
L = dp_vtoe(E);
Ei=newvect(M); Ei[I]=1;
L += (p-M)*p(util_v(y,[I]))*dp_vtoe(Ei);
for (J=0; J<M; J++) {
Ej=newvect(M); Ej[J]=1;
if (J != I)
L += q(util_v(y,[I]),util_v(y,[J]))*
(util_v(y,[I])*dp_vtoe(Ei)-(util_v(y,[J])*dp_vtoe(Ej)));
}
L -= dp_vtoe(newvect(M));
return(L);
}
def getRuleq(M,N) {
R = [];
for (I=0; I<M; I++) {
for (J=0; J<M; J++) {
if (I!=J) {
Yi = util_v(y,[I]);
Yj = util_v(y,[J]);
for (Ki=0; Ki<=N; Ki++) {
for (Kj=0; Kj<=N; Kj++) {
R = cons([diffn(q(Yi,Yj),[Ki,Kj],[Yi,Yj]),
red(diffn(qq(Yi,Yj),[Ki,Kj],[Yi,Yj]))],R);
}
}
}
}
}
return(R);
}
def getRuleq2(M,N) {
R = [];
for (I=0; I<M; I++) {
for (J=0; J<M; J++) {
if(I!=J) {
Yi = util_v(y,[I]);
Yj = util_v(y,[J]);
for (Ki=0; Ki<=N; Ki++) {
for (Kj=0; Kj<=N; Kj++) {
R = cons([diffn(q(Yi,Yj),[Ki,Kj],[Yi,Yj]), /* Ki,Kj: order of differentiation */
util_v(qq,[I,J,Ki,Kj])],R);
}
}
}
}
}
return(R);
}
def getRulep(M,N) {
R = [];
for (I=0; I<M; I++) {
Yi = util_v(y,[I]);
for (Ki=0; Ki<=N; Ki++) {
R = cons([diffn(p(Yi),[Ki],[Yi]),
red(diffn(pp(Yi),[Ki],[Yi]))],R);
}
}
return(R);
}
def getRulep2(M,N) {
R = [];
for (I=0; I<M; I++) {
Yi = util_v(y,[I]);
for (Ki=0; Ki<=N; Ki++) {
R = cons([diffn(p(Yi),[Ki],[Yi]),
util_v(pp,[I,Ki])],R);
}
}
return(R);
}
/* K! */
def mfac(K) {
N=length(K);
A=1;
for (I=0; I<N; I++) A *= fac(K[I]);
return(A);
}
def mfac2(E,K) {
N=length(K);
A=1;
for (I=0; I<N; I++)
for (J=0; J<K[I]; J++) {
A *= E[I]-J;
}
return(A);
}
/* diff(y^e,k). Returns distributed polynom. */
def mdiff(E,K) {
C = mfac2(E,K);
return(C*dp_vtoe(E-K));
}
def diffn(F,K,V) {
N = length(V);
for (I=0; I<N; I++) {
for (J=0; J<K[I]; J++) {
F = diff(F,V[I]);
}
}
return(F);
}
/* K is a vector, G is a dist poly. V is a var list
*/
def ddiff(G,K,V) {
if (G == 0) return(0);
A = 0;
while (G != 0) {
A += diffn(dp_hc(G),K,V)*dp_ht(G);
G = dp_rest(G);
}
return(A);
}
/* F is a monomial, G is a polynomial. note 2011.10.22
for debug --> test1();
*/
def mpmul(F,G,V) {
if (type(F)<=3) return(F*G);
if (dp_rest(F) != 0) error("mpmul(F,G,V), F must be a monomial.");
N = length(V);
E = dp_etov(F); C = dp_hc(F);
P=1;
for (I=0; I<N; I++) {
P *= (1+V[I])^E[I];
}
if (type(V) == 5) V=vtol(V);
P = dp_ptod(P,V);
A = 0;
while (P!= 0) {
K = dp_etov(dp_hm(P)); P = dp_rest(P);
A += (1/mfac(K))*mdiff(E,K)*ddiff(G,K,V);
}
return(C*A);
}
def test1() {
F=x*dx^2*dy;
G=x*y*dx^2+(x-y)*dx*dy+(1+x+y)^2;
A=sm1.mul(F,G,[x,y]);
V=[dx,dy];
A1=mpmul(dp_ptod(F,V),dp_ptod(G,V),[x,y]);
A1=dp_dtop(A1,V);
return(A-A1);
}
def dp_ptod0(F,Type) {
A = F+dp_hm(Type);
A = A-dp_hm(Type);
return(A);
}
/* reduce F by G */
def red1(F,G,V) {
if (F == 0) return([F,0]);
A = 0;
Q = 0;
while (F != 0) {
if (dp_redble(F,G)) {
Q0 = dp_ptod0((dp_hc(F)/dp_hc(G)),F)*dp_subd(F,G);
R = mpmul(-Q0,dp_rest(G),V);
/* R = -Q0*dp_rest(G); */
A += R; Q += Q0;
}else{
A += dp_hm(F);
}
F = dp_rest(F);
}
return([A,Q]);
}
/*
This returns error because coef of dist poly must not be rational.
*/
def test2() {
F = a*<<2,1>> + b*<<1,1>> + p*<<0,0>>;
G = c*<<1,1>> - d*<<0,0>>;
return(red1(F,G,[x,y]));
}
def test3() {
F0=(a/b)+<<1,0>>;
F =F0-<<1,0>>;
print(F);
G = c*<<1,0>>;
H = F*G;
print(type(dp_hm(H)));
return(rtostr(H)); /* ask noro 2011.10.23 */
}
def test2b() {
F = a*<<2,1>> + b*<<1,1>> + p*<<0,0>>;
G = (1/2)*<<1,1>> - d*<<0,0>>;
A=red1(F,G,[x,y]);
print(F-A[1]*G-A[0]);
return(A);
}
/* G is a set */
def redall0(F,G,V) {
N = length(G);
Q = newvect(N);
Reducible = 1;
while (Reducible) {
Reducible = 0;
for (I=0; I<N; I++) {
if (dp_redble(F,G[I])) {
A = red1(F,G[I],V);
Q[I] += A[1];
F = A[0];
if (F == 0) return([F,Q]);
Reducible = 1;
}
}
}
return([F,Q]);
}
def redall(F,G,V) {
Q = newvect(length(V));
if (F == 0) return [F,Q];
R = 0;
while (F != 0) {
A = redall0(F,G,V);
Q = Q+A[1];
R += dp_hm(A[0]);
F = dp_rest(A[0]);
}
return([R,Q]);
}
def test4() {
F = a*<<1,2>>+b*<<2,0>>;
G = [(1/2)*<<1,1>>+c*<<0,0>>,(1/3)*<<2,0>>+d*<<0,0>>];
A=redall(F,G,[x,y]);
print(F-A[1][0]*G[0]-A[1][1]*G[1]-A[0]);
return(A);
}
def test5() {
dp_ord(0);
M=2;
G = [mh(0,M),mh(1,M)];
A=redall(<<2,1>>,G,[y_0,y_1]);
return(A[0]);
}
/* red for dp */
def dpred0(F) {
A = 0;
while (F!=0) {
A += red(dp_hc(F))*dp_ht(F); /* asir does not accept */
F = dp_rest(F);
}
return(A);
}
def test6() {
M=2;
G = [mh(0,M),mh(1,M)];
V = [y_0,y_1];
S=mpmul(<<0,2>>,G[0],V)-mpmul(<<2,0>>,G[1],V);
A=redall(S,G,V);
Rule = append(getRulep(2,2),getRuleq(2,2));
/* A = dpred(base_replace(A[0],Rule)); */
A = red(dp_dtop(base_replace(A[0],Rule),[dy_1,dy_2]));
return(A);
}
def myelim(L,V) {
F=1;
for (I=0; I<length(V); I++) F = F*V[I];
V = vars(V);
A = [];
for (I=0; I<length(L); I++) {
if (vars(L[I]) == V) A=cons(L[I],A);
}
return(reverse(A));
}
/* rel_o(2) is fine, but others return no relations
2011.11.24 at Kanazawa
*/
def rel_o(N) {
T = newmat(N,N);
for (I=0; I<N; I++) for (J=0; J<N; J++) T[I][J] = util_v(t,[I,J]);
V0 = []; V1 = [];
for (I=0; I<N; I++) V0 = cons(util_v(t,[I,I]),V0);
for (I=0; I<N; I++) for (J=0; J<N; J++)
if (I != J) V1 = cons(util_v(t,[I,J]),V1);
V = append(V1,V0);
print(V);
T2 = matrix_transpose(T)*T - matrix_identity_matrix(N);
L = base_flatten(T2);
/* L = cons(matrix_det(T)-1,L); */
return myelim(nd_gr(L,V,0,2),V0);
}
/* Evaluation of normalization constant */
def fac2(N) {
if (N < -1) error("not implemented.");
if (N == -1) return 1; /* cf. google: double factorial */
if (N == 0) return 1;
return N*fac2(N-2);
}
/* cf. etrcalc.r */
extern S_dm$
extern S_ZR$
S_dm = 20$
S_ZR = 2$
extern S_fct$
extern S_dfct$
extern S_Etens$
S_fct=newvect(S_ZR+S_dm+1)$
S_dfct=newvect(S_ZR+4*S_dm+1)$
S_Etens=newvect(S_ZR+S_dm+1)$
for (I=0; I<S_ZR+S_dm+1; I++) S_Etens[I] = newmat(S_ZR+S_dm+1,S_ZR+S_dm+1)$
def comb(N,K) {
return fac(N)/(fac(K)*fac(N-K));
}
def init_table() {
S_fct[S_ZR] = 1;
for (I=1; I<=S_dm; I++) S_fct[S_ZR+I] = fac(I);
S_dfct[S_ZR-1] = S_dfct[S_ZR] = 1;
for (I=1; I<=4*S_dm; I++) S_dfct[S_ZR+I] = fac2(I);
for (K=0; K<=S_dm; K++) {
for (L=0; L<=S_dm; L++) {
for (M=0; M<=S_dm; M++) {
if (((K-L) % 2 == 0) && ((L-M) % 2 == 0)) {
E=0;
for (I=0; I<=idiv(L,2); I++) { /* L-N=2*I */
A0 = comb(L,2*I);
A1 = S_dfct[S_ZR+K+M]*S_dfct[S_ZR+2*I-1]/S_dfct[S_ZR+K+M+2*I+1];
A2 = S_dfct[S_ZR+K+L-2*I-1]*S_dfct[S_ZR+2*I-1]/S_dfct[S_ZR+K+L];
A3 = S_dfct[S_ZR+M+L-2*I-1]*S_dfct[S_ZR+2*I-1]/S_dfct[S_ZR+M+L];
E = E + A0*A1*A2*A3;
}
(S_Etens[S_ZR+K])[S_ZR+L][S_ZR+M] = E;
}else{
(S_Etens[S_ZR+K])[S_ZR+L][S_ZR+M] = 0;
}
}}}
}
def so_int_d(N) {
F=0;
if (S_fct[S_ZR] != 1) {print("Initializing table."); init_table();}
if (N > S_dm) error("Approximation degree is too big.");
Kmax = Lmax = Mmax = N-1;
for (K=0; K<=Kmax; K++)for (L=0; L<=Lmax; L++)for (M=0; M<=Mmax; M++) {{{
F = F+dp_vtoe(newvect(3,[K,L,M]))*
(S_Etens[S_ZR+K][S_ZR+L][S_ZR+M]/(S_fct[S_ZR+K]*S_fct[S_ZR+L]*S_fct[S_ZR+M]));
}}}
return F;
}
def degmin(L) {
L = dp_ptod(L,vars(L));
D = 1000; /* big */
while (L != 0) {
T = dp_etov(L);
if (T[0]+T[1]+T[2] < D) D = T[0]+T[1]+T[2];
L = dp_rest(L);
}
return(D);
}
def check1(P) {
V = [y_0,y_1,y_2];
DV = [dy_0,dy_1,dy_2];
L0 = mh(0,3);
L1 = mh(1,3);
L2 = mh(2,3);
L = [L0,L1,L2];
R = append(getRulep(3,0),getRuleq(3,0));
R = append(R,[[p,3]]);
L = map(dp_dtop,L,DV);
L = base_replace(L,R);
L = map(nm,L);
print(L);
F = dp_dtop(so_int_d(S_dm-P),V);
A=map(odiff_act,L,F,V); dp_ord(0);
return(A);
}
/*
Use odiff_act(L,F,V), ord(0); cf. wishart/Prog/w1.rr
*/
/* ---- done ---- */
/* dicimal to binary vector of size N */
def dtobv(A,N) {
B = newvect(N);
for (I=0; I<N; I++) {
B[I] = A % 2;
A = idiv(A,2);
}
return(B);
}
def bvtod(B) {
N = length(B);
A = 0;
for (I=0; I<N; I++) A += B[I]*2^I;
return(A);
}
def mhbase(M) {
A = [];
for (I=0; I<2^M; I++) {
A = cons(dp_vtoe(dtobv(I,M)),A);
}
return reverse(A);
}
def mhv(M) {
V = newvect(M);
for (I=0; I<M; I++) V[I] = util_v(y,[I]);
return V;
}
/*
F = mhbase(M);
dF/y_II = P_II F,
Returns P_II
p, q are not reduced.
*/
def mhmat0(M,II) {
V = mhv(M);
B = mhbase(M);
D = newvect(M); D[II] = 1; D = dp_vtoe(D);
DB = newvect(2^M);
for (I=0; I<2^M; I++) DB[I] = mpmul(D,B[I],V);
G = newvect(M);
for (I=0; I<M; I++) G[I] = mh(I,M);
P = newmat(2^M,2^M);
for (I=0; I<2^M; I++) {
F=redall(DB[I],G,V); F=F[0];
while (F!=0) {
J = bvtod(dp_etov(dp_ht(F)));
P[I][J] = dp_hc(F);
F = dp_rest(F);
}
}
return(P);
}
def mhmat(M,II) {
Rule = append(getRuleq(M,1),getRulep(M,1));
return map(red,base_replace(mhmat0(M,II),Rule));
}
def mhmat2(M,II) {
Rule = append(getRuleq2(M,1),getRulep2(M,1));
return base_replace(mhmat0(M,II),Rule);
}
def test7() {
P0 = mhmat(2,0);
P1 = mhmat(2,1);
/* Check the integrability conditions.
dF/dy_0 = P0 F, dF/dy_1 = P1 F
d^2F/dy_1 dy_0 = dP0/dy_1 F + P0 dF/dy_1,
d^2F/dy_0 dy_1 = dP1/dy_0 F + P1 dF/dy_0,
dP0/dy_1 + P0*P1 = dP1/dy_0 + P1*P0
*/
A0= map(diff,P0,y_1) + P0*P1;
A1= map(diff,P1,y_0) + P1*P0;
return map(red,A0-A1); /* it should return 0 */
}
/* It is easy to get up to mhmat0(6,0), but mhmat(5,0) seems to take a time. */
/*
y_0 = S0*t, y_1 = S1*t
eigen2(3,1,2);
*/
def eigen2(N,S0,S1) {
P0 = mhmat(2,0);
P1 = mhmat(2,1);
Rule = [[y_0,S0*t^(-1)],[y_1,S1*t^(-1)],[a,3/2],[c,3/2+N/2]];
P0 = base_replace(P0,Rule);
P1 = base_replace(P1,Rule);
/* print(P0);
print(P1); */
Rule2 = [[t,0]];
P0 = base_replace(P0,Rule2);
P1 = base_replace(P1,Rule2);
D = newmat(4,4);
Dii = -(S0+S1)-s;
for (I=0; I<4; I++) D[I][I] = Dii;
E0 = D + S0*P0+S1*P1;
E1 = det(E0);
if (length(vars(E1)) > 1) return(fctr(E1));
else return(pari(roots,E1));
/*
E=matrix_eigenvalues(E0);
return(E);
*/
}
def eigen3(N,S0,S1,S2) {
P0 = mhmat(3,0);
P1 = mhmat(3,1);
P2 = mhmat(3,2);
Rule = [[y_0,S0*t^(-1)],[y_1,S1*t^(-1)],[y_2,S2*t^(-1)],
[a,4/2],[c,4/2+N/2]];
P0 = base_replace(P0,Rule);
P1 = base_replace(P1,Rule);
P2 = base_replace(P2,Rule);
Rule2 = [[t,0]];
P0 = base_replace(P0,Rule2);
P1 = base_replace(P1,Rule2);
P2 = base_replace(P2,Rule2);
D = newmat(8,8);
Dii = -(S0+S1+S2)-s;
for (I=0; I<8; I++) D[I][I] = Dii;
E0 = D + S0*P0+S1*P1 + S2*P2;
E1=det(E0);
if (length(vars(E1)) > 1) return(fctr(E1));
else return(pari(roots,E1));
}
def getSer(N) {
F=so_int_d(N);
F=dp_dtop(F,[x,y,z]);
G=[F,diff(F,x),diff(F,y),diff(F,z)];
G=base_replace(G,[[x,a*t],[y,b*t],[z,c*t]]);
G=[G[0]*exp(0),G[1]*exp(0),G[2]*exp(0),G[3]*exp(0)];
G=map(eval,G);
G=map(deval,G);
return(G);
}
def getSer2(N) {
F=so_int_d(N);
F=dp_dtop(F,[x,y,z]);
G=[F,diff(F,x),diff(F,y),diff(F,z)];
G=[G[0]*exp(0),G[1]*exp(0),G[2]*exp(0),G[3]*exp(0)];
G=map(eval,G);
G=map(deval,G);
return(G);
}
def foo1(N) {
F=getSer(N);
Rule=[[a,-5.614],[b,3.079],[c,2.387],[t,1]]; /* Commets */
Val=base_replace(F,Rule);
return(Val);
}
def forC(N) {
F=so_int_d(N);
while (F!=0) {
E=dp_etov(dp_ht(F));
C=deval(eval(dp_hc(F)*exp(0)));
printf("Tnc[%a][%a][%a]=%a;\n",E[0],E[1],E[2],C);
F=dp_rest(F);
}
}
end$