File: [local] / OpenXM / src / asir-contrib / packages / src / ok_diff.rr (download)
Revision 1.2, Mon Apr 11 12:04:37 2005 UTC (19 years, 1 month ago) by takayama
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_3_1_13b, RELEASE_1_2_3_12, KNOPPIX_2006, HEAD, DEB_REL_1_2_3-9
Changes since 1.1: +8 -8
lines
Fixed a bug of ok_diff.rr ; matrix_ ==> omatrix_
Added a new document directory ok_diff (a library for differential operators
by Okutani.)
|
/* $OpenXM: OpenXM/src/asir-contrib/packages/src/ok_diff.rr,v 1.2 2005/04/11 12:04:37 takayama Exp $
Old: Diff
*/
/*------------------------------------*/
/* Package for differential equations */
/*------------------------------------*/
#define _DEBUG_Odiff_ 0
#define _OPTIMIZE_Odiff_ 1
Odiff_ActList1 = [ newvect(2, [ [(1)/(x),[x]], newvect(1,[(1)/(x)]) ]) ]$
/* [[[F,V] DF],...] */
/* DF: vector with vector elements */
/* [ [F Dv1(F) ... Dv1^N(F)]
[Dv0(F) Dv1Dv0(F) ... Dv1^(N-1)Dv0(F)]
...
[Dv0^N(F)] ] */
Odiff_ActList2 = [ newvect(2, [ [(1)/(x*y),[x,y]],
newvect(1,[ newvect(1,[(1)/(x*y)]) ]) ]) ]$
/*---------*/
/* Utility */
/*---------*/
def odiff_is_int(N) {
if (type(N) <= 1) {
if (dn(N) == 1) {
return 1;
}
}
return 0;
}
def odiff_is_neg_int(N) {
if (odiff_is_int(N) && N < 0) {
return 1;
}
return 0;
}
def odiff_is_neg0_int(N) {
if (odiff_is_int(N) && N <= 0) {
return 1;
}
return 0;
}
/* [x1,x2,...] -> [dx1,dx2,...] */
def odiff_make_dv(V) {
if (type(V) != 4) {
error("parameter type error.");
}
Dv = [];
for (I = 0; I < length(V); I++) {
Dv = append(Dv, [strtov("d" + rtostr(V[I]))]);
}
return Dv;
}
def odiff_op_toasir(L_list, V) {
Dv = odiff_make_dv(V);
L_list_asir = [];
for (I = 0; I < length(L_list); I++) {
L = L_list[I];
L_asir = 0;
for (J = 0; J < length(L); J++) {
L_asir += L[J][0]*dp_dtop(dp_vtoe(omatrix_ltov(L[J][1])), Dv);
}
L_list_asir = append(L_list_asir, [L_asir]);
}
return L_list_asir;
}
def odiff_op_fromasir(D_list, V) {
Dv = odiff_make_dv(V);
L_list = [];
for (I = 0; I < length(D_list); I++) {
D = dp_ptod(D_list[I], Dv);
L = [];
if (D == 0) {
L = append(L, [[0,[0,0]]]);
}
else {
while (D != 0) {
L = append(L, [[dp_hc(D),vtol(dp_etov(dp_ht(D)))]]);
D -= dp_hm(D);
}
}
L_list = append(L_list, [L]);
}
return L_list;
}
/* translate coefficients of type 3 to of type 2. */
def odiff_op_rat2p(L_list) {
NewL_list = [];
for (I = 0; I < length(L_list); I++) {
L = L_list[I];
Dn = [];
for (J = 0; J < length(L); J++) {
Dn = append(Dn, [dn(red(L[J][0]))]);
}
if (Dn != []) {
Lcm = Dn[0];
for (J = 1; J < length(Dn); J++) {
Lcm *= sdiv(Dn[J], gcd(Lcm, Dn[J]));
}
NewL = [];
for (J = 0; J < length(L); J++) {
NewL = append(NewL, [ [red(L[J][0]*Lcm), L[J][1]] ]);
}
NewL_list = append(NewL_list, [NewL]);
}
}
return NewL_list;
}
/* translate coefficients of type 3 to of type 2 with integer coef. */
def odiff_op_rat2zp(L_list) {
L_list = odiff_op_rat2p(L_list);
NewL_list = [];
for (I = 0; I < length(L_list); I++) {
L = L_list[I];
Dn = [];
for (J = 0; J < length(L); J++) {
Dn = append(Dn, [dn(sdiv(L[J][0], ptozp(L[J][0])))]);
}
if (Dn != []) {
Lcm = Dn[0];
for (J = 1; J < length(Dn); J++) {
Lcm *= idiv(Dn[J], igcd(Lcm, Dn[J]));
}
NewL = [];
for (J = 0; J < length(L); J++) {
NewL = append(NewL, [ [L[J][0]*Lcm, L[J][1]] ]);
}
NewL_list = append(NewL_list, [NewL]);
}
}
return NewL_list;
}
/* asir polynomial -> sm1 string */
/* note! sm1 understand only polynomial with integer coefficients. */
def odiff_poly_tosm1(Poly, V) {
if (type(Poly) > 2 || type(V) != 4) {
error("parameter type error.");
}
if (Poly == 0) {
return "0";
}
Dpoly = dp_ptod(Poly, V);
Str_poly = "";
do {
Hc = dp_hc(Dpoly);
Ht = dp_ht(Dpoly);
HtV = dp_etov(Ht);
#if _DEBUG_Odiff_
print("Dpoly: ", 0); print(Dpoly);
print("Hc: ", 0); print(Hc);
print("Ht: ", 0); print(Ht);
print("HtV: ", 0); print(HtV);
print("size(HtV): ", 0); print(size(HtV));
#endif
/* coefficient */
if (type(Hc) < 2) {
Str_poly += " + (" + rtostr(Hc) + ")";
}
else {
#if _DEBUG_Odiff_
print("Hc is polynomial.");
#endif
Str_poly += " + (" + odiff_poly_tosm1(Hc, vars(Hc)) + ")";
}
#if _DEBUG_Odiff_
print("coefficient OK.");
#endif
/* power product */
Str_mono = "";
for (J = 0; J < size(HtV)[0]; J++) {
if (HtV[J] != 0) {
Str_mono += " " + rtostr(V[J]^HtV[J]);
}
}
Str_poly += Str_mono;
#if _DEBUG_Odiff_
print("power product OK.");
#endif
Dpoly -= dp_hm(Dpoly);
} while (Dpoly != 0);
return Str_poly;
}
/* differential operators with polynomial coefficients -> sm1 string */
/* L_list: [ [[fa(V),[a1,...,an]],...], ... ], coefficient must be polynomial. */
def odiff_op_tosm1(L_list, V) {
if (type(L_list[0]) <= 2) {
L_list = odiff_op_fromasir(L_list, V);
}
L_list = odiff_op_rat2zp(L_list);
L_str_list = [];
for (I = 0; I < length(L_list); I++) {
L_str = "";
for (J = 0; J < length(L_list[I]); J++) {
Poly = red(L_list[I][J][0]);
if (Poly != 0) {
Str_Poly = "(" + odiff_poly_tosm1(Poly, V) + ")";
Str_D = "";
Index = L_list[I][J][1];
for (K = 0; K < length(Index); K++) {
if (Index[K] != 0) {
Str_D += " d" + rtostr(V[K]^Index[K]);
}
}
L_str += " + " + Str_Poly + Str_D;
}
}
#if _DEBUG_Odiff_
print("L_str: ", 0); print(L_str);
#endif
if (L_str == "") {
L_str_list = append(L_str_list, ["0"]);
}
else {
L_str_list = append(L_str_list, [L_str]);
}
}
return L_str_list;
}
def odiff_get_param(D_ideal, V) {
if (type(D_ideal[0]) <= 2) {
D_ideal = odiff_op_fromasir(D_ideal, V);
}
Vars = vars(D_ideal);
for (Param = [], I = 0; I < length(Vars); I++) {
if (!odiff_lookup(Vars[I],V)) {
Param = append(Param, [Vars[I]]);
}
}
#if _DEBUG_Odiff_
print("generic parameter: ", 0); print(Param);
#endif
return Param;
}
/*-------------------------------*/
/* In the case of many variables */
/*-------------------------------*/
def odiff_act1(L, F, V) {
if (length(V) != 1) {
error("a length of list V is wrong.");
}
extern Odiff_ActList1;
/* search list */
for (I = 0; I < length(Odiff_ActList1); I++) {
if (Odiff_ActList1[I][0] == [F, V]) {
break;
}
}
N = I;
if (N == length(Odiff_ActList1)) {
/* make Odiff_ActList1[N] */
Odiff_ActList1 = append(Odiff_ActList1, [ newvect(2, [[F,V],newvect(1,[F])]) ]);
}
Table = Odiff_ActList1[N];
/* calculate order of L */
MaxOrd = 0;
for (I = 0; I < length(L); I++) {
if (MaxOrd < L[I][1][0]) {
MaxOrd = L[I][1][0];
}
}
/* if MaxOrd >= length(Odiff_ActList1[N][1][0]), extend Odiff_ActList1[N][1]. */
if (MaxOrd >= size(Table[1])[0]) {
New = newvect(MaxOrd+1, vtol(Table[1]));
for (I = size(Table[1])[0]; I <= MaxOrd; I++) {
New[I] = red(diff(New[I-1], V[0]));
}
Table[1] = New;
}
/* action L */
#if _DEBUG_Odiff_
print("calculate Lf... ", 0);
#endif
Lf = 0;
for (I = length(L) - 1; I >= 0; I--) {
Df = Table[1][L[I][1][0]];
Lf += L[I][0]*Df;
Lf = red(Lf);
}
#if _DEBUG_Odiff_
print("done.");
#endif
return Lf;
}
def odiff_act2(L, F, V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
extern Odiff_ActList2;
/* search list */
for (I = 0; I < length(Odiff_ActList2); I++) {
if (Odiff_ActList2[I][0] == [F, V]) {
break;
}
}
N = I;
if (N == length(Odiff_ActList2)) {
/* make Odiff_ActList2[N] */
Odiff_ActList2 = append(Odiff_ActList2,
[ newvect(2, [[F, V],newvect(1,[newvect(1,[F])])]) ]);
}
Table = Odiff_ActList2[N]; /* a vetor of size 2 */
/* calculate order of L */
MaxOrd = 0;
for (I = 0; I < length(L); I++) {
if (MaxOrd < L[I][1][0] + L[I][1][1]) {
MaxOrd = L[I][1][0] + L[I][1][1];
}
}
/* if MaxOrd >= length(Odiff_ActList2[N][1][0]), extend Odiff_ActList2[N][1]. */
if (MaxOrd >= size(Table[1])[0]) {
NewTable = newvect(MaxOrd+1);
for (I = 0; I <= MaxOrd; I++) {
if (I < size(Table[1])[0]) {
NewTable[I] = newvect(MaxOrd + 1 - I, vtol(Table[1][I]));
for (J = size(Table[1][I])[0]; J <= MaxOrd - I; J++) {
NewTable[I][J] = red(diff(NewTable[I][J-1], V[1]));
}
}
else {
NewTable[I] = newvect(MaxOrd + 1 - I, [red(diff(NewTable[I-1][0], V[0]))]);
for (J = 1; J <= MaxOrd - I; J++) {
NewTable[I][J] = red(diff(NewTable[I][J-1], V[1]));
}
}
}
Table[1] = NewTable;
}
/* action L */
Lf = 0;
for (I = length(L) - 1; I >= 0; I--) {
Df = Table[1][L[I][1][0]][L[I][1][1]];
Lf += L[I][0]*Df;
Lf = red(Lf);
}
return Lf;
}
/* D: list, [a1,...,an] */
/* F: rational expr */
/* V: list, [x1,...,xn] */
def odiff_act_mono(D, F, V) {
Var = [];
for (I = 0; I < length(V); I++) {
for (J = 0; J < D[I]; J++) {
Var = append(Var, [V[I]]);
}
}
return diff(F, Var);
}
def odiff_contain_uc(Obj) {
List = vars(Obj);
for (I = 0; I < length(List); I++) {
if (vtype(List[I]) == 1) {
return 1;
}
}
return 0;
}
/* L: [ [fa(V),[a1,...,an]],... ] */
/* F: rational expr */
/* V: list, [x1,...,xn] */
/* the return value of this function is reduced by using 'red()'. */
def odiff_act(L, F, V) {
if (type(L) <= 2) {
L = odiff_op_fromasir([L], V)[0];
}
#if _OPTIMIZE_Odiff_
if (length(V) == 1 && type(F) == 3) {
return odiff_act1(L, F, V);
}
/* if (length(V) == 2 && type(F) == 3 && !odiff_contain_uc(F)) { */
if (length(V) == 2 && type(F) == 3) {
return odiff_act2(L, F, V);
}
#endif
Lf = 0;
for (I = length(L) - 1; I >= 0; I--) {
Lf += L[I][0]*odiff_act_mono(L[I][1], F, V);
Lf = red(Lf);
}
return Lf;
}
def odiff_acts(L_list, F, V) {
return map(odiff_act, L_list, F, V);
}
/* Dim: a number of variables */
/* Deg: degree */
def odiff_make_basis(Dim, Deg) {
if(Dim <= 0 || Deg < 0) {
error("parameter range error.");
}
if (Dim == 1) {
Basis = [];
for (I = 0; I <= Deg; I++) {
Basis = append(Basis, [[[I]]]);
}
return Basis;
}
if (Dim > 1) {
Basis = [];
B = odiff_make_basis(Dim - 1,Deg);
for (I = 0; I <= Deg; I++) {
HBasis = [];
for (J = 0; J <= I; J++) {
for (K = 0; K < length(B[J]); K++) {
HBasis = append(HBasis, [append(B[J][K], [I - J])]);
}
}
Basis = append(Basis, [HBasis]);
}
return Basis;
}
}
def odiff_get_homog(Poly, Deg, V) {
Basis = odiff_make_basis(length(V), Deg);
HPoly = 0;
for (I = 0; I < length(Basis[Deg]); I++) {
Coef = Poly;
PP = 1;
for (J = 0; J < length(Basis[Deg][I]); J++) {
Coef = coef(Coef, Basis[Deg][I][J], V[J]);
PP *= V[J]^Basis[Deg][I][J];
}
HPoly += Coef*PP;
}
return HPoly;
}
/* Sol: list, [U, UC] */
/* Value: rat or list */
def odiff_assign_uc(Sol,Value) {
if (type(Value) >= 0 && type(Value) <= 3) {
U = Sol[0];
for (I = 0; I < length(Sol[1]); I++) {
U = subst(U, Sol[1][I], Value);
}
return U;
}
if (type(Value) == 4) {
U = Sol[0];
for (I = 0; I < length(Sol[1]); I++) {
U = subst(U, Sol[1][I], Value[I]);
}
return U;
}
error("parameter type error");
}
def odiff_rat_solve(L_list, Dn, N, V) {
/* OK: calculate the degree of U */
Deg_U = N;
#if _DEBUG_Odiff_
print("Deg_U: ", 0); print(Deg_U);
#endif
/* OK: make a basis of U */
Basis_U = odiff_make_basis(length(V), Deg_U);
#if _DEBUG_Odiff_
print("Basis_U: ", 0); print(Basis_U);
#endif
/* OK: make coefficients of U */
Coef_U = [];
for (I = 0; I <= Deg_U; I++) {
HCoef_U = [];
for (J = 0; J < length(Basis_U[I]); J++) {
HCoef_U = append(HCoef_U, [uc()]);
}
Coef_U = append(Coef_U, [HCoef_U]);
}
#if _DEBUG_Odiff_
print("Coef_U: ", 0); print(Coef_U);
#endif
/* OK: make power series U */
U = 0;
for (I = 0; I <= Deg_U; I++) {
for (J = 0; J < length(Basis_U[I]); J++) {
U += Coef_U[I][J]*dp_dtop(dp_vtoe(omatrix_ltov(Basis_U[I][J])), V);
}
}
#if _DEBUG_Odiff_
print("U: ", 0); print(U);
#endif
/* OK: calculate Lu_list */
List = odiff_acts(L_list, U/Dn, V);
Lu_list = [];
for (I = 0; I < length(List); I++) {
Lu_list = append(Lu_list, [nm(List[I])]);
}
#if _DEBUG_Odiff_
print("Lu_list: ", 0); print(Lu_list);
#endif
/* OK: new version */
Coef_Lu_list = [];
for (I = 0; I < length(Lu_list); I++) {
if ((Dp_Lu = dp_ptod(Lu_list[I], V)) != 0) {
while (Dp_Lu != 0) {
Coef_Lu_list = append(Coef_Lu_list, [dp_hc(Dp_Lu)]);
Dp_Lu -= dp_hm(Dp_Lu);
}
}
}
if (Coef_Lu_list == []) {
Coef_Lu_list = [0];
}
#if _DEBUG_Odiff_
print("Coef_Lu_list: ", 0); print(Coef_Lu_list);
#endif
/* OK: calculate the size of matrix A and vector X */
Size_Row = length(Coef_Lu_list);
Size_Col = fac(length(V) + Deg_U)/(fac(length(V))*fac(Deg_U));
#if _DEBUG_Odiff_
print("Size of A: [", 0);
print(Size_Row, 0); print(",", 0); print(Size_Col, 0);
print("]");
#endif
/* OK: make a vector X */
X = newvect(Size_Col);
K = 0;
for (I = 0; I <= Deg_U; I++) {
for (J = 0; J < length(Coef_U[I]); J++) {
X[K] = Coef_U[I][J];
K++;
}
}
#if _DEBUG_Odiff_
print("X: ", 0); print(X);
#endif
/* OK: make a matrix A */
A = newmat(Size_Row, Size_Col);
for (Row = 0; Row < Size_Row; Row++) {
Coef_Nm = nm(Coef_Lu_list[Row]);
Coef_Dn = dn(Coef_Lu_list[Row]);
for (Col = 0; Col < Size_Col; Col++) {
A[Row][Col] = coef(Coef_Nm, 1, X[Col])/Coef_Dn;
}
}
#if _DEBUG_Odiff_
print("A: ", 0); print(A);
#endif
/* OK: solve A*X = 0 */
Sol = omatrix_solve(A, X, newvect(Size_Row));
#if _DEBUG_Odiff_
print("Sol: ", 0); print(Sol);
#endif
/* OK: assign Sol to U */
Size_Basis = fac(length(V)+N)/(fac(length(V))*fac(N));
for (I = 0; I < Size_Col; I++) {
if (I < Size_Basis) {
U = subst(U, Sol[I][0], Sol[I][1]);
}
else {
U = subst(U, Sol[I][0], 0);
}
}
#if _DEBUG_Odiff_
print("U: ", 0); print(U);
#endif
/* unknown coefficient */
Vars_U = vars(U);
UC = [];
for (I = 0; I < length(Vars_U); I++) {
if (vtype(Vars_U[I]) == 1) {
UC = append(UC, [Vars_U[I]]);
}
}
#if _DEBUG_Odiff_
print("UC: ", 0); print(UC);
#endif
return [red(U/Dn), UC];
}
def odiff_poly_solve(L_list, N, V) {
return odiff_rat_solve(L_list, 1, N, V);
}
def odiff_pochhammer(A, K) {
A_K = 1;
for (I = 0; I < K; I++) {
A_K *= (A+I);
}
return A_K;
}
/*--------------------------*/
/* the case of one variable */
/*--------------------------*/
def odiff_op_hg1(A,B,C,V) {
if (length(V) != 1) {
error("a length of list V is too long.");
}
return [ [ [V[0]*(1-V[0]),[2]],
[C-(A+B+1)*V[0],[1]],
[-A*B,[0]] ] ];
}
def odiff_act_hg1(A,B,C,F,V) {
return odiff_acts(odiff_op_hg1(A,B,C,V),F,V);
}
def odiff_sol_hg1_type1(Deg,A,B,C,I,V) {
for (Sol = 0, M = 0; M <= Deg; M++) {
A_m = odiff_pochhammer(A, M);
B_m = odiff_pochhammer(B, M);
C_m = odiff_pochhammer(C, M);
I_m = odiff_pochhammer(I, M);
Sol += (A_m*B_m)/(C_m*I_m)*V[0]^M;
}
return Sol;
}
def odiff_poly_solve_hg1(A,B,C,V) {
A_is_neg0_int = odiff_is_neg0_int(A);
B_is_neg0_int = odiff_is_neg0_int(B);
/* polynomial solution does not exist. */
if (!(A_is_neg0_int || B_is_neg0_int)) {
return [0,[]];
}
/* polynomial solution exists. */
if (!odiff_is_neg0_int(C)) {
if (A_is_neg0_int && B_is_neg0_int) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_hg1_type1(Max,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (A_is_neg0_int && !B_is_neg0_int) {
Sol = odiff_sol_hg1_type1(-A,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (!A_is_neg0_int && B_is_neg0_int) {
Sol = odiff_sol_hg1_type1(-B,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
else {
if (!A_is_neg0_int && C <= B && B <= 0) {
Sol = odiff_sol_hg1_type1(-B,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (!A_is_neg0_int && B <= C-1) {
Sol = V[0]^(1-C)*odiff_sol_hg1_type1(-B+C-1,A-C+1,B-C+1,1,2-C,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C <= A && A <= 0 && !B_is_neg0_int) {
Sol = odiff_sol_hg1_type1(-A,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C <= A && A <= 0 && C <= B && B <= 0) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_hg1_type1(Max,A,B,C,1,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C <= A && A <= 0 && B <= C-1) {
Sol1 = odiff_sol_hg1_type1(-A,A,B,C,1,V);
Sol2 = V[0]^(1-C)*odiff_sol_hg1_type1(-B+C-1,A-C+1,B-C+1,1,2-C,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (A <= C-1 && !B_is_neg0_int) {
Sol = V[0]^(1-C)*odiff_sol_hg1_type1(-A+C-1,A-C+1,B-C+1,1,2-C,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (A <= C-1 && C <= B && B <= 0) {
Sol1 = odiff_sol_hg1_type1(-B,A,B,C,1,V);
Sol2 = V[0]^(1-C)*odiff_sol_hg1_type1(-A+C-1,A-C+1,B-C+1,1,2-C,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (A <= C-1 && B <= C-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[0]^(1-C)*odiff_sol_hg1_type1(Max+C-1,A-C+1,B-C+1,1,2-C,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
/* return odiff_poly_solve(odiff_op_hg1(A,B,C,V),N,V); */
}
def odiff_poly_solve_hg1_check(A,B,C,V) {
Sol = odiff_poly_solve_hg1(A,B,C,V);
if (odiff_act_hg1(A,B,C,Sol[0],V) != [0]) {
error("odiff_poly_solve_hg1 contains bugs.");
return 0;
}
return 1;
}
def odiff_op_hg1_rat(A,B,C,K1,K2,V) {
A3 = A+B-2*K1-2*K2+1;
A2 = -A-B-C+4*K1+2*K2-1;
A1 = C-2*K1;
B2 = (A-K1-K2)*(B-K1-K2);
B1 = -(A-K1)*(B-K1)+(K1+K2)*C-K1*(K1+2*K2+1);
B0 = K1*(K1-C+1);
return [ [ [V[0]^2*(1-V[0])^2,[2]],
[A3*V[0]^3+A2*V[0]^2+A1*V[0],[1]],
[B2*V[0]^2+B1*V[0]+B0,[0]] ] ];
}
def odiff_act_hg1_rat(A,B,C,K1,K2,F,V) {
return odiff_acts(odiff_op_hg1_rat(A,B,C,K1,K2,V),F,V);
}
def odiff_rat_solve_hg1(A,B,C,V) {
A_is_int = odiff_is_int(A);
B_is_int = odiff_is_int(B);
if (!(A_is_int || B_is_int)) {
return [0,[]];
}
Dim_poly = length(odiff_poly_solve_hg1(A,B,C,V)[1]);
if (Dim_poly == 2) {
return [0,[]];
}
Dim_rat = 0;
Sol1 = 0;
C_is_int = odiff_is_int(C);
if (C_is_int && C >= 2) {
if (A_is_int && 1 <= A && A <= C-1) {
if (B_is_int && 1 <= B && B <= C-1) {
Deg = -((A > B) ? (A):(B))+C-1;
Sol1 = odiff_sol_hg1_type1(Deg,A-C+1,B-C+1,1,2-C,V)/V[0]^(C-1);
Dim_rat++;
}
else {
Sol1 = odiff_sol_hg1_type1(-A+C-1,A-C+1,B-C+1,1,2-C,V)/V[0]^(C-1);
Dim_rat++;
}
}
else {
if (B_is_int && 1 <= B && B <= C-1) {
Sol1 = odiff_sol_hg1_type1(-B+C-1,A-C+1,B-C+1,1,2-C,V)/V[0]^(C-1);
Dim_rat++;
}
}
}
Sol2 = 0;
CC = A+B-C+1;
CC_is_int = odiff_is_int(CC);
if (CC_is_int && CC >= 2) {
if (A_is_int && 1 <= A && A <= CC-1) {
if (B_is_int && 1 <= B && B <= CC-1) {
Deg = -((A > B) ? (A):(B))+CC-1;
Sol2 = odiff_sol_hg1_type1(Deg,A-CC+1,B-CC+1,1,2-CC,V)/V[0]^(CC-1);
Sol2 = subst(Sol2,V[0],1-V[0]);
Dim_rat++;
}
else {
Sol2 = odiff_sol_hg1_type1(-A+CC-1,A-CC+1,B-CC+1,1,2-CC,V)/V[0]^(CC-1);
Sol2 = subst(Sol2,V[0],1-V[0]);
Dim_rat++;
}
}
else {
if (B_is_int && 1 <= B && B <= CC-1) {
Sol2 = odiff_sol_hg1_type1(-B+CC-1,A-CC+1,B-CC+1,1,2-CC,V)/V[0]^(CC-1);
Sol2 = subst(Sol2,V[0],1-V[0]);
Dim_rat++;
}
}
}
if (Sol1 != 0 && Sol2 != 0) {
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (Sol1 != 0 && Sol2 == 0) {
UC1 = uc();
return [UC1*Sol1,[UC1]];
}
if (Sol1 == 0 && Sol2 != 0) {
UC2 = uc();
return [UC2*Sol2,[UC2]];
}
if (Sol1 == 0 && Sol2 == 0) {
return [0,[]];
}
}
/*---------------------------*/
/* the case of two variables */
/*---------------------------*/
/* appell1 */
def odiff_op_appell1(A,B1,B2,C,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
return [ [ [V[0]*(1-V[0]),[2,0]],
[(1-V[0])*V[1],[1,1]],
[C-(A+B1+1)*V[0],[1,0]],
[-B1*V[1],[0,1]],
[-A*B1,[0,0]] ],
[ [V[1]*(1-V[1]),[0,2]],
[(1-V[1])*V[0],[1,1]],
[C-(A+B2+1)*V[1],[0,1]],
[-B2*V[0],[1,0]],
[-A*B2,[0,0]] ],
[ [V[0]-V[1],[1,1]],
[-B2,[1,0]],
[B1,[0,1]] ] ];
}
def odiff_act_appell1(A,B1,B2,C,F,V) {
return odiff_acts(odiff_op_appell1(A,B1,B2,C,V),F,V);
}
def odiff_poly_solve_appell1(A,B1,B2,C,N,V) {
return odiff_poly_solve(odiff_op_appell1(A,B1,B2,C,V),N,V);
}
def odiff_poly_solve_appell1_check(Sol,A,B1,B2,C,V) {
Lf_list = odiff_act_appell1(A,B1,B2,C,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
def odiff_poly_solve_appell1s(R_A,R_B1,R_B2,R_C,N,V) {
Sols = [];
for (C = R_C[0]; C <= R_C[1]; C += R_C[2]) {
for (A = R_A[0]; A <= R_A[1]; A += R_A[2]) {
for (B1 = R_B1[0]; B1 <= R_B1[1]; B1 += R_B1[2]) {
for (B2 = R_B2[0]; B2 <= R_B2[1]; B2 += R_B2[2]) {
Sols = append(Sols, [[A,B1,B2,C,odiff_poly_solve_appell1(A,B1,B2,C,N,V)]]);
}
}
}
}
return Sols;
}
/* appell2 */
def odiff_op_appell2(A,B1,B2,C1,C2,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
return [ [ [V[0]*(1-V[0]),[2,0]],
[-V[0]*V[1],[1,1]],
[C1-(A+B1+1)*V[0],[1,0]],
[-B1*V[1],[0,1]],
[-A*B1,[0,0]] ],
[ [V[1]*(1-V[1]),[0,2]],
[-V[0]*V[1],[1,1]],
[C2-(A+B2+1)*V[1],[0,1]],
[-B2*V[0],[1,0]],
[-A*B2,[0,0]] ],
[ [V[0]*V[1],[2,1]],
[-V[0]*V[1],[1,2]],
[B2*V[0],[2,0]],
[-B1*V[1],[0,2]],
[C1*V[1]-C2*V[0],[1,1]],
[B2*C1,[1,0]],
[-B1*C2,[0,1]] ] ];
}
def odiff_act_appell2(A,B1,B2,C1,C2,F,V) {
return odiff_acts(odiff_op_appell2(A,B1,B2,C1,C2,V),F,V);
}
def odiff_poly_solve_appell2(A,B1,B2,C1,C2,N,V) {
return odiff_poly_solve(odiff_op_appell2(A,B1,B2,C1,C2,V),N,V);
}
def odiff_poly_solve_appell2_check(Sol,A,B1,B2,C1,C2,V) {
Lf_list = odiff_act_appell2(A,B1,B2,C1,C2,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
/* appell3 */
def odiff_op_appell3(A1,A2,B1,B2,C,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
return [ [ [V[0]*(1-V[0]),[2,0]],
[V[1],[1,1]],
[C-(A1+B1+1)*V[0],[1,0]],
[-A1*B1,[0,0]] ],
[ [V[1]*(1-V[1]),[0,2]],
[V[0],[1,1]],
[C-(A2+B2+1)*V[1],[0,1]],
[-A2*B2,[0,0]] ] ];
}
def odiff_act_appell3(A1,A2,B1,B2,C,F,V) {
return odiff_acts(odiff_op_appell3(A1,A2,B1,B2,C,V),F,V);
}
def odiff_poly_solve_appell3(A1,A2,B1,B2,C,N,V) {
return odiff_poly_solve(odiff_op_appell3(A1,A2,B1,B2,C,V),N,V);
}
def odiff_poly_solve_appell3_check(Sol,A1,A2,B1,B2,C,V) {
Lf_list = odiff_act_appell3(A1,A2,B1,B2,C,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
/* appell4 */
def odiff_pseries_appell4(A,B,C1,C2,N,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
Basis = odiff_make_basis(2, N);
Pseries = 0;
for (I = 0; I < length(Basis); I++) {
for (J = 0; J < length(Basis[I]); J++) {
A_K = odiff_pochhammer(A, Basis[I][J][0] + Basis[I][J][1]);
B_K = odiff_pochhammer(B, Basis[I][J][0] + Basis[I][J][1]);
C1_M = odiff_pochhammer(C1, Basis[I][J][0]);
C2_N = odiff_pochhammer(C2, Basis[I][J][1]);
Coef = A_K*B_K/(C1_M*C2_N*fac(Basis[I][J][0])*fac(Basis[I][J][1]));
Pseries += red(Coef)*dp_dtop(dp_vtoe(omatrix_ltov(Basis[I][J])), V);
}
}
return Pseries;
}
def odiff_op_appell4(A,B,C1,C2,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
return [ [ [V[0]*(1-V[0]),[2,0]],
[-2*V[0]*V[1],[1,1]],
[-V[1]^2,[0,2]],
[(C1-(A+B+1)*V[0]),[1,0]],
[-(A+B+1)*V[1],[0,1]],
[-A*B,[0,0]] ],
[ [V[1]*(1-V[1]),[0,2]],
[-2*V[0]*V[1],[1,1]],
[-V[0]^2,[2,0]],
[(C2-(A+B+1)*V[1]),[0,1]],
[-(A+B+1)*V[0],[1,0]],
[-A*B,[0,0]] ] ];
}
def odiff_act_appell4(A,B,C1,C2,F,V) {
return odiff_acts(odiff_op_appell4(A,B,C1,C2,V),F,V);
}
def odiff_sol_appell4_type1(Deg,A,B,C1,C2,V) {
Basis = odiff_make_basis(2, Deg);
for (Sol = 0, I = 0; I < length(Basis); I++) {
for (J = 0; J < length(Basis[I]); J++) {
A_mn = odiff_pochhammer(A, Basis[I][J][0]+Basis[I][J][1]);
B_mn = odiff_pochhammer(B, Basis[I][J][0]+Basis[I][J][1]);
C1_m = odiff_pochhammer(C1, Basis[I][J][0]);
C2_n = odiff_pochhammer(C2, Basis[I][J][1]);
N1_m = odiff_pochhammer(1, Basis[I][J][0]);
N2_n = odiff_pochhammer(1, Basis[I][J][1]);
Sol += (A_mn*B_mn)/(C1_m*C2_n*N1_m*N2_n)*V[0]^Basis[I][J][0]*V[1]^Basis[I][J][1];
}
}
return Sol;
}
def odiff_sol_appell4_type2(Deg,A,B,C1,C2,N1,N2,V) {
Basis = odiff_make_basis(2, Deg);
for (Sol = 0, I = 0; I < length(Basis); I++) {
for (J = 0; J < length(Basis[I]); J++) {
A_mn = odiff_pochhammer(A, Basis[I][J][0]+Basis[I][J][1]);
B_mn = odiff_pochhammer(B, Basis[I][J][0]+Basis[I][J][1]);
C1_m = odiff_pochhammer(C1, Basis[I][J][0]);
C2_n = odiff_pochhammer(C2, Basis[I][J][1]);
N1_m = odiff_pochhammer(N1, Basis[I][J][0]);
N2_n = odiff_pochhammer(N2, Basis[I][J][1]);
Sol += (A_mn*B_mn)/(C1_m*C2_n*N1_m*N2_n)*V[0]^Basis[I][J][0]*V[1]^Basis[I][J][1];
}
}
return Sol;
}
def odiff_poly_solve_appell4(A,B,C1,C2,V) {
A_is_neg0_int = odiff_is_neg0_int(A);
B_is_neg0_int = odiff_is_neg0_int(B);
C1_is_neg0_int = odiff_is_neg0_int(C1);
C2_is_neg0_int = odiff_is_neg0_int(C2);
if (!(A_is_neg0_int || B_is_neg0_int)) {
return [0,[]];
}
if (!C1_is_neg0_int && !C2_is_neg0_int) {
if (A_is_neg0_int && !B_is_neg0_int) {
Sol = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (!A_is_neg0_int && B_is_neg0_int) {
Sol = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (A_is_neg0_int && B_is_neg0_int) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_appell4_type1(Max,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
if (C1_is_neg0_int && !C2_is_neg0_int) {
if (!A_is_neg0_int) {
if (C1 <= B && B <= 0) {
Sol = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (B <= C1-1) {
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
if (C1 <= A && A <= 0) {
if (!B_is_neg0_int) {
Sol = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1 <= B && B <= 0) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_appell4_type1(Max,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (B <= C1-1) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
}
if (A <= C1-1) {
if (!B_is_neg0_int) {
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1 <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C1-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(Max+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
}
if (!C1_is_neg0_int && C2_is_neg0_int) {
if (!A_is_neg0_int) {
if (C2 <= B && B <= 0) {
Sol = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (B <= C2-1) {
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
if (C2 <= A && A <= 0) {
if (!B_is_neg0_int) {
Sol = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C2 <= B && B <= 0) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_appell4_type1(Max,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (B <= C2-1) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
}
if (A <= C2-1) {
if (!B_is_neg0_int) {
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C2 <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C2-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(Max+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
}
if (C1_is_neg0_int && C2_is_neg0_int) {
Min_C = (C1 < C2) ? C1:C2;
Max_C = (C1 > C2) ? C1:C2;
if (!A_is_neg0_int) {
if (Max_C <= B && B <= 0) {
Sol = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1 > C2 && C2 <= B && B <= C1-1) {
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1 < C2 && C1 <= B && B <= C2-1) {
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C1+C2-2) {
Sol = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
if (Max_C <= A && A <= 0) {
if (!B_is_neg0_int) {
Sol = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (Max_C <= B && B <= 0) {
Max = (-A <= -B) ? (-A):(-B);
Sol = odiff_sol_appell4_type1(Max,A,B,C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1 > C2 && C2 <= B && B <= C1-1) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 < C2 && C1 <= B && B <= C2-1) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol3 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc(); UC3 = uc();
return [UC1*Sol1+UC2*Sol2+UC3*Sol3,[UC1,UC2,UC3]];
}
if (B <= C1+C2-2) {
Sol1 = odiff_sol_appell4_type1(-A,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
}
if (C1 > C2 && C2 <= A && A <= C1-1) {
if (!B_is_neg0_int) {
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (Max_C <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 > C2 && C2 <= B && B <= C1-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[0]^(1-C1)*odiff_sol_appell4_type1(Max+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C1+C2-2) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
}
if (C1 < C2 && C1 <= A && A <= C2-1) {
if (!B_is_neg0_int) {
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (Max_C <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 < C2 && C1 <= B && B <= C2-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[1]^(1-C2)*odiff_sol_appell4_type1(Max+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C1+C2-2) {
Sol1 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
}
if (C1+C2-1 <= A && A <= Min_C-1) {
if (!B_is_neg0_int) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (Max_C <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol3 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc(); UC3 = uc();
return [UC1*Sol1+UC2*Sol2+UC3*Sol3,[UC1,UC2,UC3]];
}
if (C1 > C2 && C2 <= B && B <= C1-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 < C2 && C1 <= B && B <= C2-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Max = (-A <= -B) ? (-A):(-B);
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(Max+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(Max+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (B <= C1+C2-2) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-A+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
Sol3 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc(); UC3 = uc();
return [UC1*Sol1+UC2*Sol2+UC3*Sol3,[UC1,UC2,UC3]];
}
}
if (A <= C1+C2-2) {
if (!B_is_neg0_int) {
Sol = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
if (Max_C <= B && B <= 0) {
Sol1 = odiff_sol_appell4_type1(-B,A,B,C1,C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 > C2 && C2 <= B && B <= C1-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1 < C2 && C1 <= B && B <= C2-1) {
Sol1 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
Sol2 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc();
return [UC1*Sol1+UC2*Sol2,[UC1,UC2]];
}
if (C1+C2-1 <= B && B <= Min_C-1) {
Sol1 = V[0]^(1-C1)*odiff_sol_appell4_type1(-B+C1-1,A-C1+1,B-C1+1,2-C1,C2,V);
Sol2 = V[1]^(1-C2)*odiff_sol_appell4_type1(-B+C2-1,A-C2+1,B-C2+1,C1,2-C2,V);
Sol3 = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(-A+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC1 = uc(); UC2 = uc(); UC3 = uc();
return [UC1*Sol1+UC2*Sol2+UC3*Sol3,[UC1,UC2,UC3]];
}
if (B <= C1+C2-2) {
Max = (-A <= -B) ? (-A):(-B);
Sol = V[0]^(1-C1)*V[1]^(1-C2)*odiff_sol_appell4_type1(Max+C1+C2-2,A-C1-C2+2,B-C1-C2+2,2-C1,2-C2,V);
UC = uc();
return [UC*Sol,[UC]];
}
}
}
error("odiff_poly_solve_appell4 contains bugs. parameter is "+rtostr([A,B,C1,C2])+".");
}
def odiff_poly_solve_appell4_check(A,B,C1,C2,V) {
Sol = odiff_poly_solve_appell4(A,B,C1,C2,V);
Lf_list = odiff_act_appell4(A,B,C1,C2,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
error("odiff_poly_solve_appell4 contains bugs."
+"the parameter is "+rtostr([A,B,C1,C2])+".");
return 0;
}
}
return 1;
}
def odiff_op_appell4_rat(A,B,C1,C2,K1,K2,K3,V) {
C1Gxx = V[0]^2*(1-V[0])*((V[0]-1)^2-2*(V[0]+1)*V[1]+V[1]^2);
C1Gxy = -2*V[0]^2*V[1]*((V[0]-1)^2-2*(V[0]+1)*V[1]+V[1]^2);
C1Gyy = -V[0]*V[1]^2*((V[0]-1)^2-2*(V[0]+1)*V[1]+V[1]^2);
C1Gx = V[0]*( (x-1)^2*(C1-2*K1-(1+A+B)*V[0]
+2*(K1+K2+2*K3)*V[0])
-2*(-2*K1+2*(K1+K2+2*K3)*V[0]^2+C1*(V[0]+1)
-V[0]*(1+A+B-2*K2+(1+A+B)*V[0]))*V[1]
+(C1-2*K1-(1+A+B)*V[0]+2*(K1+K2+2*K3)*V[0])*V[1]^2 );
C1Gy = -V[0]*V[1]*( (V[0]-1)*(-1+2*K1+2*K2+A*(V[0]-1)+B*(V[0]-1)
+V[0]-2*(K1+K2+2*K3)*V[0])
-2*(1+A+B-2*K1-2*K2-2*K3+V[0]
+(A+B-2*(K1+K2+2*K3))*V[0])*V[1]
+(1+A+B-2*K1-2*K2-4*K3)*V[1]^2 );
C1G = ( (V[0]-1)*(-K1^2*(V[0]-1)^2
+K1*(V[0]-1)*(1-C1+(A+B-2*K2-4*K3)*V[0])
+V[0]*((A-K2)*(B-K2)-2*C1*K3
+(A-K2-2*K3)*(-B+K2+2*K3)*V[0]))
+2*(-K1*(1-C1+K1)+A*B*V[0]
-(K1*(1+A+B-C1-2*K2)+(A+B-K2)*K2+(1+A+B-C1-2*K2)*K3)*V[0]
+(-A+K1+K2+2*K3)*(-B+K1+K2+2*K3)*V[0]^2)*V[1]
-(-K1*(1-C1+K1)+(-A+K1+K2+2*K3)*(-B+K1+K2+2*K3)*V[0])*V[1]^2 );
C2Gxx = -V[0]^2*V[1]*((-1+V[0])^2-2*(1+V[0])*V[1]+V[1]^2);
C2Gxy = -2*V[0]*V[1]^2*((-1+V[0])^2-2*(1+V[0])*V[1]+V[1]^2);
C2Gx = -V[1]*V[0]*( (-1+V[0])*(-1-A-B+2*K1+2*K2+V[0]
+(A+B-2*(K1+K2+2*K3))*V[0])
-2*(1+A+B-2*K1-2*K2-2*K3+V[0]
+(A+B-2*(K1+K2+2*K3))*V[0])*V[1]
+(1+A+B-2*K1-2*K2-4*K3)*V[1]^2 );
C2Gyy = -V[1]*(-1+V[1])*V[1]*((-1+V[0])^2-2*(1+V[0])*V[1]+V[1]^2);
C2Gy = V[1]*( -2*K2*(-1+V[0])^2
-(1+A+B-2*K1)*V[1]
+(-(1+A+B-2*K1)*(-2+V[0])*V[0]
+4*K3*(1+V[0]^2)
+2*K2*(3+V[0]^2))*V[1]
+2*(1+A+B-2*K1-3*K2-4*K3+V[0]
+(A+B-2*(K1+K2+2*K3))*V[0])*V[1]^2
-(1+A+B-2*K1-2*K2-4*K3)*V[1]^3
+C2*((-1+V[0])^2-2*(1+V[0])*V[1]+V[1]^2) );
C2G = ( K2*(1-C2+K2)*(-1+V[0])^2
-(K1^2+2*K2+2*K1*K2+3*K2^2+4*K2*K3-2*C2*(K2+K3)
-2*((-1+C2)*(K2+K3)+K1*(K1+2*(K2+K3)))*V[0]
+(K1+K2+2*K3)^2*V[0]^2
+A*(-1+V[0])*(K1+K2+B*(-1+V[0])-(K1+K2+2*K3)*V[0])
-B*(-1+V[0])*(-K1-K2+(K1+K2+2*K3)*V[0]))*V[1]
+(2*A*B-2*A*K1-2*B*K1+2*K1^2+K2-2*A*K2-2*B*K2-C2*K2+4*K1*K2
+3*K2^2-2*A*K3-2*B*K3-2*C2*K3+4*K1*K3+8*K2*K3+4*K3^2
+2*(-A+K1+K2+2*K3)*(-B+K1+K2+2*K3)*V[0])*V[1]^2
-(-A+K1+K2+2*K3)*(-B+K1+K2+2*K3)*V[1]^3 );
return [ [ [C1Gxx,[2,0]],[C1Gxy,[1,1]],[C1Gyy,[0,2]],
[C1Gx,[1,0]],[C1Gy,[0,1]],[C1G,[0,0]] ],
[ [C2Gyy,[0,2]],[C2Gxy,[1,1]],[C2Gxx,[2,0]],
[C2Gx,[1,0]],[C2Gy,[0,1]],[C2G,[0,0]] ] ];
}
def odiff_act_appell4_rat(A,B,C1,C2,K1,K2,K3,Nm,V) {
return odiff_acts(odiff_op_appell4_rat(A,B,C1,C2,K1,K2,K3,V),Nm,V);
}
def odiff_rat_solve_appell4(A,B,C1,C2,K1,K2,K3,N,V) {
Dn = V[0]^K1*V[1]^K2*((V[0]-V[1])^2-2*(V[0]+V[1])+1)^K3;
Sol = odiff_poly_solve(odiff_op_appell4_rat(A,B,C1,C2,K1,K2,K3,V),N,V);
return [red(Sol[0]/Dn), Sol[1]];
}
/* selberg2 */
def odiff_op_selberg2(A,B,C,S,V) {
if (length(V) != 2) {
error("a length of list V is wrong.");
}
return [ [ [V[0]*(1-V[0])*(V[0]-V[1]),[2,0]],
[(C-1/S-(A+B+1-1/S)*V[0])*(V[0]-V[1])+1/S*V[0]*(1-V[0]),[1,0]],
[-1/S*V[1]*(1-V[1]),[0,1]],
[-A*B*(V[0]-V[1]),[0,0]] ],
[ [V[1]*(1-V[1])*(V[1]-V[0]),[0,2]],
[(C-1/S-(A+B+1-1/S)*V[1])*(V[1]-V[0])+1/S*V[1]*(1-V[1]),[0,1]],
[-1/S*V[0]*(1-V[0]),[1,0]],
[-A*B*(V[1]-V[0]),[0,0]] ] ];
}
def odiff_act_selberg2(A,B,C,S,F,V) {
return odiff_acts(odiff_op_selberg2(A,B,C,S,V),F,V);
}
def odiff_poly_solve_selberg2(A,B,C,S,N,V) {
return odiff_poly_solve(odiff_op_selberg2(A,B,C,S,V),N,V);
}
def odiff_poly_solve_selberg2_check(Sol,A,B,C,S,V) {
Lf_list = odiff_act_selberg2(A,B,C,S,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
def odiff_poly_solve_selberg2s(R_A,R_B,R_C,R_S,N,V) {
Sols = [];
for (S = R_S[0]; S <= R_S[1]; S += R_S[2]) {
if (S != 0) {
for (C = R_C[0]; C <= R_C[1]; C += R_C[2]) {
for (A = R_A[0]; A <= R_A[1]; A += R_A[2]) {
for (B = R_B[0]; B <= R_B[1]; B += R_B[2]) {
Sols = append(Sols, [[A,B,C,S,odiff_poly_solve_selberg2(A,B,C,S,N,V)]]);
}
}
}
}
}
return Sols;
}
def odiff_op_selberg2_1(A,B,K1,K2,V) {
if (A+B+K1+K2+1 == 0) {
error("A+B+K1+K2+1 == 0");
}
return odiff_op_selberg2(A,B,A+B+K2+1,1/(A+B+K1+K2+1),V);
}
def odiff_act_selberg2_1(A,B,K1,K2,F,V) {
return odiff_acts(odiff_op_selberg2_1(A,B,K1,K2,V),F,V);
}
def odiff_poly_solve_selberg2_1 (A,B,K1,K2,N,V) {
return odiff_poly_solve(odiff_op_selberg2_1(A,B,K1,K2,V),N,V);
}
def odiff_poly_solve_selberg2_1_check(Sol,A,B,K1,K2,V) {
Lf_list = odiff_act_selberg2_1(A,B,K1,K2,Sol[0],V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
def odiff_poly_solve_selberg2s_1 (R_A,R_B,R_K1,R_K2,N,V) {
Sols = [];
for (K1 = R_K1[0]; K1 <= R_K1[1]; K1 += R_K1[2]) {
for (K2 = R_K2[0]; K2 <= R_K2[1]; K2 += R_K2[2]) {
for (A = R_A[0]; A <= R_A[1]; A += R_A[2]) {
for (B = R_B[0]; B <= R_B[1]; B += R_B[2]) {
if (A+B+K1+K2+1 != 0) {
Sol = odiff_poly_solve_selberg2_1(A,B,K1,K2,N,V);
if (odiff_poly_solve_selberg2_1_check(Sol,A,B,K1,K2,V) == 1) {
Sols = append(Sols,[[A,B,K1,K2,Sol]]);
}
}
}
}
}
}
return Sols;
}
def odiff_poly_solve_selberg2s_1_half (R_AB,R_K1,R_K2,N,V) {
Sols = [];
for (K1 = R_K1[0]; K1 <= R_K1[1]; K1 += R_K1[2]) {
for (K2 = R_K2[0]; K2 <= R_K2[1]; K2 += R_K2[2]) {
for (A = R_AB[0]; A <= R_AB[1]; A += R_AB[2]) {
for (B = R_AB[0]; B <= R_AB[1]; B += R_AB[2]) {
if (A <= B && A+B+K1+K2+1 != 0) {
Sol = odiff_poly_solve_selberg2_1(A,B,K1,K2,N,V);
if (odiff_poly_solve_selberg2_1_check(Sol,A,B,K1,K2,V) == 1) {
Sols = append(Sols,[[A,B,K1,K2,Sol]]);
}
}
}
}
}
}
return Sols;
}
def odiff_lookup(X, List) {
if (type(List) != 4) {
error("parameter type error.");
}
for (I = 0; I < length(List); I++) {
if (X == List[I]) {
return 1;
}
}
return 0;
}
/*---------------------------*/
/* GKZ hypergeometric system */
/*---------------------------*/
/* return generators of toric ideal I_{A} */
/* [[u1,v1],[u2,v2],...] <- { x^u1-x^v1, x^u2-x^v2,...} */
def odiff_toric_ideal(A) {
if (type(A) == 6) {
A = omatrix_mtol(A);
}
if (type(A) != 4) {
error("parameter type error.");
}
A = omatrix_mtol(omatrix_trans(A));
Vx = [];
for (I = 0; I < length(A); I++) {
Vx = append(Vx, [strtov("x"+rtostr(I+1))]);
}
Vt = [];
for (I = 0; I < length(A[0]); I++) {
Vt = append(Vt, [strtov("t"+rtostr(I+1))]);
}
Vt0 = [t0];
V = append(Vt0, Vt);
V = append(V, Vx);
Ap = newmat(length(A), length(A[0]));
Am = newmat(length(A), length(A[0]));
for (I = 0; I < length(A); I++) {
for (J = 0; J < length(A[0]); J++) {
if (A[I][J] > 0) {
Ap[I][J] = A[I][J];
}
else {
Am[I][J] = -A[I][J];
}
}
}
FF = [];
for (I = 0; I < length(Vx); I++) {
F = Vx[I]*dp_dtop(dp_vtoe(Am[I]), Vt) - dp_dtop(dp_vtoe(Ap[I]), Vt);
FF = append(FF, [F]);
}
F = Vt0[0];
for (I = 0; I < length(Vt); I++) {
F *= Vt[I];
}
F -= 1;
FF = append(FF, [F]);
#if _DEBUG_Odiff_
print("FF: ", 0); print(FF);
#endif
GG = gr(FF, V, 2);
#if _DEBUG_Odiff_
print("V: ", 0); print(V);
print("GG: ", 0); print(GG);
#endif
Vt_all = append(Vt0, Vt);
Toric = [];
for (I = 0; I < length(GG); I++) {
Vars = vars(GG[I]);
Flag = 1;
for (J = 0; J < length(Vt_all); J++) {
if (odiff_lookup(Vt_all[J], Vars)) {
Flag = 0;
break;
}
}
if (Flag == 1) {
Toric = append(Toric, [GG[I]]);
}
}
#if _DEBUG_Odiff_
print("toric ideal: ", 0); print(Toric);
#endif
Toric_exp = [];
for (I = 0; I < length(Toric); I++) {
Dpoly = dp_ptod(Toric[I], Vx);
if (dp_hc(Dpoly) > 0) {
Dpoly_p = dp_hm(Dpoly);
Dpoly_m = Dpoly_p - Dpoly;
}
else {
Dpoly_m = dp_hm(Dpoly);
Dpoly_p = Dpoly_m - Dpoly;
}
Toric_exp = append(Toric_exp, [ [vtol(dp_etov(Dpoly_p)),vtol(dp_etov(Dpoly_m))] ]);
}
return Toric_exp;
}
/* A: matrix or list, Elements must be integer. */
/* B: vetor or list */
def odiff_op_gkz(A, B, V) {
if (type(A) == 6) {
A = omatrix_mtol(A);
}
if (type(B) == 5) {
B = vtol(B);
}
if (type(A) != 4 || type(B) != 4 || type(V) != 4) {
error("parameter type error.");
}
if ( length(A[0]) != length(V) || length(A) != length(B)) {
error("a length is wrong.");
}
L_list = [];
for (I = 0; I < length(A); I++) {
L = [];
for (J = 0; J < length(A[I]); J++) {
if (A[I][J] != 0) {
L = append(L, [[A[I][J]*V[J], vtol(dp_etov(dp_ptod(V[J],V)))]]);
}
}
L = append(L, [[-B[I], vtol(dp_etov(dp_ptod(1,V)))]]);
L_list = append(L_list, [L]);
}
/* toric */
Toric = odiff_toric_ideal(A);
for (I = 0; I < length(Toric); I++) {
L = [[1,Toric[I][0]],[-1,Toric[I][1]]];
L_list = append(L_list, [L]);
}
return L_list;
}
def odiff_act_gkz(A, B, F, V) {
return odiff_acts(odiff_op_gkz(A, B, V), F, V);
}
def odiff_poly_solve_gkz(A, B, N, V) {
return odiff_poly_solve(odiff_op_gkz(A, B, V), N, V);
}
def odiff_poly_solve_gkz_check(Sol, A, B, V) {
Lf_list = odiff_act_gkz(A, B, Sol[0], V);
for (I = 0; I < length(Lf_list); I++) {
if (Lf_list[I] != 0) {
return 0;
}
}
return 1;
}
end$
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ok_diff.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