import("nk_restriction.rr");
module ns_twistedlog;
localf mkd,mkd_,mke,mke_;
localf dp_dtol,dp_etop,ltop;
localf gauss,excp_gauss,leadInd;
localf t_in,m_in,c_in,reducible;
localf normal_form_weyl,w_normal_form_weyl,sp_weyl,buchberger_weyl;
localf normal_form,w_normal_form,sp,buchberger;
localf v_homo,v_order;
localf adjoint;
localf twisted_log_cohomology,twisted_deRham;
localf difference_equation;
localf differential_equation;
localf monomial_sum,monomial_quotient;
localf combination,homogeneous;
localf monomial_hilbert_poly,hilbelt_function,hilbelt_poly,affine_hilbelt_poly;
localf holonomic;
localf mkt,mky,mkx;
localf example6_5;
localf example6_6;
localf example6_9;
localf example6_10;
localf example7_1;
localf example7_2;
localf example7_4;
localf example7_5;
localf example7_6;
localf time_f1,time_f1_a;
localf time_f2,time_f2_a;
localf time_f3,time_f3_a;
localf time_f4,time_f4_a;
localf time_g;
localf time_h,time_h_a;
def mkd(L){
return strtov("d" + rtostr(L));
}
def mkd_(L,N){
return strtov("d" + rtostr(L) + "_" + rtostr(N));
}
def mke(L){
return strtov("e" + rtostr(L));
}
def mke_(L,N){
return strtov("e" + rtostr(L) + "_" + rtostr(N));
}
def dp_dtol(F,V){
N = length(V);
VF = 0;
while (F != 0) {
E = dp_etov(dp_ht(F));
M = 0;
for (I = 0; I < N; I++)
M += rtostr(V[I]) + rtostr(E[I]);
VF += dp_hc(F) * strtov(M);
F = dp_rest(F);
}
return VF;
}
def dp_etop(F,V,M){
N = length(V);
NM = N*M;
P = 0;
while (F != 0) {
E = dp_etov(dp_ht(F));
if (E[NM] != 0) {
P += dp_hc(F);
F = dp_rest(F);
} else {
for (I = 0; I < NM; I++)
if (E[I] != 0)
break;
CN = idiv(I,M);
CM = I % M;
P += dp_hc(F) * V[CN]^(CM + 1);
F = dp_rest(F);
}
}
return P;
}
def ltop(P,V,M){
N = length(V);
W = newvect(N);
for (I = 0; I < N; I++)
W[I] = 1;
W = vtol(W);
B = nk_restriction.gen(W,M);
B = map(dp_vtoe,map(ltov,B));
NB = length(B);
BX = map(dp_dtop,B,V);
BXX = map(dp_dtol,B,V);
DP = dp_ptod(P,BXX);
PP = 0;
while (DP != 0) {
EP = dp_etov(dp_ht(DP));
for (I = 0; I < NB; I++)
if (EP[I] != 0)
break;
PP += dp_hc(DP) * BX[I];
DP = dp_rest(DP);
}
return PP;
}
def gauss(A)
{
MN = size(A);
M = MN[0];
N = MN[1];
LeadInd = newvect(M);
for( I = 0; I < M; I++ )
LeadInd[I] = leadInd(A[I],N);
for( J = 0; J < N; J++ ){
I0 = -1;
for( I = 0; I < M && I0 == -1; I++ ){
if( LeadInd[I] == J ){
C = A[I][J];
I0 = I;
for( K = J; K < N; K++ )
A[I][K] = red(A[I][K]/C);
}
}
if( I0 >= 0 ){
for( I = 0; I < M; I++ ){
C = A[I][J];
if( I != I0 && C != 0 ){
for( K = J; K < N; K++ )
A[I][K] = red(A[I][K]-C*A[I0][K]);
LeadInd[I] = leadInd(A[I],N);
}
}
}
}
return ([A,LeadInd]);
}
def excp_gauss(A)
{
ExcpList = [];
MN = size(A);
M = MN[0];
N = MN[1];
LeadInd = newvect(M);
for( I = 0; I < M; I++ )
LeadInd[I] = leadInd(A[I],N);
for( J = 0; J < N; J++ ){
I0 = -1;
for( I = 0; I < M && I0 == -1; I++ ){
if( LeadInd[I] == J ){
C = A[I][J];
ExcpList = cons(C,ExcpList);
I0 = I;
for( K = J; K < N; K++ )
A[I][K] = red(A[I][K]/C);
}
}
if( I0 >= 0 ){
for( I = 0; I < M; I++ ){
C = A[I][J];
if( I != I0 && C != 0 ){
for( K = J; K < N; K++ )
A[I][K] = red(A[I][K]-C*A[I0][K]);
LeadInd[I] = leadInd(A[I],N);
}
}
}
}
return ([A,LeadInd,ExcpList]);
}
def t_in(F, VL, Order)
{
OldOrder = dp_ord();
dp_ord(Order);
NM = nm(F);
DNM = dp_ptod(NM,VL);
DIN = dp_hm(DNM);
IN = dp_dtop(DIN,VL);
IN = IN / dn(F);
IN = red(IN);
dp_ord(Order);
return IN;
}
def m_in(F, VL, Order)
{
OldOrder = dp_ord();
dp_ord(Order);
NM = nm(F);
DNM = dp_ptod(NM,VL);
DIN = dp_ht(DNM);
IN = dp_dtop(DIN,VL);
dp_ord(Order);
return IN;
}
def c_in(F, VL, Order)
{
OldOrder = dp_ord();
dp_ord(Order);
NM = nm(F);
DNM = dp_ptod(NM,VL);
LC = dp_hc(DNM);
LC = LC / dn(F);
LC = red(LC);
dp_ord(OldOrder);
return LC;
}
def reducible(R, G, VL, Order)
{
M = length(G);
InR = m_in(R, VL, Order);
for (I = 0; I < M; I++) {
InG = m_in(G[I], VL, Order);
T = tdiv(InR, InG);
if (T != 0)
return I;
}
return -1;
}
def normal_form(F, G, VL, Order)
{
M = length(G);
Q = newvect(M);
R = 0;
LC = [c_in(F, VL, Order)];
while (F != 0) {
L = w_normal_form(F, G, VL, Order);
RR = L[0];
QQ = L[1];
LC = append(LC, L[2]);
F = RR - t_in(RR, VL, Order);
LC = cons(c_in(F, VL, Order), LC);
R += t_in(RR, VL, Order);
Q += QQ;
R = red(R);
Q = map(red, Q);
}
return [R, Q, LC];
}
def w_normal_form(F, G, VL, Order)
{
M = length(G);
Q = newvect(M);
R = F;
LC = [c_in(F, VL, Order)];
while ((Index = reducible(R, G, VL, Order)) != -1) {
S = G[Index];
D = tdiv(m_in(R, VL, Order), m_in(S, VL, Order));
X = c_in(R, VL, Order) / c_in(S, VL, Order);
T = X*D*S;
Q[Index] += X*D;
R -= T;
Q[Index] = red(Q[Index]);
R = red(R);
LC = cons(c_in(F, VL, Order),LC);
}
return [R, Q, LC];
}
def sp(F, G, VL, Order)
{
InF = m_in(F, VL, Order);
InG = m_in(G, VL, Order);
LCF = c_in(F, VL, Order);
LCG = c_in(G, VL, Order);
LCM = lcm(InF, InG);
MF = tdiv(LCM, InF);
MG = tdiv(LCM, InG);
Sp = MF * F - red(LCF/LCG) * MG * G;
return red(Sp);
}
def buchberger(L, VL, Order)
{
LC = [];
N = length(L);
for (I = 0; I < N; I++)
LC = cons(c_in(L[I], VL, Order), LC);
Jlist = [];
for (I = 0; I < N; I++)
Jlist = append(Jlist,[[I, J]]);
while (Jlist != []) {
Index = car(Jlist);
Jlist = cdr(Jlist);
Sp = sp(L[Index[0]], L[Index[1]], VL, Order);
LC = cons(c_in(SP, VL, Order), LC);
Result = normal_form(Sp, L, VL, Order);
R = Result[0];
LC = append(LC, Result[2]);
if (R != 0) {
L = append(L, [R]);
N = length(L);
for (I = 0; I < N-1; I++)
Jlist = append(Jlist, [[I, N-1]]);
}
}
return [L,LC];
}
def normal_form_weyl(F, G, VL, Order)
{
M = length(G);
Q = newvect(M);
R = 0;
LC = [c_in(F, VL, Order)];
while (F != 0) {
L = w_normal_form_weyl(F, G, VL, Order);
RR = L[0];
QQ = L[1];
LC = append(LC, L[2]);
F = RR - t_in(RR, VL, Order);
LC = cons(c_in(F, VL, Order), LC);
R += t_in(RR, VL, Order);
Q += QQ;
R = red(R);
Q = map(red, Q);
}
return [R, Q, LC];
}
def w_normal_form_weyl(F, G, VL, Order)
{
M = length(G);
Q = newvect(M);
R = F;
LC = [c_in(F, VL, Order)];
while ((Index = reducible(R, G, VL, Order)) != -1) {
S = G[Index];
D = tdiv(m_in(R, VL, Order), m_in(S, VL, Order));
X = c_in(R, VL, Order) / c_in(S, VL, Order);
T = X * nk_restriction.weyl_mul(D, nm(S), VL) / dn(S);
Q[Index] += X*D;
R -= T;
Q[Index] = red(Q[Index]);
R = red(R);
LC = cons(c_in(F, VL, Order),LC);
}
return [R, Q, LC];
}
def sp_weyl(F, G, VL, Order)
{
InF = m_in(F, VL, Order);
InG = m_in(G, VL, Order);
LCF = c_in(F, VL, Order);
LCG = c_in(G, VL, Order);
LCM = lcm(InF, InG);
MF = tdiv(LCM, InF);
MG = tdiv(LCM, InG);
Sp = nk_restriction.weyl_mul(MF, nm(F), VL) / dn(F) - red(LCF/LCG) * nk_restriction.weyl_mul(MG, nm(G), VL) / dn(G);
return red(Sp);
}
def buchberger_weyl(L, VL, Order)
{
LC = [];
N = length(L);
for (I = 0; I < N; I++)
LC = cons(c_in(L[I], VL, Order), LC);
Jlist = [];
for (I = 0; I < N; I++)
Jlist = append(Jlist,[[I, J]]);
while (Jlist != []) {
Index = car(Jlist);
Jlist = cdr(Jlist);
Sp = sp_weyl(L[Index[0]], L[Index[1]], VL, Order);
LC = cons(c_in(SP, VL, Order), LC);
Result = normal_form_weyl(Sp, L, VL, Order);
R = Result[0];
LC = append(LC, Result[2]);
if (R != 0) {
L = append(L, [R]);
N = length(L);
for (I = 0; I < N-1; I++)
Jlist = append(Jlist, [[I, N-1]]);
}
}
return [L,LC];
}
def v_homo(F,V,DV,W)
{
N = length(V);
VDV = append(V,DV);
VDVH = append(VDV,[h]);
DF = dp_ptod(F,VDV);
F = dp_ptod(F,VDVH);
WL = [];
MIN = [];
while (DF != 0) {
Index = dp_etov(dp_ht(DF));
WW = 0;
for (I = 0; I < N; I++)
WW += W[I] * Index[I];
for (I = 0; I < N; I++)
WW -= W[I] * Index[I+N];
if (MIN == []) {
MIN = WW;
} else if (WW < MIN) {
MIN = WW;
}
WL = cons(WW,WL);
DF = dp_rest(DF);
}
WL = ltov(reverse(WL));
NL = length(WL);
for (I = 0; I < NL; I++)
WL[I] -= MIN;
HF = 0;
for (I = 0; I < NL; I++) {
LM = dp_etov(dp_ht(F));
LC = dp_hc(F);
LM[2*N] = WL[I];
HF += LC * dp_vtoe(LM);
F = dp_rest(F);
}
return dp_dtop(HF,VDVH);
}
def v_order(N)
{
M = newmat(2*N+3,2*N+2);
M[0][N] = 1;
for (I = 0; I < 2*(N+1); I++)
M[1][I] = 1;
for (I = 0; I < 2*N+1; I++)
M[I+2][2*N+1-I] = -1;
return M;
}
def leadInd(V,M){
Index = -1;
for( I = 0; I < M && Index == -1; I++ )
if( V[I] != 0 )
Index = I;
return(Index);
}
def adjoint(P,VL,DVL)
{
N = length(VL);
VLL = append(VL,DVL);
DP = dp_ptod(P,VLL);
S = 0;
EX = newvect(2*N);
EDX = newvect(2*N);
while(DP != 0) {
M = dp_etov(dp_ht(DP));
C = dp_hc(DP);
Sum = 0;
for( I = 0; I < N; I++ )
EX[I] = M[I];
for( I = 0; I < N; I++ ){
EDX[N+I] = M[N+I];
Sum += M[N+I];
}
S += C * (-1)^Sum * dp_weyl_mul(dp_vtoe(EDX),dp_vtoe(EX));
DP = dp_rest(DP);
}
return dp_dtop(S,VLL);
}
def twisted_log_cohomology(F,P,V)
{
if ( type(Exp=getopt(exp)) == 2) ; else Exp = 0;
if ( type(Check=getopt(check)) == -1 ) Check = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(Excp=getopt(excp)) == -1 ) Excp = -1;
N = length(F);
if( length(P) != N ){
print("error : number of parameters is not equal length of factorization");
return;
}
Vnum = length(V);
DV = map(mkd,V);
FF = 1;
for( I = 0; I < N; I++ )
FF *= F[I];
L = newvect(Vnum + 1);
for (I = 0; I < Vnum; I++)
L[I] = (-1)^(I+1) * diff(FF,V[I]);
L[Vnum] = FF;
L = vtol(L);
L = newsyz.module_syz(L,V,1,0);
L = L[0];
NL = length(L);
DL = [];
for( I = 0; I < NL; I++ ){
M = 0;
for (K = 0; K < Vnum; K++)
M += (-1)^(K+1) * L[I][K] * DV[K];
M -= L[I][Vnum];
for (J = 0; J < N; J++) {
MJ = 0;
for (K = 0; K < Vnum; K++)
MJ += (-1)^K * L[I][K] * diff(F[J],V[K]);
MJ = P[J] * sdiv(MJ,F[J]);
M += MJ;
}
if ( Exp ) {
ME = 0;
for (K = 0; K < Vnum; K++)
ME += (-1)^K * L[I][K] * diff(Exp,V[K]);
M += ME;
}
DL = cons(M,DL);
}
if (Check)
holonomic(DL,V,DV);
FDL = map(nk_restriction.fourier_trans,DL,V,DV);
Param = vars(P);
W = newvect(Vnum);
for (I = 0; I < Vnum; I++)
W[I] = 1;
W = vtol(W);
if (Excp == -1) {
if (S0 == -1) {
if (Param == [])
Result = nk_restriction.restriction(FDL,V,DV,W);
else
Result = nk_restriction.restriction(FDL,V,DV,W | param = Param);
if( Result == 0 )
return 0;
Generator = Result[0];
Gen = Result[1];
} else {
WW = newvect(2 *Vnum);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
VDV = append(V,DV);
VDVH = append(VDV,[h]);
FDLH = map(dp_dtop,map(dp_homo,map(dp_ptod,FDL,VDV)),VDVH);
GBH = nd_weyl_gr(FDLH,VDVH,0,11);
GB = map(subst,GBH,h,1);
NN = length(GB);
VarV = append(V, DV);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VarV, WW);
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Gen = nk_restriction.gen(W,S0);
}
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
Index = (gauss(M))[1];
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
dp_ord(OldOrder);
print("dimension : " + rtostr(length(Base)));
return map(dp_dtop,Base,V);
} else {
NotInt = [];
while (P != []) {
A = car(P);
P = cdr(P);
if (vars(A) != [])
NotInt = cons(A,NotInt);
}
NotInt = reverse(NotInt);
HFDL = map(v_homo,FDL,V,DV,W);
VDV = append(V,DV);
VH = append(V,[h]);
DVH = append(DV,[dh]);
VDVH = append(VH, DVH);
BB = buchberger_weyl(HFDL, VDVH, v_order(Vnum));
GB = map(nm,map(subst,BB[0],h,1));
ExcpList = BB[1];
Vnum2 = 2*Vnum;
dp_weyl_set_weight(W);
M = newmat(Vnum2,Vnum2);
for ( J = 0; J < Vnum2; J++ )
M[0][J] = 1;
for ( I = 1; I < Vnum2; I++ )
M[I][Vnum2-I] = -1;
MW = newmat(Vnum2+1,Vnum2);
for ( J = 0; J < Vnum; J++ )
MW[0][J] = -W[J];
for ( ; J < Vnum2; J++ )
MW[0][J] = W[J-Vnum];
for ( I = 1; I <= Vnum2; I++ )
for ( J = 0; J < Vnum2; J++ )
MW[I][J] = M[I-1][J];
GIN = map(initial_part,GB,VDV,MW,W);
for ( I = 0, T = 0; I < Vnum; I++ )
T += W[I]*V[I]*DV[I];
BF = nk_restriction.weyl_minipoly_by_elim2(GIN,V,DV,Param,T);
L = fctr(BF);
BFF = L;
L = cdr(L);
NL = length(L);
S0 = [];
NonPosInt = [];
for (I = 0; I < NL; I++) {
T = L[I][0];
if ( vars(T) == [s] ) {
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
} else {
NonPosInt = cons(-coef(T, 0, s) / coef(T, 1, s), NonPosInt);
}
}
print("generic bfct : " + rtostr(BFF), 1);
print("S0 : " + rtostr(S0), 1);
WW = newvect(Vnum2);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
NN = length(GB);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VDV, WW);
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VDV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Gen = nk_restriction.gen(W,S0);
print("B_{S0} length : " + rtostr(length(Gen)),1);
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
M = excp_gauss(M);
Index = M[1];
ExcpList = append(M[2],ExcpList);
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
dp_ord(OldOrder);
print("dimension : " + rtostr(length(Base)));
ExcpSet = [];
while (ExcpList != []) {
M = car(ExcpList);
ExcpList = cdr(ExcpList);
if (type(NM = nm(M)) == 2 && !member(NM,ExcpSet) && !member(-NM,ExcpSet))
ExcpSet = cons(nm(M),ExcpSet);
if (type(DN = dn(M)) == 2 && !member(DN,ExcpSet) && !member(-DN,ExcpSet))
ExcpSet = cons(dn(M),ExcpSet);
}
FExcpSet = [];
while (ExcpSet != []) {
ES = car(ExcpSet);
ExcpSet = cdr(ExcpSet);
FExcpSet = append(FExcpSet,fctr(ES));
}
while (FExcpSet != []) {
M = car(FExcpSet)[0];
FExcpSet = cdr(FExcpSet);
if (type(M) == 2 && !member(M,ExcpSet))
ExcpSet = cons(M,ExcpSet);
}
return ["Basis",map(dp_dtop,Base,V),"Not integer",NotInt,"Not non-negative integer",NonPosInt, "Not zero",ExcpSet];
}
}
def twisted_deRham(F,A,V)
{
DV = map(mkd,V);
Vnum = length(V);
L = ann(F);
BF = bfct(F);
BFF = cdr(fctr(BF));
N = length(BFF);
R = newvect(N);
R0 = newvect(N);
for (I = 0; I < N; I++) {
R[I] = -subst(BFF[I][0],s,0)/coef(BFF[I][0],1,s);
R0[I] = 1;
}
R1 = A * R0 - R;
Crt = -1;
for (I = 0; I < N; I++)
if (R1[I] > 0 && dn(R1[I]) == 1) {
Crt = I;
break;
}
if (Crt == -1)
L = map(subst,L,s,A);
else {
LL = map(subst,L,s,V[I]);
LL = cons(F^R1[I],LL);
LL = newsyz.module_syz(LL,append(V,DV),1,0 | weyl = 1);
LL = LL[0];
L = [];
for (; LL != []; LL = cdr(LL))
L = cons((car(LL))[0],L);
}
FL = map(nk_restriction.fourier_trans,L,V,DV);
W = newvect(Vnum);
for (I = 0; I < Vnum; I++)
W[I] = 1;
W = vtol(W);
Result = nk_restriction.restriction(FL,V,DV,W);
if( Result == 0 )
return 0;
Generator = Result[0];
Gen = Result[1];
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OrdOrder = dp_ord();
dp_ord(Mat);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
Index = (gauss(M))[1];
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
dp_ord(OldOrder);
print("dimension : " + rtostr(length(Base)));
return map(dp_dtop,Base,V);
}
def difference_equation(F,P,V)
{
if ( type(IH = getopt(inhomo)) == -1 ) IH = 0;
if ( type(Exp = getopt(exp)) == 2 ) ; else Exp = 0;
Shift = getopt(shift);
if (nk_restriction.subset([Shift],P)) ; else Shift = 0;
if ( type(Order = getopt(order)) == 1 ) ; else Order = -1;
if ( type(Check = getopt(check)) == -1 ) Check = 0;
if ( type(Excp = getopt(excp)) == -1 ) Excp = 0;
if ( Excp && ( IH || Shift || (Order != -1)) ) {
print("error : we can not use excp with inhomo and shift and order");
return ;
}
N = length(F);
Vnum = length(V);
if( length(P) != N ){
print("error : number of parameters is not equal number of factors");
return;
}
DV = map(mkd,V);
FF = 1;
for( I = 0; I < N; I++ )
FF *= F[I];
L = newvect(Vnum + 1);
for (I = 0; I < Vnum; I++)
L[I] = (-1)^(I+1) * diff(FF,V[I]);
L[Vnum] = FF;
L = vtol(L);
L = newsyz.module_syz(L,V,1,0);
L = L[0];
NL = length(L);
DL = [];
for( I = 0; I < NL; I++ ){
M = 0;
for (K = 0; K < Vnum; K++)
M += (-1)^(K+1) * L[I][K] * DV[K];
M -= L[I][Vnum];
for (J = 0; J < N; J++) {
MJ = 0;
for (K = 0; K < Vnum; K++)
MJ += (-1)^K * L[I][K] * diff(F[J],V[K]);
MJ = P[J] * sdiv(MJ,F[J]);
M += MJ;
}
if ( Exp ) {
ME = 0;
for (K = 0; K < Vnum; K++)
ME += (-1)^K * L[I][K] * diff(Exp,V[K]);
M += ME;
}
DL = cons(M,DL);
}
if (Check)
holonomic(DL,V,DV);
FDL = map(nk_restriction.fourier_trans,DL,V,DV);
PV = newvect(N);
PL = [];
for (I = 0; I < N; I++)
if (type(P[I]) == 2) {
PV[I] = var(P[I]);
PL = append(PL,[PV[I]]);
}
W = newvect(Vnum);
for (I = 0; I < Vnum; I++)
W[I] = 1;
W = vtol(W);
WW = newvect(2 *Vnum);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
if ( !Excp ) {
if (Order == -1) {
Result = nk_restriction.generic_bfct_and_gr(FDL, V, DV, W | param = PL);
BF = Result[0];
GB = Result[1];
} else {
VDV = append(V,DV);
VDVH = append(VDV,[h]);
FDLH = map(dp_dtop,map(dp_homo,map(dp_ptod,FDL,VDV)),VDVH);
GBH = nd_weyl_gr(FDLH,VDVH,0,11);
GB = map(subst,GBH,h,1);
}
NN = length(GB);
VarV = append(V, DV);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VarV, WW);
if (Order == -1) {
L = fctr(BF);
BFF = L;
L = cdr(L);
S0 = [];
for (I = 0; I < length(L); I++) {
T = L[I][0];
if ( vars(T) == [s] ) {
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
}
}
if (S0 == [] || S0 < 0)
return 0;
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = nk_restriction.gen(W,S0);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
Index = (gauss(M))[1];
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
Order = length(Base);
}
dp_ord(OldOrder);
print("Order : " + rtostr(Order));
if ( Shift ) {
for (S = 0; S < N; S++)
if (var(PV[S]) == Shift)
break;
Max = dp_td(dp_ptod(F[S], V)) * Order;
} else {
Max = 0;
for (I = 0; I < N; I++)
if (Max < dp_td(dp_ptod(F[I], V)))
Max = dp_td(dp_ptod(F[I], V));
Max = Max * Order;
}
Generator = [];
GB = map(nk_restriction.inv_fourier_trans, GB, V, DV);
GB = map(adjoint, GB, V, DV);
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, Max - OrderGB[I]);
while (B != []) {
Index = car(B);
X = 1;
for (J = 0; J < Vnum; J++) {
X *= V[J]^Index[J];
}
P = nk_restriction.weyl_mul(GB[I], X, VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, DV[J], 0);
if (P != 0) {
P = dp_dtol(dp_ptod(P,V),V);
if ( IH ) {
XX = dp_dtol(dp_ptod(X,V),V);
Generator = cons([P,I,XX], Generator);
} else
Generator = cons(P,Generator);
}
B = cdr(B);
}
}
if ( IH ) {
GenVect = [];
while (Generator != []){
Vect = newvect(NN + 1);
Init = car(Generator);
Vect[0] = Init[0];
Vect[Init[1] + 1] = Init[2];
GenVect = cons(vtol(Vect),GenVect);
Generator = cdr(Generator);
}
}
Elim = map(ltov,nk_restriction.gen(W, Max));
Elim = map(dp_dtol,map(dp_vtoe,Elim),V);
EV = [];
if ( Shift ) {
for (J = 1; J <= Order; J++)
EV = append(EV,[mke_(PV[S],J)]);
} else {
for (I = 0; I < N; I++)
if (PV[I] != 0)
for (J = 1; J <= Order; J++)
EV = append(EV,[mke_(PV[I],J)]);
}
EV = append(EV,[mke_(0,0)]);
Vars = append(Elim, EV);
if ( IH ){
Vect = newvect(NN + 1);
Vect[0] = mke_(0,0) - dp_dtol(dp_ptod(1,V),V);
R = [vtol(Vect)];
if( Shift ) {
for (J = 1; J <= Order; J++) {
Vect = newvect(NN + 1);
Vect[0] = mke_(PV[S],J) - dp_dtol(dp_ptod(F[S]^J,V),V);
R = cons(vtol(Vect),R);
}
} else {
for (I = 0; I < N; I++)
if (PV[I] != 0) {
for (J = 1; J <= Order; J++) {
Vect = newvect(NN + 1);
Vect[0] = mke_(PV[I],J) - dp_dtol(dp_ptod(F[I]^J,V),V);
R = cons(vtol(Vect),R);
}
}
}
G = append(R,GenVect);
G = nd_gr(G,Vars,0,[1,2]);
} else {
R = [mke_(0,0) - dp_dtol(dp_ptod(1,V),V)];
if( Shift ) {
for (J = 1; J <= Order; J++)
R = cons(mke_(PV[S],J) - dp_dtol(dp_ptod(F[S]^J,V),V),R);
} else {
for (I = 0; I < N; I++)
if (PV[I] != 0) {
for (J = 1; J <= Order; J++)
R = cons(mke_(PV[I],J) - dp_dtol(dp_ptod(F[I]^J,V),V),R);
}
}
G = append(R,Generator);
G = nd_gr(G,Vars,0,2);
}
EVP = append(EV,PL);
Eq = [];
if( Shift )
EP = [mke(Shift)];
else
EP = map(mke,PL);
while (G != []) {
Init = car(G);
if ( IH ) {
List = [];
if (nk_restriction.subset(vars(Init[0]), EVP) && Init[0] != 0) {
PP = dp_etop(dp_ptod(Init[0],EV),EP,Order);
List = cons(PP,List);
Crt = 0;
for (I = 0; I < NN; I ++)
if(Init[I + 1] != 0) {
Crt = 1;
List = append(List, [[GB[I],ltop(Init[I + 1],V,Max)]]);
}
if (Crt == 0)
List = append(List, [[]]);
Eq = cons(List,Eq);
}
} else {
if (nk_restriction.subset(vars(Init), EVP)){
Init = dp_etop(dp_ptod(Init,EV),EP,Order);
Eq = cons(Init,Eq);
}
}
G = cdr(G);
}
return Eq;
} else {
NotInt = [];
while (P != []) {
A = car(P);
P = cdr(P);
if (vars(A) != [])
NotInt = cons(A,NotInt);
}
NotInt = reverse(NotInt);
HFDL = map(v_homo,FDL,V,DV,W);
VDV = append(V,DV);
VH = append(V,[h]);
DVH = append(DV,[dh]);
VDVH = append(VH, DVH);
BB = buchberger_weyl(HFDL, VDVH, v_order(Vnum));
GB = map(nm,map(subst,BB[0],h,1));
ExcpList = BB[1];
Vnum2 = 2*Vnum;
dp_weyl_set_weight(W);
M = newmat(Vnum2,Vnum2);
for ( J = 0; J < Vnum2; J++ )
M[0][J] = 1;
for ( I = 1; I < Vnum2; I++ )
M[I][Vnum2-I] = -1;
MW = newmat(Vnum2+1,Vnum2);
for ( J = 0; J < Vnum; J++ )
MW[0][J] = -W[J];
for ( ; J < Vnum2; J++ )
MW[0][J] = W[J-Vnum];
for ( I = 1; I <= Vnum2; I++ )
for ( J = 0; J < Vnum2; J++ )
MW[I][J] = M[I-1][J];
GIN = map(initial_part,GB,VDV,MW,W);
for ( I = 0, T = 0; I < Vnum; I++ )
T += W[I]*V[I]*DV[I];
BF = nk_restriction.weyl_minipoly_by_elim2(GIN,V,DV,PL,T);
L = fctr(BF);
BFF = L;
L = cdr(L);
NL = length(L);
S0 = [];
NotPosInt = [];
for (I = 0; I < NL; I++) {
T = L[I][0];
if ( vars(T) == [s] ) {
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
} else {
NotPosInt = cons(-coef(T, 0, s) / coef(T, 1, s), NotPosInt);
}
}
WW = newvect(Vnum2);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
NN = length(GB);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VDV, WW);
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VDV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Gen = nk_restriction.gen(W,S0);
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
M = excp_gauss(M);
Index = M[1];
ExcpList = append(M[2],ExcpList);
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
Order = length(Base);
dp_ord(OldOrder);
print("Order : " + rtostr(Order));
Max = 0;
for (I = 0; I < N; I++)
if (Max < dp_td(dp_ptod(F[I], V)))
Max = dp_td(dp_ptod(F[I], V));
Max = Max * Order;
VarV = append(V,DV);
Generator = [];
GB = map(nk_restriction.inv_fourier_trans, GB, V, DV);
GB = map(adjoint, GB, V, DV);
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, Max - OrderGB[I]);
while (B != []) {
Index = car(B);
X = 1;
for (J = 0; J < Vnum; J++) {
X *= V[J]^Index[J];
}
P = nk_restriction.weyl_mul(GB[I], X, VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, DV[J], 0);
if (P != 0) {
P = dp_dtol(dp_ptod(P,V),V);
Generator = cons(P,Generator);
}
B = cdr(B);
}
}
Elim = map(ltov,nk_restriction.gen(W, Max));
Elim = map(dp_dtol,map(dp_vtoe,Elim),V);
EV = [];
for (I = 0; I < N; I++)
if (PV[I] != 0)
for (J = 1; J <= Order; J++)
EV = append(EV,[mke_(PV[I],J)]);
EV = append(EV,[mke_(0,0)]);
Vars = append(Elim, EV);
R = [mke_(0,0) - dp_dtol(dp_ptod(1,V),V)];
for (I = 0; I < N; I++)
if (PV[I] != 0) {
for (J = 1; J <= Order; J++)
R = cons(mke_(PV[I],J) - dp_dtol(dp_ptod(F[I]^J,V),V),R);
}
G = append(R,Generator);
G = buchberger(G,Vars,2);
ExcpList = append(ExcpList,G[1]);
G = map(nm,G[0]);
EVP = append(EV,PL);
Eq = [];
EP = map(mke,PL);
while (G != []) {
Init = car(G);
if (nk_restriction.subset(vars(Init), EVP)){
Init = dp_etop(dp_ptod(Init,EV),EP,Order);
Eq = cons(Init,Eq);
}
G = cdr(G);
}
ExcpSet = [];
while (ExcpList != []) {
M = car(ExcpList);
ExcpList = cdr(ExcpList);
if (type(NM = nm(M)) == 2 && !member(NM,ExcpSet) && !member(-NM,ExcpSet))
ExcpSet = cons(nm(M),ExcpSet);
if (type(DN = dn(M)) == 2 && !member(DN,ExcpSet) && !member(-DN,ExcpSet))
ExcpSet = cons(dn(M),ExcpSet);
}
FExcpSet = [];
while (ExcpSet != []) {
ES = car(ExcpSet);
ExcpSet = cdr(ExcpSet);
FExcpSet = append(FExcpSet,fctr(ES));
}
while (FExcpSet != []) {
M = car(FExcpSet)[0];
FExcpSet = cdr(FExcpSet);
if (type(M) == 2 && !member(M,ExcpSet))
ExcpSet = cons(M,ExcpSet);
}
return ["Operator",Eq,"Not integer",NotInt,"Not non-negative integer",NotPosInt,"Not zero",ExcpSet];
}
}
def differential_equation(Exp,F,P,V,X)
{
if ( type(IH = getopt(inhomo)) == -1 ) IH = 0;
Diff = getopt(diff);
if (nk_restriction.subset([Diff],X)) ; else Diff = 0;
if ( type(Order = getopt(order)) == 1 ) ; else Order = -1;
if ( type(Check = getopt(check)) == -1) Check = 0;
if ( type(Excp = getopt(excp)) == -1 ) Excp = 0;
if ( Excp && ( IH || Shift || (Order != -1)) ) {
print("error : we can not use option excp with inhomo and shift and order");
return ;
}
Param = P;
N = length(F);
Vnum = length(V);
Xnum = length(X);
if( length(P) != N ){
print("error : number of parameter is not equal number of factors");
return;
}
DV = map(mkd,V);
FF = 1;
for( I = 0; I < N; I++ )
FF *= F[I];
L = newvect(Vnum + 1);
for (I = 0; I < Vnum; I++)
L[I] = (-1)^(I+1) * diff(FF,V[I]);
L[Vnum] = FF;
L = vtol(L);
L = newsyz.module_syz(L,V,1,0);
L = L[0];
NL = length(L);
DL = [];
for( I = 0; I < NL; I++ ){
M = 0;
for (K = 0; K < Vnum; K++)
M += (-1)^(K+1) * L[I][K] * DV[K];
M -= L[I][Vnum];
for (J = 0; J < N; J++) {
MJ = 0;
for (K = 0; K < Vnum; K++)
MJ += (-1)^K * L[I][K] * diff(F[J],V[K]);
MJ = P[J] * sdiv(MJ,F[J]);
M += MJ;
}
ME = 0;
for (K = 0; K < Vnum; K++)
ME += (-1)^K * L[I][K] * diff(Exp,V[K]);
M += ME;
DL = cons(M,DL);
}
if (Check)
holonomic(DL,V,DV);
FDL = map(nk_restriction.fourier_trans,DL,V,DV);
W = newvect(Vnum);
for (I = 0; I < Vnum; I++)
W[I] = 1;
W = vtol(W);
PV = newvect(N);
PL = [];
for (I = 0; I < N; I++)
if (type(P[I]) == 2) {
PV[I] = var(P[I]);
PL = append(PL,[PV[I]]);
}
WW = newvect(2 * Vnum);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
if ( !Excp ) {
if (Order == -1) {
Result = nk_restriction.generic_bfct_and_gr(FDL, V, DV, W | param = PL);
BF = Result[0];
GB = Result[1];
} else {
VDV = append(V,DV);
VDVH = append(VDV,[h]);
FDLH = map(dp_dtop,map(dp_homo,map(dp_ptod,FDL,VDV)),VDVH);
GBH = nd_weyl_gr(FDLH,VDVH,0,11);
GB = map(subst,GBH,h,1);
}
NN = length(GB);
VarV = append(V, DV);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VarV, WW);
if (Order == -1) {
L = fctr(BF);
BFF = L;
L = cdr(L);
S0 = [];
for (I = 0; I < length(L); I++) {
T = L[I][0];
if ( vars(T) == [s] ) {
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
}
}
if (S0 == [] || S0 < 0)
return 0;
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = nk_restriction.gen(W,S0);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
Index = (gauss(M))[1];
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
dp_ord(OldOrder);
Order = length(Base);
}
print("Order : " + rtostr(Order));
if ( Diff ) {
for (S = 0; S < Xnum; S++)
if (var(X[S]) == Diff)
break;
Max = dp_td(dp_ptod(diff(Exp,X[S]), V)) * Order;
} else {
Max = 0;
for (I = 0; I < Xnum; I++)
if (Max < dp_td(dp_ptod(diff(Exp,X[I]), V)))
Max = dp_td(dp_ptod(diff(Exp,X[I]), V));
Max = Max * Order;
}
Generator = [];
GB = map(nk_restriction.inv_fourier_trans, GB, V, DV);
GB = map(adjoint, GB, V, DV);
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, Max - OrderGB[I]);
while (B != []) {
Index = car(B);
T = 1;
for (J = 0; J < Vnum; J++) {
T *= V[J]^Index[J];
}
P = nk_restriction.weyl_mul(GB[I], T, VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, DV[J], 0);
if (P != 0) {
P = dp_dtol(dp_ptod(P,V),V);
if ( IH ) {
TT = dp_dtol(dp_ptod(T,V),V);
Generator = cons([P,I,TT], Generator);
} else
Generator = cons(P,Generator);
}
B = cdr(B);
}
}
if ( IH ) {
GenVect = [];
while (Generator != []){
Vect = newvect(NN + 1);
Init = car(Generator);
Vect[0] = Init[0];
Vect[Init[1] + 1] = Init[2];
GenVect = cons(vtol(Vect),GenVect);
Generator = cdr(Generator);
}
}
Elim = map(ltov,nk_restriction.gen(W,Max));
Elim = map(dp_dtol,map(dp_vtoe,Elim),V);
DX = [];
if ( Diff ) {
for (J = 1; J <= Order; J++)
DX = append(DX,[mkd_(X[S],J)]);
} else {
for (I = 0; I < Xnum; I++)
for (J = 1; J <= Order; J++)
DX = append(DX,[mkd_(X[I],J)]);
}
DX = append(DX,[mkd_(0,0)]);
Vars = append(Elim, DX);
if ( IH ){
Vect = newvect(NN + 1);
Vect[0] = mkd_(0,0) - dp_dtol(dp_ptod(1,V),V);
R = [vtol(Vect)];
if( Diff ) {
H = diff(Exp,X[S]);
for (J = 1; J <= Order; J++) {
Vect = newvect(NN + 1);
Vect[0] = mkd_(X[S],J) - dp_dtol(dp_ptod(H,V),V);
R = cons(vtol(Vect),R);
H = diff(Exp,X[S]) * H + diff(H,X[S]);
}
} else {
for (I = 0; I < Xnum; I++) {
H = diff(Exp,X[I]);
for (J = 1; J <= Order; J++) {
Vect = newvect(NN + 1);
Vect[0] = mkd_(X[I],J) - dp_dtol(dp_ptod(H,V),V);
R = cons(vtol(Vect),R);
H = diff(Exp,X[I]) * H + diff(H,X[I]);
}
}
}
G = append(R,GenVect);
G = nd_gr(G,Vars,0,[1,2]);
} else {
R = [mkd_(0,0) - dp_dtol(dp_ptod(1,V),V)];
if( Diff ) {
H = diff(Exp,X[S]);
for (J = 1; J <= Order; J++) {
R = cons(mkd_(X[S],J) - dp_dtol(dp_ptod(H,V),V),R);
H = diff(Exp,X[S]) * H + diff(H,X[S]);
}
} else {
for (I = 0; I < Xnum; I++) {
H = diff(Exp,X[I]);
for (J = 1; J <= Order; J++) {
R = cons(mkd_(X[I],J) - dp_dtol(dp_ptod(H,V),V),R);
H = diff(Exp,X[I]) * H + diff(H,X[I]);
}
}
}
G = append(R,Generator);
G = nd_gr(G,Vars,0,2);
}
DXXP = append(append(DX,X),vars(Param));
Eq = [];
if( Diff )
DDX = [mkd(Diff)];
else
DDX = map(mkd,X);
while (G != []) {
Init = car(G);
if ( IH ) {
List = [];
if (nk_restriction.subset(vars(Init[0]), DXXP) && Init[0] != 0) {
PP = dp_etop(dp_ptod(Init[0],DX),DDX,Order);
List = cons(PP,List);
Crt = 0;
for (I = 0; I < NN; I ++)
if(Init[I + 1] != 0) {
Crt = 1;
List = append(List, [[GB[I],ltop(Init[I + 1],V,Max)]]);
}
if (Crt == 0)
List = append(List, [[]]);
Eq = cons(List,Eq);
}
} else {
if (nk_restriction.subset(vars(Init), DXXP)){
Init = dp_etop(dp_ptod(Init,DX),DDX,Order);
Eq = cons(Init,Eq);
}
}
G = cdr(G);
}
return Eq;
} else {
NotInt = [];
while (P != []) {
A = car(P);
P = cdr(P);
if (vars(A) != [])
NotInt = cons(A,NotInt);
}
NotInt = reverse(NotInt);
HFDL = map(v_homo,FDL,V,DV,W);
VDV = append(V,DV);
VH = append(V,[h]);
DVH = append(DV,[dh]);
VDVH = append(VH, DVH);
BB = buchberger_weyl(HFDL, VDVH, v_order(Vnum));
GB = map(nm,map(subst,BB[0],h,1));
ExcpList = BB[1];
Vnum2 = 2*Vnum;
dp_weyl_set_weight(W);
M = newmat(Vnum2,Vnum2);
for ( J = 0; J < Vnum2; J++ )
M[0][J] = 1;
for ( I = 1; I < Vnum2; I++ )
M[I][Vnum2-I] = -1;
MW = newmat(Vnum2+1,Vnum2);
for ( J = 0; J < Vnum; J++ )
MW[0][J] = -W[J];
for ( ; J < Vnum2; J++ )
MW[0][J] = W[J-Vnum];
for ( I = 1; I <= Vnum2; I++ )
for ( J = 0; J < Vnum2; J++ )
MW[I][J] = M[I-1][J];
GIN = map(initial_part,GB,VDV,MW,W);
for ( I = 0, T = 0; I < Vnum; I++ )
T += W[I]*V[I]*DV[I];
BF = nk_restriction.weyl_minipoly_by_elim2(GIN,V,DV,PL,T);
L = fctr(BF);
BFF = L;
L = cdr(L);
NL = length(L);
S0 = [];
NotPosInt = [];
for (I = 0; I < NL; I++) {
T = L[I][0];
if ( vars(T) == [s] ) {
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
} else {
NotPosInt = cons(-coef(T, 0, s) / coef(T, 1, s), NotPosInt);
}
}
WW = newvect(Vnum2);
for (I = 0; I < Vnum; I++)
WW[I] = -1;
for (I = 0; I < Vnum; I++)
WW[Vnum + I] = 1;
WW = vtol(WW);
NN = length(GB);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = nk_restriction.ord_w(GB[I], VDV, WW);
Generator = [];
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < Vnum; J++) {
D *= DV[J]^Index[J];
}
P = nk_restriction.weyl_mul(D, GB[I], VDV);
for (J = 0; J < Vnum; J++)
P = subst(P, V[J], 0);
if (P != 0)
Generator = cons(P, Generator);
B = cdr(B);
}
}
Gen = nk_restriction.gen(W,S0);
Mat = newmat(Vnum,Vnum);
for (I = 0; I < Vnum; I++)
Mat[I][I] = 1;
OldOrder = dp_ord();
dp_ord(Mat);
Gen = map(dp_vtoe,map(ltov,Gen));
Gennum = length(Gen);
M = [];
while( Generator != [] ){
DP = dp_ptod(car(Generator),DV);
LP = newvect(Gennum);
for( I = 0; I < Gennum; I++ ){
if( DP == 0 )
break;
if( dp_ht(DP) == Gen[I] ){
LP[I] = dp_hc(DP);
DP = dp_rest(DP);
}
}
M = cons(vtol(LP),M);
Generator = cdr(Generator);
}
M = newmat(length(M),Gennum,M);
M = excp_gauss(M);
Index = M[1];
ExcpList = append(M[2],ExcpList);
Crt = newvect(Gennum);
for( I = 0; I < length(Index); I++ )
if(Index[I] > -1)
Crt[Index[I]] = 1;
Base = [];
for( I = 0; I < Gennum; I++ ){
if( Crt[I] == 0 )
Base = append(Base,[Gen[I]]);
}
Order = length(Base);
dp_ord(OldOrder);
print("Order : " + rtostr(Order));
Max = 0;
for (I = 0; I < N; I++)
if (Max < dp_td(dp_ptod(F[I], V)))
Max = dp_td(dp_ptod(F[I], V));
Max = Max * Order;
VarV = append(V,DV);
Generator = [];
GB = map(nk_restriction.inv_fourier_trans, GB, V, DV);
GB = map(adjoint, GB, V, DV);
for (I = 0; I < NN; I++ ) {
B = nk_restriction.gen(W, Max - OrderGB[I]);
while (B != []) {
Index = car(B);
T = 1;
for (J = 0; J < Vnum; J++) {
T *= V[J]^Index[J];
}
P = nk_restriction.weyl_mul(GB[I], T, VarV);
for (J = 0; J < Vnum; J++)
P = subst(P, DV[J], 0);
if (P != 0) {
P = dp_dtol(dp_ptod(P,V),V);
Generator = cons(P,Generator);
}
B = cdr(B);
}
}
Elim = map(ltov,nk_restriction.gen(W, Max));
Elim = map(dp_dtol,map(dp_vtoe,Elim),V);
DX = [];
for (I = 0; I < Xnum; I++)
for (J = 1; J <= Order; J++)
DX = append(DX,[mkd_(X[I],J)]);
DX = append(DX,[mkd_(0,0)]);
Vars = append(Elim, DX);
R = [mkd_(0,0) - dp_dtol(dp_ptod(1,V),V)];
for (I = 0; I < Xnum; I++) {
H = diff(Exp,X[I]);
for (J = 1; J <= Order; J++) {
R = cons(mkd_(X[I],J) - dp_dtol(dp_ptod(H,V),V),R);
H = diff(Exp,X[I]) * H + diff(H,X[I]);
}
}
G = append(R,Generator);
G = buchberger(G,Vars,2);
ExcpList = append(ExcpList,G[1]);
G = map(nm,G[0]);
DXXP = append(append(DX,X),vars(Param));
Eq = [];
DDX = map(mkd,X);
while (G != []) {
Init = car(G);
if (nk_restriction.subset(vars(Init), DXXP)){
Init = dp_etop(dp_ptod(Init,DX),DDX,Order);
Eq = cons(Init,Eq);
}
G = cdr(G);
}
ExcpSet = [];
while (ExcpList != []) {
M = car(ExcpList);
ExcpList = cdr(ExcpList);
if (type(NM = nm(M)) == 2 && !member(NM,ExcpSet) && !member(-NM,ExcpSet))
ExcpSet = cons(nm(M),ExcpSet);
if (type(DN = dn(M)) == 2 && !member(DN,ExcpSet) && !member(-DN,ExcpSet))
ExcpSet = cons(dn(M),ExcpSet);
}
FExcpSet = [];
while (ExcpSet != []) {
ES = car(ExcpSet);
ExcpSet = cdr(ExcpSet);
FExcpSet = append(FExcpSet,fctr(ES));
}
while (FExcpSet != []) {
M = car(FExcpSet)[0];
FExcpSet = cdr(FExcpSet);
if (type(M) == 2 && !member(M,ExcpSet))
ExcpSet = cons(M,ExcpSet);
}
return ["Operator",Eq,"Not integer",NotInt,"Not non-negative integer",NotPosInt,"Not zero",ExcpSet];
}
}
def monomial_sum(I,M)
{
J = [M];
while(I != []){
if(srem(car(I),M) != 0)
J = cons(I[0],J);
I = cdr(I);
}
J = reverse(J);
return J;
}
def monomial_quotient(I,M)
{
L = length(I);
I = newvect(L,I);
for( N = 0; N < L; N++ ){
if(deg(I[N],M) > 0)
I[N] = sdiv(I[N],M);
if(I[N] == 1)
return [1];
}
I = vtol(I);
return I;
}
def combination(F,N){
G = 1;
if(N == -1){
if(F == -1)
return 1;
else
return 0;
}else if(N == 0)
return 1;
else{
while(N > 0){
G *= F/N;
F--;
N--;
}
}
return G;
}
def homogeneous(I,Order){
I=nd_gr(I,Order,0,0);
L=length(I);
I=newvect(L,I);
for(I0=0;I0<L;I0++){
I[I0]=dp_homo(dp_ptod(I[I0],Order));
}
Order=append(Order,[h]);
for(I0=0;I0<L;I0++){
I[I0]=dp_dtop(I[I0],Order);
}
return(vtol(I));
}
def monomial_hilbert_poly(I,Order)
{
I = nd_gr(I,Order,0,0);
S = length(I);
T = length(Order);
if(I == [0])
return 1;
if(I == [1])
return 0;
for( M = 0; M < S; M++ ){
if(dp_td(dp_ptod(I[M],Order)) > 1)
break;
}
if( M == S ){
F = (1-t)^S;
return F;
}
Crt = 0;
for( N = 0; N < T; N++ ){
for( M = 0; M < S; M++ ){
if((dp_td(dp_ptod(I[M],Order)) > 1) && (deg(I[M],Order[N]) >= 1)){
Crt = 1;
break;
}
}
if(Crt == 1)
break;
}
I1 = monomial_sum(I,Order[N]);
I2 = monomial_quotient(I,Order[N]);
F1 = monomial_hilbert_poly(I1,Order);
F2 = monomial_hilbert_poly(I2,Order);
F = F1+t*F2;
return F;
}
def hilbelt_function(I,Order){
List = [];
Plist = [];
I = nd_gr(I,Order,0,0);
L = length(I);
I = newvect(L,I);
for( I0 = 0; I0 < L; I0++ ){
I[I0] = dp_hm(dp_ptod(I[I0],Order));
I[I0] = dp_dtop(I[I0],Order);
}
I = vtol(I);
return monomial_hilbert_poly(I,Order);
}
def hilbelt_poly(I,Order){
F = hilbelt_function(I,Order);
N = length(Order);
G = F;
while((F = tdiv(G,1-t)) != 0){
G = F;
N--;
}
G = dp_ptod(G,[t]);
P = 0;
D = dp_td(G);
for( I0 = D; I0 >= 0; I0-- ){
P += dp_hc(G) * combination(N-1+x-I0,N-1);
G = dp_rest(G);
}
return P;
}
def affine_hilbelt_poly(I,Order){
I = homogeneous(I,Order);
Order = append(Order,[h]);
return hilbelt_poly(I,Order);
}
def holonomic(I,V,DV){
L = length(V);
VDV = append(V,DV);
M = newmat(2*L+1,2*L);
for( J0 = L; J0 < 2*L ;J0++ )
M[0][J0] = 1;
for( J0 = 0; J0 < 2*L; J0++ )
M[1][J0] = 1;
for( I0 = 2; I0 < 2*L+1; I0++)
M[I0][2*L-I0+1] = -1;
G = nd_weyl_gr(I,VDV,0,M);
N = length(G);
OldOrder = dp_ord();
dp_ord(M);
G = map(dp_ptod,G,VDV);
DIN = map(dp_ht,G);
IN = map(dp_dtop,DIN,VDV);
P = affine_hilbelt_poly(IN,VDV);
print("Hilbert polynomial : ",0);
print(P);
D = deg(P,x);
dp_ord(OldOrder);
if(D == L){
print("holonomic : Yes");
IN = ltov(IN);
for( I0 = 0; I0 < N; I0++ ){
for( K = 0; K < L; K++ )
IN[I0] = subst(IN[I0],V[K],1);
}
IN = vtol(IN);
if(dp_td(dp_ptod(IN[0],DV)) == 0){
R = 0;
}else{
SM = dp_mbase(map(dp_ptod,IN,DV));
R = length(SM);
}
print("holonomic rank : ",0);
print(R);
return map(dp_dtop,SM,DV);
}else{
print("holonomic : No");
return -1;
}
}
def mkt(L){
return strtov("t" + rtostr(L));
}
def mkx(I,J){
return strtov("x" + rtostr(I) + rtostr(J));
}
def mky(L){
return strtov("y" + rtostr(L));
}
def example6_5()
{
return difference_equation([x,y,1-x-y],[a,b,c],[x,y]|inhomo=1);
}
def example6_6()
{
print("System:");
print(difference_equation([x,y,1-x-y,1-2*x],[a,b,c,d],[x,y]));
print("Equation with respect to E_a:");
print(difference_equation([x,y,1-x-y,1-2*x],[a,b,c,d],[x,y]|shift = a));
print("Equation with respect to E_b:");
print(difference_equation([x,y,1-x-y,1-2*x],[a,b,c,d],[x,y]|shift = b));
print("Equation with respect to E_c:");
print(difference_equation([x,y,1-x-y,1-2*x],[a,b,c,d],[x,y]|shift = c));
print("Equation with respect to E_d:");
print(difference_equation([x,y,1-x-y,1-2*x],[a,b,c,d],[x,y]|shift = d));
}
def example6_9()
{
return difference_equation([x+y*z,y+z*x,z+x*y],[a,b,c],[x,y]);
}
/* (1,,1,1,-1,-1,-1)-filteration に関する graded module M_S の次元を計算 */
def example6_10(S)
{
return twisted_log_cohomology([x*z+y,x^4+y^5+x*y^4],[a,b],[x,y]|S0 = S);
}
def example7_1()
{
return differential_equation(t1*x+t2*y,[x^2+y^2-1],[-1],[x,y],[t1,t2]);
}
def example7_2()
{
return differential_equation(t1*x^2+t2*y^2,[x+y],[a],[x,y],[t1,t2]);
}
def example7_4()
{
print("System");
print(differential_equation(x1*t^2+x2*t^3,[t],[a],[t],[x1,x2]));
print("Equation with respect to x1");
print(differential_equation(x1*t^2+x2*t^3,[t],[a],[t],[x1,x2]|diff = x1));
print("Equation with respect to x2");
print(differential_equation(x1*t^2+x2*t^3,[t],[a],[t],[x1,x2]|diff = x2));
}
def example7_5()
{
return differential_equation(x1*t2*t3+x2*t1^2*t3+x3*t3,[t1,t2,t3],[b1,b2,b3],[t1,t2,t3],[x1,x2,x3]);
}
/* S^N 上の Fisher-Bingham integral に付随する twisted log cohomology を計算*/
def example7_6(N)
{
H = -1;
for (I = 1; I <= N+1; I++)
H += mkt(I)^2;
Exp = 0;
for (I = 1; I <= N+1; I++)
Exp += mkx(I,I)*mkt(I)^2 + mky(I)*mkt(I);
for (J = 1; J <= N+1; J++)
for (I = 1; I < J; I++)
Exp += 2*mkx(I,J)*mkt(I)*mkt(J);
T = [];
Y = [];
X = [];
for (I = 1; I <= N+1; I++){
T = cons(mkt(I),T);
Y = cons(mky(I),Y);
}
T = reverse(T);
Y = reverse(Y);
for (J = 1; J <= N+1; J++)
for (I = 1; I < J; I++)
X = cons(mkx(I,J),X);
X = reverse(X);
for (I = 0; I <= N ;I++)
X = cons(mkx(N+1-I,N+1-I),X);
for (I = 0; I <= N; I++)
Exp = subst(Exp,X[I],1/(I+1));
for (I = N+1; I < length(X); I++)
Exp = subst(Exp,X[I],1/2);
for (I = 0; I <= N; I++)
Exp = subst(Exp,Y[I],1/2);
T0 = time();
Base = twisted_log_cohomology([H],[0],T | exp = Exp);
T1 = time();
return ["Dimension",length(Base),"Time",T1[0]-T0[0]];
}
/*以下は timing data 用の関数*/
def time_f1()
{
F = [x,y,1-x-y,1-2*x];
P = [a,b,c,d];
V = [x,y];
T0 = time();
difference_equation(F,P,V);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f2()
{
F = [x,y,1-x-y,1-2*x-3*y];
P = [a,b,c,d];
V = [x,y];
T0 = time();
difference_equation(F,P,V);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f3()
{
F = [x,y,1-x-y,1-2*x-3*y,1-x];
P = [a,b,c,d,e];
V = [x,y];
T0 = time();
difference_equation(F,P,V);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f4()
{
F = [x,y,1-x-y,1-2*x-3*y,1-2*y];
P = [a,b,c,d,e];
V = [x,y];
T0 = time();
difference_equation(F,P,V);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_g(P,Q)
{
F = x^P+y^Q+x*y^(Q-1);
T0 = time();
E = ns_twistedlog.difference_equation([F],[a],[x,y]);
T1 = time();
return ["Integrand",[F,a],"Time",T1[0]-T0[0]];
}
def time_h(P,Q)
{
F = x^P+y^Q+x*y^(Q-1);
G = x+y;
T0 = time();
E = ns_twistedlog.difference_equation([F,G],[a,b],[x,y]);
T1 = time();
return ["Integrand",[[F,G],[a,b]],"Time",T1[0]-T0[0]];
}
def time_f1_a()
{
F = [x,y,1-x-y,1-2*x];
P = [a,b,c,d];
V = [x,y];
T0 = time();
difference_equation(F,P,V|shift = P[0]);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f2_a()
{
F = [x,y,1-x-y,1-2*x-3*y];
P = [a,b,c,d];
V = [x,y];
T0 = time();
difference_equation(F,P,V|shift = P[0]);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f3_a()
{
F = [x,y,1-x-y,1-2*x-3*y,1-x];
P = [a,b,c,d,e];
V = [x,y];
T0 = time();
difference_equation(F,P,V|shift = P[0]);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_f4_a()
{
F = [x,y,1-x-y,1-2*x-3*y,1-2*y];
P = [a,b,c,d,e];
V = [x,y];
T0 = time();
difference_equation(F,P,V|shift = P[0]);
T1 = time();
return ["Integrand",[F,P],"Time",T1[0]-T0[0]];
}
def time_h_a(P,Q)
{
F = x^P+y^Q+x*y^(Q-1);
G = x+y;
T0 = time();
E = ns_twistedlog.difference_equation([F,G],[a,b],[x,y]|shift = a);
T1 = time();
return ["Integrand",[[F,G],[a,b]],"Time",T1[0]-T0[0]];
}
endmodule$
end$![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ns_twistedlog.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / src