File: [local] / OpenXM / src / asir-contrib / packages / src / nk_restriction.rr (download)
Revision 1.23, Sat May 14 06:37:43 2016 UTC (8 years ago) by nakayama
Branch: MAIN
CVS Tags: HEAD
Changes since 1.22: +4 -3
lines
nk_module_gr2.rr : compute module integration, restriction.
nk_restriction.rr : import nk_module_gr2.rr
|
/*
#define DEBUG_RESTRICTION
*/
import("bfct")$
import("noro_module_syz.rr")$ /* for newsyz.module_syz() */
import("nk_module_gr2.rr")$ /* module_integration, module_restriction */
module nk_restriction;
extern NK_RESTRICTION_NOLOG;
localf version;
localf generic_bfct_and_gr;
localf initial_ideal;
localf initial_w, initial_part_w;
localf restriction;
localf max_int_root;
localf test_max_int_root;
localf gen;
localf cons_rev;
localf test_gen;
localf ord_w;
localf weyl_mul;
localf test_restriction;
localf test_restriction2;
localf test_restriction3;
localf restriction_ideal;
localf test_restriction_ideal;
localf restriction_ideal_internal;
localf test_restriction_ideal_internal;
localf test_nd_weyl_gr;
localf fourier_trans;
localf test_fourier_trans;
localf inv_fourier_trans;
localf test_inv_fourier_trans;
localf integration;
localf test_integration;
localf test_integration2;
localf test_integration3;
localf test_integration4;
localf test_integration5;
localf integration_ideal;
localf integration_ideal_internal;
localf integration_ideal_;
localf integration_ideal_internal_;
localf test_integration_ideal;
localf test_integration_ideal2;
localf test_integration_ideal5;
localf test_integration_ideal5_;
localf test_integration_ideal6;
localf test_integration_ideal7;
localf test_integration_ideal8;
localf coord_conv;
localf check_coord_conv, check_coord_conv2;
localf dummy_var;
localf ann_mul;
localf test_ann_mul, test_ann_mul2, test_ann_mul3;
localf subst_vect;
localf fano2;
localf fano3;
localf p_to_mat;
localf test_p_to_mat;
localf fano_polytope_a;
localf test_fano_polytope_a;
localf fano_polytope;
localf a_to_poly;
localf test_a_to_poly;
localf var_list;
localf test_var_list;
localf d_var;
localf int_fano;
/* inhomogeneous */
localf inhomo_part;
localf trans_inhomo, trans_inhomo_sub, pick_gen_relation, is_nonzero_const;
localf separate_var, make_elim_order_mat, make_weight_mat;
localf action_exp, action_exp_1;
localf subset, zero_inhomo;
localf integration_exp;
localf test_integration_exp;
localf time_print, time_round;
localf ost_integration_ideal, ost_integration_ideal2;
localf test_ost_integration_ideal;
localf ost_integration_exp, ost_integration_exp2;
localf ann_heaviside;
/* difference */
localf mellin_trans, mellin_trans2, inv_mellin_trans, inv_mellin_trans2;
localf mellin1, dsubst, theta, normal;
localf integration_ideal_shift;
localf inhomo_inv_mellin, inhomo_inv_mellin_sub;
localf car_func;
localf test_integration_ideal_shift, test_integration_ideal_shift2;
localf test_integration_ideal_shift3, test_integration_ideal_shift4;
localf ost_sum, test_ost_sum, test_ost_sum2;
/* parameter */
localf weyl_minipoly_by_elim, weyl_minipoly_by_elim2, weyl_minipoly_by_elim3;
def version()
{
T = "$Revision: 1.23 $";
S = str_chr(T,0," ")+1;
return sub_str(T,S,str_chr(T,S," ")-1);
}
/* source from lib/bfct generic_bfct */
/* generic b-fucntion に加えて F の <_(-W,W) についての GB を返す */
/* F : ideal, V : var list(x), DV : var list(dx), W : weight vect */
def generic_bfct_and_gr(F,V,DV,W)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( member(s,append(vars(F),vars(Param))) )
error("generic_bfct_and_gr : the variable 's' is reserved.");
N = length(V);
N2 = N*2;
/* If W is a list, convert it to a vector */
if ( type(W) == 4 )
W = newvect(length(W),W);
/* dp_weyl_set_weight とは...? */
dp_weyl_set_weight(W);
/* create a term order M in D<x,d> (DRL) */
M = newmat(N2,N2);
for ( J = 0; J < N2; J++ )
M[0][J] = 1;
for ( I = 1; I < N2; I++ )
M[I][N2-I] = -1;
VDV = append(V,DV);
/* create a non-term order MW in D<x,d> */
MW = newmat(N2+1,N2);
for ( J = 0; J < N; J++ )
MW[0][J] = -W[J];
for ( ; J < N2; J++ )
MW[0][J] = W[J-N];
for ( I = 1; I <= N2; I++ )
for ( J = 0; J < N2; J++ )
MW[I][J] = M[I-1][J];
/* create a homogenized term order MWH in D<x,d,h> */
MWH = newmat(N2+2,N2+1);
for ( J = 0; J <= N2; J++ )
MWH[0][J] = 1;
for ( I = 1; I <= N2+1; I++ )
for ( J = 0; J < N2; J++ )
MWH[I][J] = MW[I-1][J];
/* homogenize F */
VDVH = append(VDV,[h]);
FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
/* compute a groebner basis of FH w.r.t. MWH */
dp_gr_flags(["Top",1,"NoRA",1]);
if ( !NK_RESTRICTION_NOLOG ) T0 = time(); /*tstart();*/
if ( Param && Param != [] && ::version() < 20100206 ) {
/* これ以前の nd_weyl_gr は係数体が有理関数体では動かない */
if ( !NK_RESTRICTION_NOLOG ) print("-- dp_weyl_gr_main :", 2);
GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
} else {
if ( !NK_RESTRICTION_NOLOG ) print("-- nd_weyl_gr :", 2);
GH = nd_weyl_gr(FH,VDVH,0,11);
}
if ( !NK_RESTRICTION_NOLOG ) time_print(T0); /*tstop();*/
dp_gr_flags(["Top",0,"NoRA",0]);
/* dehomigenize GH */
G = map(subst,GH,h,1);
/* G is a groebner basis w.r.t. a non term order MW */
/* take the initial part w.r.t. (-W,W) */
GIN = map(initial_part,G,VDV,MW,W);
/* GIN is a groebner basis w.r.t. a term order M */
/* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
/* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
for ( I = 0, T = 0; I < N; I++ )
T += W[I]*V[I]*DV[I];
if ( S0 >= 0 ) {
if ( !NK_RESTRICTION_NOLOG ) print("skip the computation of generic b-function");
return [s-S0, G];
}
if ( !NK_RESTRICTION_NOLOG ) T0 = time(); /*tstart();*/
if ( Param ) {
/* パラメータ付きは消去法で計算 */
B = weyl_minipoly_by_elim2(GIN,V,DV,Param,T);
if ( !NK_RESTRICTION_NOLOG ) print("-- weyl_minipoly_by_elim :", 2);
} else {
B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
if ( !NK_RESTRICTION_NOLOG ) print("-- weyl_minipoly :", 2);
}
if ( !NK_RESTRICTION_NOLOG ) time_print(T0); /*tstop();*/
return [B, G];
}
/* source from lib/bfct generic_bfct */
/* F の <_(-W,W) についての GB とその initial part を返す */
/* F : ideal, V : var list(x), DV : var list(dx), W : weight vect */
def initial_ideal(F,V,DV,W)
{
N = length(V);
N2 = N*2;
/* If W is a list, convert it to a vector */
if ( type(W) == 4 )
W = newvect(length(W),W);
/* dp_weyl_set_weight とは...? */
dp_weyl_set_weight(W);
/* create a term order M in D<x,d> (DRL) */
M = newmat(N2,N2);
for ( J = 0; J < N2; J++ )
M[0][J] = 1;
for ( I = 1; I < N2; I++ )
M[I][N2-I] = -1;
VDV = append(V,DV);
/* create a non-term order MW in D<x,d> */
MW = newmat(N2+1,N2);
for ( J = 0; J < N; J++ )
MW[0][J] = -W[J];
for ( ; J < N2; J++ )
MW[0][J] = W[J-N];
for ( I = 1; I <= N2; I++ )
for ( J = 0; J < N2; J++ )
MW[I][J] = M[I-1][J];
/* create a homogenized term order MWH in D<x,d,h> */
MWH = newmat(N2+2,N2+1);
for ( J = 0; J <= N2; J++ )
MWH[0][J] = 1;
for ( I = 1; I <= N2+1; I++ )
for ( J = 0; J < N2; J++ )
MWH[I][J] = MW[I-1][J];
/* homogenize F */
VDVH = append(VDV,[h]);
FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
/* compute a groebner basis of FH w.r.t. MWH */
dp_gr_flags(["Top",1,"NoRA",1]);
GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
dp_gr_flags(["Top",0,"NoRA",0]);
/* dehomigenize GH */
G = map(subst,GH,h,1);
/* G is a groebner basis w.r.t. a non term order MW */
/* take the initial part w.r.t. (-W,W) */
GIN = map(initial_part,G,VDV,MW,W);
/* GIN is a groebner basis w.r.t. a term order M */
/* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
return [G, GIN];
}
def initial_w(B,V,W)
{
MW = make_weight_mat(W);
GB = nd_weyl_gr(B,V,0,MW);
return map(initial_part_w,GB,V,MW,W);
}
def initial_part_w(F,V,MW,W)
{
N2 = length(V);
N = N2/2;
dp_ord(MW);
DF = dp_ptod(F,V);
R = dp_hm(DF);
DF = dp_rest(DF);
E = dp_etov(R);
for ( I = 0, TW = 0; I < N2; I++ )
TW += W[I]*E[I];
RW = TW;
for ( ; DF; DF = dp_rest(DF) ) {
E = dp_etov(DF);
for ( I = 0, TW = 0; I < N2; I++ )
TW += W[I]*E[I];
if ( TW == RW )
R += dp_hm(DF);
else if ( TW < RW )
break;
else
error("initial_part : cannot happen");
}
return dp_dtop(R,V);
}
/* 消去法による最小多項式計算 (パラメータを含む場合) */
def weyl_minipoly_by_elim(GIN,V,DV,Param,T)
{
TMP_S = s; TMP_DS = ds;
SV = cons(TMP_S,V); DSV = cons(TMP_DS,DV);
DParam = map(d_var,Param);
VV = append(append(SV,Param),append(DSV,DParam));
NP = length(Param);
NV = length(V);
N = NP+NV+1;
M = make_elim_order_mat(0,N);
G = cons(TMP_S-T,GIN);
for ( K = 0; K < 2; K++ ) {
for ( I = 0; I < NV; I++ ) {
IND = 2*NV+NP+1-I-K*(NV+NP+1);
M[0][IND] = 1;
GB = nd_weyl_gr(G,VV,0,M);
G = [];
while ( GB != [] ) {
if ( member(VV[IND],vars(car(GB))) )
break;
G = cons(car(GB),G);
GB = cdr(GB);
}
G = reverse(G);
}
}
if ( !subset(vars(G),cons(TMP_S,Param)) )
error("weyl_minipoly_by_elim : invalid elimination");
if ( G == [] )
error("weyl_minipoly_by_elim : b-function does not exist");
if ( length(G) != 1 ) {
if ( !NK_RESTRICTION_NOLOG ) print(G);
G = nd_gr(G,[TMP_S],0,0);
if ( !NK_RESTRICTION_NOLOG ) print(G);
}
if ( length(G) != 1 )
error("weyl_minipoly_by_elim : invalid b-function");
return car(G);
}
/* 消去法による最小多項式計算 (パラメータを含む場合) */
def weyl_minipoly_by_elim2(GIN,V,DV,Param,T)
{
TMP_S = s; TMP_DS = ds;
SV = cons(TMP_S,V); DSV = cons(TMP_DS,DV);
DParam = map(d_var,Param);
VV = append(SV,DSV);
NP = length(Param);
NV = length(V);
N = NV+1;
M = make_elim_order_mat(0,N);
G = cons(TMP_S-T,GIN);
for ( K = 0; K < 2; K++ ) {
for ( I = 0; I < NV; I++ ) {
IND = 2*NV+1-I-K*(NV+1);
M[0][IND] = 1;
if ( Param != [] && ::version() < 20100206 )
GB = dp_weyl_gr_main(G,VV,0,0,M);
else
GB = nd_weyl_gr(G,VV,0,M);
G = [];
while ( GB != [] ) {
if ( member(VV[IND],vars(car(GB))) )
break;
G = cons(car(GB),G);
GB = cdr(GB);
}
G = reverse(G);
}
}
if ( !subset(vars(G),cons(TMP_S,Param)) )
error("weyl_minipoly_by_elim : invalid elimination");
if ( G == [] )
error("weyl_minipoly_by_elim : b-function does not exist");
if ( length(G) != 1 )
error("weyl_minipoly_by_elim : invalid b-function");
return car(G);
}
/* 消去法による最小多項式計算 (パラメータを含む場合) */
def weyl_minipoly_by_elim3(GIN,V,DV,Param,T)
{
TMP_S = s; TMP_DS = ds;
SV = cons(TMP_S,V); DSV = cons(TMP_DS,DV);
DParam = map(d_var,Param);
VV = append(append(SV,Param),append(DSV,DParam));
NP = length(Param);
NV = length(V);
N = NP+NV+1;
M = make_elim_order_mat(0,N);
G = cons(TMP_S-T,GIN);
for ( K = 0; K < 2; K++ ) {
for ( I = 0; I < NV; I++ ) {
IND = 2*NV+NP+1-I-K*(NV+NP+1);
M[0][IND] = 1;
GB = dp_weyl_gr_main(G,VV,0,0,M);
G = [];
while ( GB != [] ) {
if ( member(VV[IND],vars(car(GB))) )
break;
G = cons(car(GB),G);
GB = cdr(GB);
}
G = reverse(G);
}
}
if ( !subset(vars(G),cons(TMP_S,Param)) )
error("weyl_minipoly_by_elim : invalid elimination");
if ( G == [] )
error("weyl_minipoly_by_elim : b-function does not exist");
if ( length(G) != 1 ) {
if ( !NK_RESTRICTION_NOLOG ) print(G);
G = nd_gr(G,[TMP_S],0,0);
if ( !NK_RESTRICTION_NOLOG ) print(G);
}
if ( length(G) != 1 )
error("weyl_minipoly_by_elim : invalid b-function");
return car(G);
}
/*
* restriction module の計算([SST, Alg5.2.8])
* Id : ideal
* VL : var list (x)
* DVL: var list (dx)
* W : weight vector
* VL と長さは同じ
* W[0], ..., W[M] が正で、それ以降の要素は 0 としておく
*
* Option
* param : パラメータとみなす変数のリスト
* s0 : s0 >= 0 のとき generic b-function の計算を行わず
* s-s0 を generic b-function として計算する
*/
def restriction(Id, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( !NK_RESTRICTION_NOLOG ) tstart();
Result = generic_bfct_and_gr(Id, VL, DVL, W|param=Param, s0=S0);
if ( !NK_RESTRICTION_NOLOG ) print("-- generic_bfct_and_gr :",2);
if ( !NK_RESTRICTION_NOLOG ) tstop();
BF = Result[0];
GB = Result[1];
if ( !NK_RESTRICTION_NOLOG ) tstart();
/* BF の最大整数根 S0 を取り出す */
L = fctr(BF);
BFF = L;
/* L の最初の要素は係数なので捨てる */
L = cdr(L);
N = length(L);
S0 = [];
for (I = 0; I < N; I++) {
T = L[I][0];
if ( vars(T) == [s] ) { /* パラメータは generic とする */
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
}
}
if ( !NK_RESTRICTION_NOLOG ) {
print("generic bfct : " + rtostr(BFF), 1);
print("S0 : " + rtostr(S0), 1);
}
/* b-function が非負整数根を持たない場合([SST, Cor5.2.7]) */
if (S0 == [] || S0 < 0) {
if ( IH )
S0 = 0;
else
return 0;
}
for (M = 0; M < length(W); M++)
if (W[M] <= 0)
break;
/* W の正の要素の部分だけとりだしたもの WW */
WW = newvect(M);
for (I = 0; I < M; I++)
WW[I] = W[I];
WW = vtol(WW);
/* B_{S0} の生成 */
Base = gen(WW, S0);
if ( !NK_RESTRICTION_NOLOG ) print("B_{S0} length : " + rtostr(length(Base)),1);
V = newvect(length(W), W);
MinusW = vtol(-V);
WVect = append(MinusW, W);
VarL = append(VL, DVL);
NN = length(GB);
OrderGB = newvect(NN);
for (I = 0; I < NN; I++)
OrderGB[I] = ord_w(GB[I], VarL, WVect);
Generator = [];
RGen = [];
for (I = 0; I < NN; I++) {
/*
* B_{S0 - OrderGB[I]}
* i.e. WW[0]*I[0]+...+WW[M-1]*I[M-1] <=
* S0 - OrderGB[I] なる (I[0],...,I[M-1]) の生成
*/
B = gen(WW, S0 - OrderGB[I]);
while (B != []) {
Index = car(B);
D = 1;
for (J = 0; J < M; J++)
D *= DVL[J]^Index[J];
P = weyl_mul(D, GB[I], VarL);
Q = P;
for (J = 0; J < M; J++)
P = subst(P, VL[J], 0);
if (P != 0) {
Generator = cons(P, Generator);
/* inhomogeneous part の保持 */
if ( IH ) RGen = cons(P-Q, RGen);
}
B = cdr(B);
}
}
if ( !NK_RESTRICTION_NOLOG ) {
print("-- fctr(BF) + base :", 2);
tstop();
}
if ( !IH ) return [Generator, Base];
return [Generator, Base, RGen];
}
/* F は s の多項式 */
/* ただし F は Q 上 1 次式に因数分解できるものとする */
def max_int_root(F)
{
L = fctr(F);
/* L の最初の要素は係数なので捨てる */
L = cdr(L);
N = length(L);
S0 = [];
for (I = 0; I < N; I++) {
T = L[I][0];
Root = -coef(T, 0, s) / coef(T, 1, s);
if (dn(Root) == 1) {
if (S0 == []) {
S0 = Root;
} else if (Root > S0) {
S0 = Root;
}
}
}
return S0;
}
def test_max_int_root()
{
print(max_int_root(3*(s+3)^2*(s+1)^3));
print(max_int_root((s+1/2)*(s+1/3)*(s+3)^2*(s+1)^3));
print(max_int_root((s+1/2)*(s+1/3)*(s+1/5)^2*(s+1/10)^3));
}
/*
* I[0]*W[0] + ... + I[M-1]*W[M-1] <= S0 なる
* [I[0], ..., I[M-1]] を返す
*/
def gen(W, S0)
{
if (W == [])
return [[]];
Weight = car(W);
LL = [];
for (I = 0; I <= S0 / Weight; I++) {
L = gen(cdr(W), S0 - Weight * I);
L = map(cons_rev, L, I);
LL = append(L, LL);
}
return LL;
}
/* map 適用のため */
def cons_rev(L, I)
{
return cons(I, L);
}
def test_gen()
{
print(gen([1],3));
print(gen([1,1,1],3));
print(gen([1,1,1],5));
print(gen([1,2,3],5));
}
/* from localb.rr */
def ord_w(P, VL, W)
{
if (P == 0) {
print("ord_w(0)");
return;
}
N = length(VL);
DP = dp_ptod(P, VL);
V = dp_etov(DP);
DP = dp_rest(DP);
for (Max = 0, I = 0; I < N; I++)
Max += V[I] * W[I];
while (DP != 0) {
V = dp_etov(DP);
DP = dp_rest(DP);
for (Weight = 0, I = 0; I < N; I++)
Weight += V[I] * W[I];
if (Max < Weight)
Max = Weight;
}
return Max;
}
/* from localb.rr */
def weyl_mul(P, Q, VL)
{
DP = dp_ptod(P, VL);
DQ = dp_ptod(Q, VL);
return dp_dtop(dp_weyl_mul(DP, DQ), VL);
}
/* restriction のよい計算例はないか... */
def test_restriction()
{
L = [x*dx-1, y*dy-1];
print("restriction");
R1 = restriction(L, [x,y], [dx,dy], [1,0]);
print(R1);
print("sm1.restriction");
R2 = sm1.restriction(L, [x,y], [x]);
print(R2);
}
def test_restriction2()
{
/* from [SST, p202] */
L = [dx^2, dy^2];
print("restriction");
R1 = restriction(L, [x,y], [dx,dy], [1,0]);
print(R1);
print("sm1.restriction");
R2 = sm1.restriction(L, [x,y], [x]);
print(R2);
}
def test_restriction3()
{
/* Appell F_1 */
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
L = [x*(1-x)*dx^2+y*(1-x)*dx*dy+(c-(a+b+1)*x)*dx-b*y*dy-a*b,
y*(1-y)*dy^2+x*(1-y)*dx*dy+(c-(a+bb+1)*y)*dy-bb*x*dx-a*bb];
return restriction(L,[x,y],[dx,dy],[1,0]|param=[a,b,aa,bb,c],s0=S0);
}
/*
* restriction ideal の計算 [SST, Alg5.2.12]
* Id : D-ideal
* VL : var list (x)
* DVL: var list (dx)
* W : weight vector
* VL と長さは同じ
* W[0], ..., W[M] が正で、それ以降の要素は 0 としておく
*
* Option
* param : パラメータとみなす変数のリスト
* s0 : s0 >= 0 のとき generic b-function の計算を行わず
* s-s0 を generic b-function として計算する
* ht : homo, trace-lifting のフラグ
* ht
* 0 = 00(2) 以下 : homo なし, trace なし
* 1 = 01(2) : homo なし, trace あり
* 2 = 10(2) : homo あり, trace なし default (::version() >= 20100415)
* 3 = 11(2) 以上 : homo あり, trace あり
*/
def restriction_ideal(Id, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
if ( ORD && Param )
error("The options `param' and `ord' must be used exclusively.");
L = restriction(Id, VL, DVL, W|inhomo=IH,param=Param,s0=S0);
if ( !L ) return [1];
if ( !NK_RESTRICTION_NOLOG ) tstart();
GG = restriction_ideal_internal(L, VL, DVL, W|inhomo=IH,param=Param,ht=HTrace,ord=ORD);
if ( !NK_RESTRICTION_NOLOG ) {
print("-- restriction_ideal_internal :", 2);
tstop();
}
return GG;
}
def test_restriction_ideal()
{
/* Appell F_1 */
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
L = [x*(1-x)*dx^2+y*(1-x)*dx*dy+(c-(a+b+1)*x)*dx-b*y*dy-a*b,
y*(1-y)*dy^2+x*(1-y)*dx*dy+(c-(a+bb+1)*y)*dy-bb*x*dx-a*bb];
return restriction_ideal(L,[x,y],[dx,dy],[1,0]|inhomo=IH,param=[a,b,aa,bb,c],s0=S0);
}
/* L : result of restriction */
def restriction_ideal_internal(L, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
/*
* restriction module の結果をベクトル表示に直す
* Base と同じ順に dx^Base[I] の係数を I 成分に格納
* restriction の結果の Base は最後の要素が [0,...,0]
* だからベクトル表示した場合、最後の成分が 1 の係数
*/
Gen = L[0];
Base = L[1];
if ( IH ) RGen = L[2];
M = 0;
for (M = 0; M < length(W); M++)
if (W[M] <= 0)
break;
N = length(Base);
P = length(Gen);
GenV = newvect(P);
for (I = 0; I < P; I++) {
T = Gen[I];
V = newvect(N);
for (J = 0; J < N; J++) {
Exp = Base[J];
/* T の dx_1^Exp[0] .. dx_M-1^Exp[M-1] の係数 を求める */
C = T;
for (K = 0; K < M; K++) {
C = coef(C, Exp[K], DVL[K]);
}
V[J] = C;
}
GenV[I] = vtol(V);
}
#ifdef DEBUG_RESTRICTION
print("Generators :");
print(GenV);
#endif
TMP = separate_var(VL,DVL,W|param=Param,ord=ORD);
VY = TMP[0]; DVY = TMP[1]; VLY = append(VY,DVY);
VZ = TMP[2]; DVZ = TMP[3];
OrdM = TMP[4];
GenV = vtol(GenV);
if ( !IH ) {
if ( HTrace <= 0 )
G = nd_weyl_gr(GenV, VLY, 0, [1,OrdM]);
else if ( HTrace == 1 )
G = nd_weyl_gr_trace(GenV, VLY, 0, 0, [1,OrdM]);
else if ( HTrace == 2 )
G = nd_weyl_gr(GenV, VLY, 0, [1,OrdM]|homo=1);
else
G = nd_weyl_gr_trace(GenV, VLY, 1, 0, [1,OrdM]);
if ( G == []) G = [vtol(newvect(N))];
} else {
/* 生成関係式を得るために sygyzy を計算 */
/* OpenXM/src/asir-contrib/testing/noro/module_syz.rr が必要 */
VDVL = append(VL,DVL);
if ( GenV == [] ) GenV = [[0]];
SYZ = newsyz.module_syz(GenV, VLY, 1, OrdM|weyl=1);
G = SYZ[1];
QUO = SYZ[2];
DRGen = map(dp_ptod,RGen,VDVL);
RL = [];
}
/* G の元の中で、最後の成分(i.e. 1 の係数) のみ 0 でないものを取り出す */
Len = length(G);
GG = [];
for (I = 0; I < Len; I++) {
T = G[I];
for (J = 0; J < N - 1; J++)
if (T[J] != 0)
break;
if (J == N - 1 && T[N - 1] != 0) {
if ( !IH )
GG = cons(T[N - 1], GG);
else {
R = inhomo_part(T[N-1],QUO[I],DRGen,VZ,DVZ,VDVL|rest=1);
GG = cons(R[0],GG);
RL = cons(R[1],RL);
}
}
}
if ( GG == [] )
GG = [0];
#ifdef DEBUG_RESTRICTION
print("GG :");
print(GG);
#endif
if ( !IH )
return GG;
if ( RL == [] )
RL = [zero_inhomo(VZ)];
return [GG,RL];
}
def test_restriction_ideal_internal()
{
L = restriction([x*dx-1,y*dy-1],VL=[x,y],DVL=[dx,dy],W=[1,0]);
print(L);
restriction_ideal_internal(L, VL, DVL, W);
}
def test_nd_weyl_gr()
{
/* internal error
nd_weyl_gr([newvect(2,[1,x]),newvect(2,[x,y])], [x,y,dx,dy],0,[1,0]);
*/
/* invarid argument
nd_weyl_gr(newvect(2,[[1,x],[x,y]]), [x,y,dx,dy], 0, [1,0]);
*/
/* correct */
nd_weyl_gr([[1,x],[x,y]],[x,y,dx,dy],0,[1,0]);
}
/* 微分作用素 P(x,dx) を Fourier 変換したものを返す */
/* i.e. x_i -> - dx_i, dx_i -> x_i */
def fourier_trans(P, VL, DVL)
{
N = length(VL);
VLL = append(VL, DVL);
DP = dp_ptod(P, VLL);
S = 0;
VE = newvect(2 * N);
while (DP != 0) {
C = dp_hc(DP);
V = dp_etov(DP);
/* x^a -> (-1)^|a| dx^a に変換(dp で) */
M = 0;
for (I = 0; I < N; I++)
VE[I] = 0;
for (I = 0; I < N; I++) {
VE[N + I] = V[I];
M += V[I];
}
DPD = dp_vtoe(VE);
if (M % 2 == 1)
DPD = -DPD;
/* dx^a -> x^a に変換 */
for (I = 0; I < N; I++)
VE[I] = V[N + I];
for (I = 0; I < N; I++)
VE[N + I] = 0;
DPX = dp_vtoe(VE);
T = C * dp_weyl_mul(DPD, DPX);
S += T;
DP = dp_rest(DP);
}
return dp_dtop(S, VLL);
}
def test_fourier_trans()
{
print(fourier_trans(x*dx, [x], [dx]));
print(fourier_trans(x^2*dx^2, [x], [dx]));
print(fourier_trans(x*dx*y+1, [x], [dx]));
}
/* 微分作用素 P(x,dx) を Fourier 変換したものを返す */
/* i.e. x_i -> dx_i, dx_i -> - x_i */
def inv_fourier_trans(P, VL, DVL)
{
N = length(VL);
VLL = append(VL, DVL);
DP = dp_ptod(P, VLL);
S = 0;
VE = newvect(2 * N);
while (DP != 0) {
C = dp_hc(DP);
V = dp_etov(DP);
/* x^a -> dx^a に変換(dp で) */
M = 0;
for (I = 0; I < N; I++)
VE[I] = 0;
for (I = 0; I < N; I++)
VE[N + I] = V[I];
DPD = dp_vtoe(VE);
/* dx^a -> (-1)^|a| x^a に変換 */
for (I = 0; I < N; I++) {
VE[I] = V[N + I];
M += V[N + I];
}
for (I = 0; I < N; I++)
VE[N + I] = 0;
DPX = dp_vtoe(VE);
if (M % 2 == 1)
DPX = -DPX;
T = C * dp_weyl_mul(DPD, DPX);
S += T;
DP = dp_rest(DP);
}
return dp_dtop(S, VLL);
}
def test_inv_fourier_trans()
{
print(inv_fourier_trans(fourier_trans(x^2*dx^2*y+x*dx*dy+y, [x], [dx]),[x],[dx]));
}
/*
* integration module の計算([SST, Alg5.5.4])
* Id : ideal
* VL : var list (x)
* DVL: var list (dx)
* W : weight vector
* VL と長さは同じ
* W[0], ..., W[M] が正で、それ以降の要素は 0 としておく
*
* Option
* param : パラメータとみなす変数のリスト
* s0 : s0 >= 0 のとき generic b-function の計算を行わず
* s-s0 を generic b-function として計算する
*/
def integration(Id, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
for (M = 0; M < length(W); M++)
if (W[M] <= 0)
break;
VLM = newvect(M);
DVLM = newvect(M);
for (I = 0; I < M; I++) {
VLM[I] = VL[I];
DVLM[I] = DVL[I];
}
VLM = vtol(VLM);
DVLM = vtol(DVLM);
FId = map(fourier_trans, Id, VLM, DVLM);
R = restriction(FId, VL, DVL, W|inhomo=IH,param=Param, s0=S0);
if ( !R ) return 0;
#ifdef DEBUG_RESTRICTION
print("restriction of FId :");
print(R);
#endif
Generator = map(inv_fourier_trans, R[0], VLM, DVLM);
if ( !IH ) return [Generator, R[1]];
return [Generator, R[1], R[2]];
}
def test_integration()
{
/* [SST, Ex5.5.2, 5.5.6] */
/* integration of the annihilating ideal of 1/(t^2-x) w.r.t. t */
Id = [2*t*dx+dt, t*dt+2*x*dx+2];
return integration(Id, [t,x],[dt,dx],[1,0]);
}
def test_integration2()
{
/* [SST, Ex5.6.13] */
/* integration of the annihilating ideal of e^(t^3-x*t) w.r.t. t */
Id = [dt - (3*t^2-x), dx + t];
return integration(Id, [t,x],[dt,dx],[1,0]);
}
def test_integration3()
{
/* [SST, Ex5.5.9] */
/* 計算が合わない?? */
F = x^3+t1^3+t2^3;
Ann = ann(F);
Ann0 = map(subst, Ann, s, -2);
Int = integration(Ann0, [t1,t2,x],[dt1,dt2,dx],[1,1,0]);
return Int;
}
def test_integration4()
{
/* [SST, Ex5.5.16] */
/* test_integration3 と同じイデアルだが、重みベクトルが異なる */
F = x1^3+x2^3+x3^3;
Ann = ann(F);
Ann0 = map(subst, Ann, s, -2);
Int = integration(Ann0, [x1,x2,x3],[dx1,dx2,dx3],[1,2,3]);
return Int;
}
def test_integration5()
{
/* [SST, Ex5.6.4] */
G = x1*t1^2*t2+x2*t1^2*t2^2+x3*t1*t2^2+x4+x5*t1*t2;
Ann = ann(G);
/* G の bfct は s + 1 */
Ann0 = map(subst, Ann, s, -1);
Int = integration(Ann0, [t1,t2,x1,x2,x3,x4,x5],[dt1,dt2,dx1,dx2,dx3,dx4,dx5],[1,1,0,0,0,0,0]);
return Int;
}
/*
* integration ideal の計算 [SST, Alg5.5.4]
* Id : D-ideal
* VL : var list (x)
* DVL: var list (dx)
* W : weight vector
* VL と長さは同じ
* W[0], ..., W[M] が正で、それ以降の要素は 0 としておく
*
* Option
* param : パラメータとみなす変数のリスト
* s0 : s0 >= 0 のとき generic b-function の計算を行わず
* s-s0 を generic b-function として計算する
* inhomo: 0 でないとき inhomogeneous part の情報もあわせて返す
* ht : homo, trace-lifting のフラグ
* ht
* 0 = 00(2) 以下 : homo なし, trace なし
* 1 = 01(2) : homo なし, trace あり
* 2 = 10(2) : homo あり, trace なし default (::version() >= 20100415)
* 3 = 11(2) 以上 : homo あり, trace あり
*/
def integration_ideal(Id, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
if ( ORD && Param )
error("The options `param' and `ord' must be used exclusively.");
L = integration(Id, VL, DVL, W|inhomo=IH,param=Param, s0=S0);
if ( !L ) return [1];
#ifdef DEBUG_RESTRICTION
print("integration :");
print(L);
#endif
if ( !NK_RESTRICTION_NOLOG ) tstart();
GG = integration_ideal_internal(L, VL, DVL, W|inhomo=IH,param=Param,ht=HTrace, ord=ORD);
if ( !NK_RESTRICTION_NOLOG ) {
print("-- integration_ideal_internal :", 2);
tstop();
}
return GG;
}
/* L : result of integration */
def integration_ideal_internal(L, VL, DVL, W)
{
if ( type(IH=getopt(inhomo)) == -1) IH = 0;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
/*
* integration module の結果をベクトル表示に直す
* Base と同じ順に x^Base[I] の係数を I 成分に格納
* restriction の結果の Base は最後の要素が [0,...,0]
* だからベクトル表示した場合、最後の成分が 1 の係数
*/
Gen = L[0];
Base = L[1];
if ( IH ) RGen = L[2];
M = 0;
for (M = 0; M < length(W); M++)
if (W[M] <= 0)
break;
N = length(Base);
P = length(Gen);
GenV = newvect(P);
for (I = 0; I < P; I++) {
T = Gen[I];
V = newvect(N);
for (J = 0; J < N; J++) {
Exp = Base[J];
/* T の x_1^Exp[0] .. x_M-1^Exp[M-1] の係数 を求める */
C = T;
for (K = 0; K < M; K++) {
C = coef(C, Exp[K], VL[K]);
}
V[J] = C;
}
GenV[I] = vtol(V);
}
#ifdef DEBUG_RESTRICTION
print("Generators :");
print(GenV);
#endif
TMP = separate_var(VL,DVL,W|param=Param,ord=ORD);
VY = TMP[0]; DVY = TMP[1]; VLY = append(VY,DVY);
VZ = TMP[2]; DVZ = TMP[3];
OrdM = TMP[4];
GenV = vtol(GenV);
if ( !IH ) {
if ( HTrace <= 0 )
G = nd_weyl_gr(GenV, VLY, 0, [1,OrdM]);
else if ( HTrace == 1 )
G = nd_weyl_gr_trace(GenV, VLY, 0, 0, [1,OrdM]);
else if ( HTrace == 2 )
G = nd_weyl_gr(GenV, VLY, 0, [1,OrdM]|homo=1);
else
G = nd_weyl_gr_trace(GenV, VLY, 1, 0, [1,OrdM]);
if ( G == []) G = [vtol(newvect(N))];
} else {
/* 生成関係式を得るために sygyzy を計算 */
/* OpenXM/src/asir-contrib/testing/noro/module_syz.rr が必要 */
VDVL = append(VL,DVL);
if ( GenV == [] ) GenV = [[0]];
SYZ = newsyz.module_syz(GenV, VLY, 1, OrdM|weyl=1);
G = SYZ[1];
QUO = SYZ[2];
DRGen = map(dp_ptod,RGen,VDVL);
RL = [];
}
#if 0
if (ParamType != 4) {
NN = length(VL) - M;
VLY = newvect(NN * 2);
for (I = 0; I < NN; I++) {
VLY[I] = VL[I + M];
VLY[I + NN] = DVL[I + M];
}
VLY = vtol(VLY);
GenV = vtol(GenV);
G = nd_weyl_gr(GenV, VLY, 0, [1,0]);
} else {
/*
* VLY
* [x, param, dx, dparam]
*/
NN = length(VL) - M + length(Param);
VLY = newvect(NN * 2);
for (I = 0; I < length(VL) - M; I++) {
VLY[I] = VL[I + M];
VLY[I + NN] = DVL[I + M];
}
for (;I < NN; I++) {
VLY[I] = Param[I - length(VL) + M];
VLY[I + NN] = strtov("d" + rtostr(Param[I - length(VL) + M]));
}
VLY = vtol(VLY);
GenV = vtol(GenV);
/*
* OrdM
* x,param,dx,dparam
* [1, 0, 1, 0]
* [ grevlex ]
*/
OrdM = newmat(2 * NN + 2, 2 * NN);
for (I = 0; I < length(VL) - M; I++) {
OrdM[0][I] = 1;
OrdM[0][I + NN] = 1;
}
for (I = 0; I < 2 * NN; I++)
OrdM[1][I] = 1;
for (I = 0; I < 2 * NN; I++)
OrdM[I][2 * NN - 1 - I] = -1;
G = nd_weyl_gr(GenV, VLY, 0, [1,OrdM]);
}
#endif
/* G の元の中で、最後の成分(i.e. 1 の係数) のみ 0 でないものを取り出す */
Len = length(G);
GG = [];
for (I = 0; I < Len; I++) {
T = G[I];
for (J = 0; J < N - 1; J++)
if (T[J] != 0)
break;
if (J == N - 1 && T[N - 1] != 0) {
if ( !IH )
GG = cons(T[N - 1], GG);
else {
R = inhomo_part(T[N-1],QUO[I],DRGen,VZ,DVZ,VDVL);
GG = cons(R[0],GG);
RL = cons(R[1],RL);
}
}
}
if ( GG == [] )
GG = [0];
#ifdef DEBUG_RESTRICTION
print("GG :");
print(GG);
#endif
if ( !IH )
return GG;
if ( RL == [] )
RL = [zero_inhomo(DVZ)];
return [GG,RL];
}
/* nd_weyl_gr を実行せず、加群の生成元を返す */
/* 計算は完了していない */
def integration_ideal_(Id, VL, DVL, W)
{
L = integration(Id, VL, DVL, W);
LL = integration_ideal_internal_(L, VL, DVL, W);
return LL;
}
/* nd_weyl_gr を実行せず、加群の生成元を返す */
def integration_ideal_internal_(L, VL, DVL, W)
{
/*
* integration module の結果をベクトル表示に直す
* Base と同じ順に x^Base[I] の係数を I 成分に格納
* restriction の結果の Base は最後の要素が [0,...,0]
* だからベクトル表示した場合、最後の成分が 1 の係数
*/
Gen = L[0];
Base = L[1];
M = 0;
for (M = 0; M < length(W); M++)
if (W[M] <= 0)
break;
N = length(Base);
P = length(Gen);
GenV = newvect(P);
for (I = 0; I < P; I++) {
T = Gen[I];
V = newvect(N);
for (J = 0; J < N; J++) {
Exp = Base[J];
/* T の x_1^Exp[0] .. x_M-1^Exp[M-1] の係数 を求める */
C = T;
for (K = 0; K < M; K++) {
C = coef(C, Exp[K], VL[K]);
}
V[J] = C;
}
GenV[I] = vtol(V);
}
#ifdef DEBUG_RESTRICTION
print("Generators :");
print(GenV);
#endif
NN = length(VL) - M;
VLY = newvect(NN * 2);
for (I = 0; I < NN; I++) {
VLY[I] = VL[I + M];
VLY[I + NN] = DVL[I + M];
}
VLY = vtol(VLY);
GenV = vtol(GenV);
return [GenV, VLY];
}
def test_integration_ideal()
{
/* [SST, Ex5.5.2, 5.5.6] */
/* integration of the annihilating ideal of 1/(t^2-x) w.r.t. t */
Id = [2*t*dx+dt, t*dt+2*x*dx+2];
L = integration(Id, VL=[t,x],DVL=[dt,dx],W=[1,0]);
print("integration module :");
print(L);
R = integration_ideal_internal(L, VL, DVL, W);
print("integration ideal :");
print(R);
}
def test_integration_ideal2()
{
/* [SST, Ex5.5.7] */
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
F = (x1+x2*t+x3*t^2+x4*t^3)*t^2;
BF = bfct(F);
print("BF :");
print(BF);
Ann = ann(F);
Ann0 = map(subst, Ann, s, -1/5);
print("Ann(F^s)|s=(-1/5) :");
print(Ann0);
print("integration ideal :");
Int = integration_ideal(Ann0, [t,x1,x2,x3,x4],[dt,dx1,dx2,dx3,dx4],[1,0,0,0,0]|inhomo=IH);
print(Int);
}
def test_integration_ideal5()
{
/* [SST, Ex5.6.4] */
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
G = x1*t1^2*t2+x2*t1^2*t2^2+x3*t1*t2^2+x4+x5*t1*t2;
Ann = ann(G);
/* G の bfct は s + 1 */
Ann0 = map(subst, Ann, s, -1);
Int = integration_ideal(Ann0, [t1,t2,x1,x2,x3,x4,x5],[dt1,dt2,dx1,dx2,dx3,dx4,dx5],[1,1,0,0,0,0,0]|inhomo=IH);
return Int;
}
def test_integration_ideal5_()
{
/* [SST, Ex5.6.4] */
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
G = x1*t1^2*t2+x2*t1^2*t2^2+x3*t1*t2^2+x4+x5*t1*t2;
Ann = ann(G);
/* G の bfct は s + 1 */
Ann0 = map(subst, Ann, s, -1);
Int = integration_ideal_(Ann0, [t1,t2,x1,x2,x3,x4,x5],[dt1,dt2,dx1,dx2,dx3,dx4,dx5],[1,1,0,0,0,0,0]|inhomo=IH);
return Int;
}
def test_integration_ideal6()
{
/*
* Gauss hypergeometric integration
* L = ndbf.ann_n([t,1-t,1-x*t]);
* L = map(subst, L, s0, b-1, s1, c-b-1, s2, -a);
*/
if (type(IH=getopt(inhomo)) == -1) IH = 0;
if (type(S0=getopt(s0)) == -1) S0 = -1;
L = [(-x^2+x)*dx^2+((-t+1)*dt+(-a-b-1)*x+c-1)*dx-b*a,(-t+1)*x*dx+(t^2-t)*dt+(-c+2)*t+b-1,(t*x-1)*dx+a*t];
Int = integration_ideal(L,[t,x],[dt,dx],[1,0]|param=[a,b,c],s0=S0,inhomo=IH);
return Int;
}
def test_integration_ideal7()
{
/* A-hypergeometric */
/* \int (x11+x21*t)^a1 (x12+x22*t)^a2 t^c dt */
if (type(IH=getopt(inhomo)) == -1) IH = 0;
if (type(S0=getopt(s0)) == -1) S0 = -1;
L = [-x11*dx11-dx21*x21+a1,-x12*dx12-dx22*x22+a2,-t*dt+c+dx21*x21+dx22*x22,t*dx12-dx22,t*dx11-dx21,dx21*dx12-dx22*dx11];
Int = integration_ideal(L,[t,x11,x12,x21,x22],[dt,dx11,dx12,dx21,dx22],[1,0,0,0,0]|param=[a1,a2,c],s0=S0,inhomo=IH);
return Int;
}
def test_integration_ideal8()
{
ORD = getopt(ord);
F = a*x^2-2*b*x+c;
A = ann(F);
AN = map(subst,A,s,-(-p*dp-1));
MB = cons(p*F-1,AN);
V = [x,p,a,b,c];
DV = [dx,dp,da,db,dc];
W = [1,0,0,0,0];
if ( type(ORD) == -1 ) {
ORD = newmat(9,8,[[0,0,0,0,0,1,1,1],[1,1,1,1,1,1,1,1]]);
for ( I = 0; I < 7; I++ )
ORD[2+I][7-I] = -1;
}
INT = integration_ideal(MB,V,DV,W|inhomo=1,ord=ORD);
return INT;
}
/* ann_mul */
/*
* x_i --> x_i
* y_i --> x_i - y_i
* dx_i --> dx_i + dy_i
* dy_i --> - dy_i
*/
def coord_conv(P, VL1, DVL1, VL2, DVL2)
{
N = length(VL1);
for (I = 0; I < N; I++)
P = subst(P, VL2[I], VL1[I] - VL2[I]);
for (I = 0; I < N; I++)
P = subst(P, DVL2[I], -DVL2[I]);
for (I = 0; I < N; I++)
P = subst(P, DVL1[I], DVL1[I] + DVL2[I]);
return P;
}
def check_coord_conv()
{
return coord_conv((1-tt)*dtt,[t],[dt],[tt],[dtt]);
}
def check_coord_conv2()
{
return coord_conv(t*dt,[t],[dt],[tt],[dtt]);
}
/* variable a --> a_ */
def dummy_var(A)
{
return strtov(rtostr(A) + "_");
}
def ann_mul(Id1, Id2, VL, DVL)
{
if (type(Param = getopt(param)) == -1)
Param = 0;
if (type(S0 = getopt(s0)) == -1)
S0 = -1;
VLL = map(dummy_var, VL);
DVLL = map(dummy_var, DVL);
Id2 = map(subst_vect, Id2, VL, VLL);
Id2 = map(subst_vect, Id2, DVL, DVLL);
Id1 = map(coord_conv, Id1, VL, DVL, VLL, DVLL);
Id2 = map(coord_conv, Id2, VL, DVL, VLL, DVLL);
Id = append(Id1, Id2);
V = append(VLL, VL);
DV = append(DVLL, DVL);
W = newvect(length(V));
for (I = 0; I < length(V)/2; I++)
W[I] = 1;
W = vtol(W);
if ( !NK_RESTRICTION_NOLOG ) {
print([Id, V, DV, W]);
print("option:" + rtostr([Param, S0]));
}
return restriction_ideal(Id, V, DV, W|param = Param, s0 = S0);
}
def test_ann_mul()
{
return ann_mul([t*dt], [(1-t)*dt],[t],[dt]);
}
def test_ann_mul2()
{
F = x^3+y^2;
G = x+y^2;
Ann1 = map(subst, ann(F), s, -1);
Ann2 = map(subst, ann(G), s, -1);
Ann = map(subst, ann(F*G), s, -1);
AnnMul = ann_mul(Ann1, Ann2, [x,y],[dx,dy]);
return [AnnMul,Ann];
}
def test_ann_mul3()
{
/* t^{a-1} (1-t)^{c-a-1} の ann */
I12 = ndbf.ann_n([t,1-t]);
I12 = map(subst, I12, s0, a-1, s1, c-a-1);
/* I12 は C<t,z,dt,dz> のイデアルとみなすので dz も追加 */
I12 = cons(dz, I12);
print("I12 :"); print(I12);
/* e^{tz} の ann */
I3 = [dt-z,dz-t];
I123 = ann_mul(I12, I3, [t,z], [dt,dz]|param=[a,c],s0=0);
print("I12 :"); print(I12);
return integration_ideal(I123, [t,z], [dt,dz], [1,0]|param=[a,c],s0=0);
}
/* F の VL[0] に V[0], ... を代入 */
def subst_vect(F, VL, V)
{
N = length(VL);
for (I = 0; I < N; I++)
F = subst(F, VL[I], V[I]);
return F;
}
/*----------------------------------------------------------------- */
/*-------------- test data for integration ideal ------------------ */
/*----------------------------------------------------------------- */
/*
#define DEBUG_FANO
*/
/* 2 次元 Fano 多面体で頂点の座標が 0,1,-1 のいづれかのもの + 原点 */
def fano2(I)
{
Fano2 = [
[[-1,0],[0,1],[1,-1],[0,0]],
[[-1,0],[0,-1],[0,1],[1,0],[0,0]],
[[-1,0],[-1,1],[0,1],[1,-1],[0,0]],
[[-1,0],[-1,1],[0,-1],[0,1],[1,-1],[0,0]],
[[-1,0],[-1,1],[0,-1],[0,1],[1,-1],[1,0],[0,0]]];
return Fano2[I];
}
/* 3 次元 Fano 多面体で頂点の座標が 0,1,-1 のいづれかのもの + 原点 */
def fano3(I)
{
Fano3 = [
[[-1,0,0],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[0,0,-1],[0,0,1],[0,1,1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,1,1],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[0,0,-1],[0,0,1],[0,1,-1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,0,1],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[0,-1,1],[0,0,-1],[0,0,1],[0,1,-1],[1,0,0],[0,0,0]],
[[-1,0,0],[-1,1,0],[-1,1,1],[0,0,1],[0,1,1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[0,-1,1],[0,0,-1],[0,0,1],[0,1,0],[1,0,-1],[0,0,0]],
[[-1,0,0],[-1,1,1],[0,0,-1],[0,0,1],[0,1,1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,0,-1],[-1,1,1],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,0,-1],[0,0,-1],[0,0,1],[0,1,-1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,1,1],[0,-1,-1],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,1,1],[0,-1,-1],[0,0,-1],[0,0,1],[0,1,1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,1,0],[-1,1,1],[0,-1,-1],[0,0,1],[0,1,1],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,0,1],[0,-1,1],[0,0,-1],[0,0,1],[0,1,0],[1,0,-1],[0,0,0]],
[[-1,0,0],[-1,0,1],[-1,1,1],[0,-1,0],[0,0,1],[0,1,0],[1,-1,-1],[0,0,0]],
[[-1,0,0],[-1,1,1],[1,-1,-1],[0,0,-1],[0,0,1],[0,1,1],[1,-1,-1],[1,0,0],[0,0,0]],
[[-1,0,0],[-1,0,1],[-1,1,1],[0,-1,0],[0,0,1],[0,1,0],[1,-1,-1],[1,0,-1],[0,0,0]]];
return Fano3[I];
}
/* L : 点データのリスト ex. fano2, fano3 の返り値 */
def p_to_mat(L)
{
N = length(L);
M = length(L[0]);
A = newmat(M + 1, N);
for (J = 0; J < N; J++)
A[0][J] = 1;
for (J = 0; J < N; J++)
for (I = 0; I < M; I++)
A[I + 1][J] = L[J][I];
return A;
}
def test_p_to_mat()
{
L = fano2(0);
print(L);
print(p_to_mat(L));
L = fano3(10);
print(L);
print(p_to_mat(L));
}
/* fano polytope に付随する GKZ 系の行列 A を返す */
/* D : dimension (2 or 3), I : data index (0--4(D=2), 0--17(D=3)) */
def fano_polytope_a(D, I)
{
if (D == 2) {
L = fano2(I);
} else if (D == 3) {
L = fano3(I);
} else {
print("D is 2 or 3");
return;
}
return p_to_mat(L);
}
def test_fano_polytope_a()
{
print(fano_polytope_a(2, 0));
print(fano_polytope_a(3, 11));
}
def fano_polytope(D, I)
{
A = fano_polytope_a(D, I);
P = a_to_poly(A);
return append([A], P);
}
/*
* 行列 A から多項式を生成
* p_to_mat(fano*(I)) が A
* ex. A = [[1,1,1],[-1,0,1],[0,1,-1]]
* F = x1*t1^-1 + x2*t2 + x3*t1*t2^-1
* t1*t2*F = x1*t2 x2*t1*t2^2 + x3*t1^2 <-- これを返す
*/
def a_to_poly(A)
{
M = size(A)[0];
N = size(A)[1];
XV = var_list("x", 1, N);
TV = var_list("t", 1, M-1);
S = 0;
V = newvect(M - 1);
for (J = 0; J < N; J++) {
for (I = 0; I < M - 1; I++) {
V[I] = A[I + 1][J];
V[I]++;
}
DPE = dp_vtoe(V);
E = dp_dtop(DPE, TV);
S += E * XV[J];
}
VL = append(TV, XV);
DXV = map(d_var, XV);
DTV = map(d_var, TV);
DVL = append(DTV, DXV);
return [S, VL, DVL];
}
def test_a_to_poly()
{
A = newmat(3,3,[[1,1,1],[-1,0,1],[0,1,-1]]);
print(A);
print(a_to_poly(A));
A = p_to_mat(fano3(10));
print(A);
print(a_to_poly(A));
}
def var_list(Str, From, To)
{
L = [];
for (I = From; I <= To; I++)
L = cons(strtov(Str + rtostr(I)), L);
return reverse(L);
}
def test_var_list()
{
print(var_list("a", 1, 11));
print(var_list("tt", 0, 5));
}
/* 不定元 Var に "d" を追加した不定元を生成 */
def d_var(Var)
{
return strtov("d" + rtostr(Var));
}
def int_fano(D, I)
{
L = fano_polytope(D, I);
F = L[1];
VL = L[2];
DVL = L[3];
W = newvect(length(VL));
for (I = 0; I < D; I++)
W[I] = 1;
W = vtol(W);
print("F :", 2);
print(F, 2);
tstart();
Id = ann(F);
print("-- ann :", 2);
tstop();
print("Ann F^s : ", 2);
print(Id);
#ifdef DEBUG_FANO
/* b-関数の計算が重たいので skip */
tstart();
BF = bfunction(F);
print("-- bfunction :", 2);
tstop();
print("bfct : ");
print(fctr(BF));
#endif
Id0 = map(subst, Id, s, -1);
print("Ann F^(-1) : ", 2);
print(Id0, 2);
J = integration_ideal(Id0, VL, DVL, W);
return J;
}
/**** added functions for inhomogeneous ****/
/*
* J = (I + dt_1 D + dt_2 D + ... + dt_m D) \cap D_Y
*
* Input
* TT : 積分イデアル J の生成元 (content が取れる可能性あり)
* QI : 生成関係式の係数 i.e TT = \sum QI[I]*(GenV[I]の最後の成分)
* DRGen : Rest part
* VZ : [ t_1, t_2, ..., t_m]
* DVZ : [dt_1, dt_2, ..., dt_m]
*
* Output [G,R]
* G[K] : 積分イデアル J の生成元 (content = 1)
* R[K] : [[[dt_1, ***],[dt_2, ***],...,[dt_m, ***]], 分母]
* i.e.
* G[K]-(dt_1 R[K][0][0][1] + dt_2 R[K][0][1][1]
* + ... + dt_m R[K][0][m][1])/R[K][1] \in I
*
*/
def inhomo_part(TT,QI,DRGen,VZ,DVZ,VDVL)
{
if ( type(REST=getopt(rest)) == -1 ) REST = 0;
#if 0
TMP = cont(TT);
CT = dp_hc(dp_ptod(TMP,vars(TMP)));
G = TT/CT;
#else
G = ptozp(TT);
CT = tdiv(TT,G);
#endif
DQI = map(dp_ptod,QI,VDVL);
for ( DP = 0, K = 0; K < length(DQI); K++ )
DP += dp_weyl_mul(DQI[K],DRGen[K]);
P = dp_dtop(DP,VDVL);
GCD = 1;
if ( P ) {
#if 0
TMP = cont(P);
CP = dp_hc(dp_ptod(TMP,vars(TMP)));
#else
CP = tdiv(P,ptozp(P));
#endif
GCD = igcd(CT,CP);
P /= GCD;
}
R = [];
for ( I = 0; I < length(VZ); I++ ) {
Q = subst(P,VZ[I],0);
S = sdiv(P-Q,VZ[I]);
if ( REST )
R = cons([VZ[I],S],R);
else
R = cons([DVZ[I],inv_fourier_trans(S,VZ,DVZ)],R);
P = Q;
}
if ( P )
error("inhomo_part : invalid inhomogeneous operator");
R = [reverse(R),CT/GCD];
return [G,R];
}
def trans_inhomo(P,INT,V,DV,W)
{
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
TMP = separate_var(V,DV,W);
VX = TMP[0]; DVX = TMP[1]; VDVX = append(VX,DVX);
VY = TMP[2]; DVY = TMP[3]; VDVY = append(VY,DVY);
VDV = append(V,DV);
J = INT[0]; IH = INT[1];
IND = car(car(IH));
for ( CJ = [], CIH = []; J != []; J = cdr(J), IH =cdr(IH) ) {
A = car(IH);
CJ = cons(A[1]*car(J),CJ);
CIH = cons(A[0],CIH);
}
SYZ = newsyz.module_syz(cons(P,CJ),VDVX,1,ORD|weyl=1);
GR = pick_gen_relation(SYZ[0],VDVX);
DR = vtol(newvect(length(IND)));
CONST = car(GR);
for ( GR = cdr(GR); GR != []; GR = cdr(GR), CIH = cdr(CIH) ) {
if ( HGR = car(GR) ) {
DH = dp_ptod(HGR,VDV);
DR = trans_inhomo_sub(DH,DR,car(CIH),VDV);
}
}
RR = map(dp_dtop,DR,VDV);
for ( GCD = CONST, TR = RR; TR != []; TR =cdr(TR) ) {
if ( TT = car(TR) ) {
CC = tdiv(TT,ptozp(TT));
GCD = igcd(GCD,CC);
}
}
for ( R = []; RR != []; RR = cdr(RR), IND = cdr(IND) )
R = cons([car(car(IND)),car(RR)/GCD],R);
return [reverse(R),-CONST/GCD];
}
def trans_inhomo_sub(DH,DR,IH,VDV)
{
for ( L = []; IH != []; IH = cdr(IH), DR = cdr(DR) )
L = cons(dp_weyl_mul(DH,dp_ptod(car(IH)[1],VDV)) + car(DR),L);
return reverse(L);
}
def pick_gen_relation(GR,V)
{
for ( ; GR != []; GR = cdr(GR) )
if ( is_nonzero_const(car(GR)[0],V) )
return car(GR);
error("The input is not in integration/restriction ideal.");
}
def is_nonzero_const(A,V)
{
if ( !A )
return 0;
for ( VA = vars(A); VA != []; VA = cdr(VA) )
if ( member(car(VA),V) )
return 0;
return 1;
}
def separate_var(VL,DVL,W)
{
DUMMY = dummy_var;
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(ORD=getopt(ord)) == -1 ) ORD = 0;
if ( type(Weight=getopt(weight)) == -1 ) Weight = 0;
VY = []; DVY = [];
VZ = []; DVZ = [];
W1 = []; W2 = [];
for ( I = 0; I < length(VL); I++ ) {
if ( W[I] ) {
VZ = cons(VL[I],VZ);
DVZ = cons(DVL[I],DVZ);
W1 = cons(W[I],W1);
} else {
VY = cons(VL[I],VY);
DVY = cons(DVL[I],DVY);
W2 = cons(W[I],W2);
}
}
OrdM = ORD;
if ( !Weight )
if ( VY == [] ) {
VY = [DUMMY];
DVY = [strtov("d"+rtostr(DUMMY))];
}
if ( Param ) {
DParam =[];
for ( I = 0; I < length(Param); I++ )
DParam = cons(strtov("d"+rtostr(Param[I])),DParam);
OrdM = make_elim_order_mat(length(VY),length(Param));
DVY = append(DParam,DVY);
VY = append(reverse(Param),VY);
}
R = append(map(reverse,[VY,DVY,VZ,DVZ]),[OrdM]);
if ( Weight )
return append(R,map(reverse,[W2,W1]));
return R;
}
/* 1 行目の左から順に 1 が M 個, 0 が N 個, 1が M 個, 0 が N 個 並び,
* 下は grevlex 順を表す行列を返す
*/
def make_elim_order_mat(M,N)
{
NN = M+N;
OrdM = newmat(2*NN+1,2*NN);
for ( I = 0; I < M; I++ ) {
OrdM[0][I] = 1;
OrdM[0][I+NN] = 1;
}
for ( I = 0; I < 2*NN; I++ )
OrdM[1][I] = 1;
for ( I = 2; I < 2*NN+1 ; I++ )
OrdM[I][2*NN-I+1] = -1;
return OrdM;
}
def make_weight_mat(W)
{
N = length(W);
M = newmat(N+1,N,[W]);
for ( I = 0; I < N; I++)
M[1][I] = 1;
for ( I = 0; I < N-1; I++)
M[I+2][N-I-1] = -1;
return M;
}
/* P : a operator in D<V,DV>
* CE : C*exp(E) を表すリスト [C,E]
*
* C*exp(E) に P を左から作用させた結果を返す
*/
def action_exp(P,V,DV,CE)
{
VDV = append(V,DV);
DP = dp_ptod(P,VDV);
N = length(V);
R = 0;
for ( T = DP; T; T = dp_rest(T) )
R += action_exp_1(dp_hm(T),N,V,CE);
return [R,CE[1]];
}
def action_exp_1(M,N,V,CE)
{
C = CE[0];
E = CE[1];
CO = dp_hc(M);
EX = dp_etov(M);
for ( I = 0; I < N; I++ ) {
VI = V[I];
CO *= VI^EX[I];
DEVI = diff(E,VI);
for ( J = 0, EXI = EX[I+N]; J < EXI; J++ )
C = diff(C,VI)+C*DEVI;
}
return CO*C;
}
/* A \subset B ? */
def subset(A,B)
{
while ( A != [] ) {
if ( !member(car(A),B) )
return 0;
A = cdr(A);
}
return 1;
}
def zero_inhomo(VL)
{
R = [];
for ( I = 0; I < length(VL); I++ )
R = cons([VL[I],0],R);
R = reverse(R);
RL = [R,1];
return RL;
}
/* \int C*exp(E) dt */
def integration_exp(C,E,V,DV,W)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
Id = [];
for ( I = 0; I < length(V); I++ ) {
T = red(diff(C,V[I])/C);
NT = nm(T); DT = dn(T);
Id = cons(DT*DV[I]-NT-DT*diff(E,V[I]),Id);
}
G = integration_ideal(Id,V,DV,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
return G;
}
def test_integration_exp()
{
/* CREST GB School 演習問題 2-(c) */
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
return integration_exp(1,-t-x*t^3,[t,x],[dt,dx],[1,0]|inhomo=IH);
}
def time_print(T)
{
T1 = time();
CPU = time_round(T1[0] - T[0],4);
GC = time_round(T1[1] - T[1],4);
ESP = time_round(T1[3] - T[3],4);
P = rtostr(CPU)+"sec";
if ( GC )
P = P + " + gc : " + rtostr(GC)+ "sec";
P = P + "(" + rtostr(ESP) + "sec)";
if ( !NK_RESTRICTION_NOLOG ) print(P);
return 1;
}
def time_round(Z,N)
{
if ( Z == 0 ) return 0;
B = 10^N;
if ( Z >= B ) {
for ( I = 0; Z > B; I++, Z = deval(Z/10) );
Z = rint(Z);
Z = deval(Z*10^I);
/*
Z = deval(rint(Z));
*/
} else {
B /= 10;
for ( I = 0; Z < B; I++, Z *= 10 );
Z = rint(Z);
Z = deval(Z/10^I);
}
return Z;
}
/*
* Input
* Id : an ideal of D<t,dt,x,dx> which annihilates u(x,t)
* V : var list (t,x) : [ t_1,..., t_m, x_1,..., x_n]
* DV : var list (dt,dx) : [dt_1,...,dt_m,dx_1,...,dx_n]
* W : weight : W[I] > 0 ( I < m ), W[I] = 0 ( I >= m )
* LB : list of lower bounds [a_1,...,a_m]
* UB : list of upper bounds [b_1,...,b_m]
* a_i, b_i : number, strings "inf" or "-inf", parameter
*
* Output
* ideal J \subset D<x,dx> s.t.
* J \bullet \int_{a_1}^{b_1} ... \int_{a_m}^{b_m} u(x,t) dt = 0
*/
def ost_integration_ideal(Id,V,DV,W,LB,UB)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( type(OL=getopt(ol)) == -1 ) OL = 0;
if ( (VB = vars(append(LB,UB))) != [] )
if ( Param == 0 )
Param = VB;
else
Param = append(Param,VB);
TMP = separate_var(V,DV,W);
VX = TMP[0]; DVX = TMP[1];
VY = TMP[2]; DVY = TMP[3]; VYDVY = append(VY,DVY);
M = length(VY);
WW = newvect(M);
TI = [];
DT = 1; OLL = newvect(length(VY)); /* for option ol */
while ( Id != [] ) {
T = 1;
MIN = 0; J = 0; /* for #else */
P = car(Id);
if ( OL ) {
OLL = car(OL); OL = cdr(OL);
DT = car(OLL); OLL = cdr(OLL);
}
for ( I = 0; I < M; I++ ) {
WW[I] = W[I];
MWW = -WW;
WT = append(vtol(MWW),vtol(WW));
O = ord_w(P,VYDVY,WT)-OLL[I];
WW[I] = 0;
#if 0
if ( O > 0 ) {
if ( LB[I] != "-inf" && subst(DT,VY[J],LB[J]) )
T *= (VY[I]-LB[I])^O;
if ( UB[I] != "inf" && UB[I] != "+inf" && subst(DT,VY[J],UB[J]) )
T *= (UB[I]-VY[I])^O;
}
}
#else
if ( O > 0 ) {
if ( MIN ) {
if ( O < MIN ) {
MIN = O;
J = I;
}
} else {
MIN = O;
J = I;
}
}
}
if ( LB[J] != "-inf" && subst(DT,VY[J],LB[J]) )
T *= (VY[J]-LB[J])^MIN;
if ( UB[J] != "inf" && UB[J] != "+inf" && subst(DT,VY[J],UB[J]) )
T *= (UB[J]-VY[J])^MIN;
#endif
TI = cons(T*P,TI);
Id = cdr(Id);
}
return integration_ideal(TI,V,DV,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
}
def test_ost_integration_ideal()
{
B = [dt+(3*t^2-1)*x,dx+t^3-t];
return ost_integration_ideal(B,[t,x],[dt,dx],[1,0],[0],["inf"]);
}
def ost_integration_exp(C,E,V,DV,W,LB,UB)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
Id = []; OL = [];
TMP = separate_var(V,DV,W); VY = TMP[2]; DVY = TMP[3];
for ( I = 0; I < length(V); I++ ) {
T = red(diff(C,V[I])/C);
NT = nm(T); DT = dn(T);
Id = cons(DT*DV[I]-NT-DT*diff(E,V[I]),Id);
for ( OLL = [], J = 0; J < length(VY); J++ )
OLL = cons(ord_w(DT,[VY[J],DVY[J]],[-1,1]),OLL);
OL = cons(cons(DT,reverse(OLL)),OL);
}
G = ost_integration_ideal(Id,V,DV,W,LB,UB|param=Param,s0=S0,inhomo=IH,ht=HTrace,ol=OL);
return G;
}
def ost_integration_ideal2(Id,V,DV,W,LB,UB)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
if ( !NK_RESTRICTION_NOLOG ) T = time();
L = ann_heaviside(V,DV,LB,UB);
TI = ann_mul(L,Id,V,DV); /* ann_mul.rr が必要 */
if ( !NK_RESTRICTION_NOLOG ) {
print("-- ann_mul :",2);
time_print(T);
T = time();
}
INT = integration_ideal(TI,V,DV,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
if ( !NK_RESTRICTION_NOLOG ) {
print("-- integration_ideal :",2);
time_print(T);
}
return INT;
}
def ann_heaviside(V,DV,LB,UB)
{
L = [];
for ( I = 0; I < length(LB); I++ ) {
T = 1;
if ( type(LB[I]) > 1 ) {
if ( LB[I] != "-inf" )
error("ann_heaviside: invalid lower bound.");
} else
T *= V[I]-LB[I];
if ( type(UB[I]) > 1 ) {
if ( UB[I] != "inf" && UB[I] != "+inf" )
error("ann_heaviside: invalid upper bound.");
} else
T *= V[I]-UB[I];
L = cons(T*DV[I],L);
}
for ( ; I < length(V); I++ )
L = cons(DV[I],L);
return L;
}
def ost_integration_exp2(C,E,V,DV,W,LB,UB)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
Id = [];
for ( I = 0; I < length(V); I++ ) {
T = red(diff(C,V[I])/C);
NT = nm(T); DT = dn(T);
Id = cons(DT*DV[I]-NT-DT*diff(E,V[I]),Id);
}
G = ost_integration_ideal2(Id,V,DV,W,LB,UB|param=Param,s0=S0,inhomo=IH,ht=HTrace);
return G;
}
/**** added functions for difference operators ****/
/* Mellin transform */
/* es <=> x, s+1 <=> -x*dx */
/* mellin_trans((s+1)*es,[s],[es],[x],[dx]) = -dx*x^2-x */
/* 差分作用素は常に右側にあるものとする */
def mellin_trans(F,E,EV,V,DV)
{
if ( type(Form=getopt(form)) == -1 ) Form = 0;
if ( Form )
for ( I = 0; I < length(E); I++ )
F = subst(F,E[I],E[I]+1);
F = normal(F,E,EV|right=1);
for ( I = 0; I < length(E); I ++ )
F = mellin1(F,E[I],EV[I],V[I],DV[I]);
return F;
}
def mellin_trans2(F,E,EV,V,DV)
{
return mellin_trans(F,E,EV,V,DV|form=1);
}
def mellin1(F,E,EV,V,DV)
{
R = dsubst(F,E,-V*DV-1,[V,DV]);
RR = subst(R,EV,V);
return RR;
}
def dsubst(P,X,L,XL)
{
R = 0;
LL = dp_ptod(1,XL);
DL = dp_ptod(L,XL);
for ( I = 0; I <= deg(P,X); I++ ){
C = coef(P,I,X);
SC = C*dp_dtop(LL,XL);
R += SC;
LL = dp_weyl_mul(LL,DL);
}
return R;
}
/* output: [nm,dn] dn*nm (dn は nm の左からかかる) */
def inv_mellin_trans(F,V,DV,E,EV)
{
if ( type(Form=getopt(form)) == -1 ) Form = 0;
VDV = append(V,DV);
DF = dp_ptod(F,VDV);
N = length(V);
R = 0;
while ( DF ) {
H = dp_etov(DF);
for ( C = 1, I = 0; I < N; I++ ) {
S = H[I] - H[I+N];
C *= EV[I]^S;
C *= theta(E[I],H[I+N]);
}
R += dp_hc(DF)*C;
DF = dp_rest(DF);
}
R = red(R); NM = nm(R); DN = dn(R);
if ( Form )
for ( I = 0; I < length(E); I++ )
NM = subst(NM,E[I],E[I]-1);
return [normal(NM,E,EV),dn(R)];
}
def inv_mellin_trans2(F,V,DV,E,EV)
{
return inv_mellin_trans(F,V,DV,E,EV|form=1);
}
def theta(S,N)
{
for ( R=1, I = 0; I < N; I++ )
R *= -S-1-I;
return R;
}
/* right=0: 右にある変数を差分作用素の"左"に */
/* right=1: 左にある変数を差分作用素の"右"に */
def normal(P,E,EV)
{
if ( type(Right=getopt(right)) == -1 ) Right = 0;
DP = dp_ptod(P,EV);
N = length(EV);
R = 0;
while ( DP ) {
H = dp_etov(DP);
C = dp_hc(DP);
for ( I = 0; I < N; I++ )
if ( Right )
C = subst(C,E[I],E[I]-H[I]);
else
C = subst(C,E[I],E[I]+H[I]);
R += C*dp_vtoe(H);
DP = dp_rest(DP);
}
return dp_dtop(R,EV);
}
def integration_ideal_shift(Id,V,DV,W1,E,EV,W2)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
MV = pppppppp;
for ( P = DP = [], I = 0; I < length(E); I++ ) {
P = cons(strtov(rtostr(MV)+rtostr(I+1)),P);
DP = cons(strtov("d"+rtostr(MV)+rtostr(I+1)),DP);
}
P = reverse(P); DP = reverse(DP);
MId = map(mellin_trans,Id,E,EV,P,DP);
TMP1 = separate_var(V,DV,W1|weight=1);
TMP2 = separate_var(P,DP,W2|weight=1);
NV = append(append(TMP1[2],TMP2[2]),append(TMP1[0],TMP2[0]));
NDV = append(append(TMP1[3],TMP2[3]),append(TMP1[1],TMP2[1]));
NW = append(append(TMP1[6],TMP2[6]),append(TMP1[5],TMP2[5]));
J = integration_ideal(MId,NV,NDV,NW|param=Param,s0=S0,inhomo=IH,ht=HTrace);
if ( IH ) {
IMJ = map(inv_mellin_trans,J[0],P,DP,E,EV);
return [map(car_func,IMJ),inhomo_inv_mellin(J[1],P,DP,E,EV,IMJ)];
} else
return map(car_func,map(inv_mellin_trans,J,P,DP,E,EV));
}
def inhomo_inv_mellin(IH,P,DP,E,EV,DN)
{
if ( type(Sum=getopt(sum)) == -1 ) Sum = 0;
if ( type(Form=getopt(form)) == -1 ) Form = 0;
R = [];
while ( IH != [] ) {
H = car(IH);
R = cons([map(inhomo_inv_mellin_sub,H[0],P,DP,E,EV,Sum,Form),H[1]/car(DN)[1]],R);
IH = cdr(IH);
DN = cdr(DN);
}
return reverse(R);
}
def inhomo_inv_mellin_sub(IH,P,DP,E,EV,Sum,Form)
{
F = IH[0]; L = IH[1];
if ( Sum ) {
F = car(inv_mellin_trans(F,P,DP,E,EV|form=Form))-1;
for ( I = 0; I < length(Sum) && Sum[I] > 0; I++ )
L = subst(L,P[I],P[I]-1);
}
return [F,inv_mellin_trans(L,P,DP,E,EV|form=Form)];
}
def car_func(L)
{
return car(L);
}
/* グレブナー道場 例題7.4.29 */
def test_integration_ideal_shift()
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
F = t*(1-t);
A = ann(F);
I = cons(es-F,A);
V = [t]; DV = [dt];
E = [s]; EV = [es];
W1 = [1];
W2 = [0];
return integration_ideal_shift(I,V,DV,W1,E,EV,W2|inhomo=IH);
}
/* グレブナー道場 問題7.4.29 */
def test_integration_ideal_shift2()
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
F = t1*t2*(1-t1-t2);
A = ann(F);
I = cons(es-F,A);
V = [t1,t2]; DV = [dt1,dt2];
E = [s]; EV = [es];
W1 = [1,1];
W2 = [0];
return integration_ideal_shift(I,V,DV,W1,E,EV,W2|inhomo=IH);
}
/* グレブナー道場 例題7.4.30 */
def test_integration_ideal_shift3()
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
/* A = ndbf.ann_n([t,1-t,1+2*t]); */
A = [2*dt*t^3+(-dt-2*s0-2*s1-2*s2)*t^2+(-dt+s0-s1+2*s2)*t+s0];
AN = map(ptozp,map(subst,A,s0,1/2,s1,1/2+n,s2,1/2));
I = cons(en-(1-t),AN);
J = integration_ideal_shift(I,[t],[dt],[1],[n],[en],[0]|inhomo=IH);
return J;
}
def test_integration_ideal_shift4()
{
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
/* A = ndbf.ann_n([t,1-t,1+2*t]); */
A = [2*dt*t^3+(-dt-2*s0-2*s1-2*s2)*t^2+(-dt+s0-s1+2*s2)*t+s0];
AN = map(ptozp,map(subst,A,s0,1/2+m,s1,1/2+n,s2,1/2));
I = append([en-(1-t),em-(t)],AN);
J = integration_ideal_shift(I,[t],[dt],[1],[n,m],[en,em],[0,0]|inhomo=IH);
return J;
}
def ost_sum(Id,E,EV,W)
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
MV = pppppppp;
for ( P = DP = [], I = 0; I < length(E); I++ ) {
P = cons(strtov(rtostr(MV)+rtostr(I+1)),P);
DP = cons(strtov("d"+rtostr(MV)+rtostr(I+1)),DP);
}
P = reverse(P); DP = reverse(DP);
MId = map(mellin_trans2,Id,E,EV,P,DP);
for ( I = 0; I < length(W) && W[I] > 0; I++ )
MId = subst(MId,P[I],P[I]+1);
J = restriction_ideal(MId,P,DP,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
if ( IH ) {
IMJ = map(inv_mellin_trans2,J[0],P,DP,E,EV);
return [map(car_func,IMJ),inhomo_inv_mellin(J[1],P,DP,E,EV,IMJ|sum=W,form=1)];
} else
return map(car_func,map(inv_mellin_trans2,J,P,DP,E,EV));
}
def test_ost_sum()
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
IN = [(m1+1-n1)*em1-(m1+1),(m2+1-n+n1)*em2-(m2+1),(n1+1)*(m2-n+n1+1)*en1-(m1-n1)*(n-n1)];
I = subst(IN,n,10);
E = [n1,m1,m2];
EV = [en1,em1,em2];
W = [1,0,0];
R = ost_sum(I,E,EV,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
return R;
}
/* ost ex. 6.2 */
def test_ost_sum2()
{
if ( type(Param=getopt(param)) == -1 ) Param = 0;
if ( type(S0=getopt(s0)) == -1 ) S0 = -1;
if ( type(IH=getopt(inhomo)) == -1 ) IH = 0;
if ( type(HTrace=getopt(ht)) == -1 ) HTrace = 2;
I = [(n-k+1)*en-(n+1),(k+1)*ek-(n-k)];
E = [k,n];
EV = [ek,en];
W = [1,0];
R = ost_sum(I,E,EV,W|param=Param,s0=S0,inhomo=IH,ht=HTrace);
return R;
}
endmodule$
end$
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to nk_restriction.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