/* $OpenXM: OpenXM/src/asir-contrib/packages/src/tk_hgpoly.rr,v 1.10 2018/09/25 01:35:23 takayama Exp $ */
/* Developed in h-mle/A-hg2/Prog */
#define USE_MODULE 1
import("names.rr")$
import("nk_toric.rr")$
import("ok_diff.rr")$
import("taka_diffop.rr")$
/* import("ot_hgm_ahg.rr"); by hand. mytoric() is called. */
/* To debug this file,
import("tk_fd.rr");
import("ot_hgm_ahg.rr");
and load this file.
*/
#if !USE_MODULE
extern P_A$
extern P_B$
extern TableOfFac$
extern TableOfFacSize$
extern Verbose$
#else
module tk_hgpoly;
static P_A$
static P_B$
static TableOfFac$
static TableOfFacSize$
static Verbose$
localf hgpoly$
localf check20$
localf check21$
localf check21b$
localf esum$
localf contingency3$
localf cprune$
localf feasible$
localf optip$
localf testip1$
localf testip2$
localf expectation$
localf gmvector$
localf test_exp1$
localf toric_in$
localf lfac$
localf distraction0$
localf distraction$
localf teststartc111$
localf isNonNegativeIntegerVector$
localf isNonNegativeIntegerVector0$
localf hgpoly_start$
localf teststartc111b$
localf base_subset_by_index$
localf kset$
localf multifac$
localf hgpolyk$
localf hgpolyk_v1$
localf testpolyc111$
localf upperBoundOfDegree$
localf facByTable$
localf assert_hgpolyk_1$
localf xpower$
localf dfac0$
localf dfac$
localf assert_hgpolyk_2$
localf verbose$
#endif
def verbose(F) {
Verbose = (F!=0)? 1: 0;
}
/* 2014.07.09 sst.ps (sst-ip.ps)
A must not contain negative number.
x^u/u! where u runs under the constraint Au=B.
*/
def hgpoly(A,B) {
Deg = 1+base_sum(B,0,0,length(B)-1);
if (type(getopt(deg)) >=0) Deg=eval_str(rtostr(getopt(deg)));
if (type(A) == 4) A=matrix_list_to_matrix(A);
D=size(A)[0];
N=size(A)[1];
Ap = matrix_transpose(A);
F=0;
Vx=[];
for (I=0; I<N; I++) {
Mon = 1;
for (J=0; J<D; J++) {
Mon = Mon*util_v(t,[J+1])^(Ap[I][J]);
if (Ap[I][J] < 0) error("A contains negative number.");
}
/* deg_1(Mon) >=1 */
F = F+util_v(x,[I+1])*Mon;
Vx = cons(util_v(x,[I+1]),Vx);
}
Vx=reverse(Vx);
/* print(print_input_form(poly_sort(F))); */
Fall = 1;
for (I=1; I<=Deg; I++) {
Fall += F^I; /* deg_1(each term of F^p) >= p */
}
// printf("Fall=%a\n",Fall);
P = coef(Fall,B[0],util_v(t,[1]));
for (I=1; I<length(B); I++) {
P = coef(P,B[I],util_v(t,[I+1]));
}
Pdist=dp_ptod(P,Vx);
if ( P ==0 ) return 0;
Pnew=0;
while (Pdist != 0) {
U = dp_ht(Pdist); U=dp_etov(U);
Degree=base_sum(U,0,0,length(U)-1);
Degree=fac(Degree);
Pnew += dp_hm(Pdist)/Degree;
Pdist = dp_rest(Pdist);
}
return [dp_dtop(Pnew,Vx),Pnew];
}
def check20() {
A=[[1,1,1,1],[0,1,0,1],[0,0,1,1]];
//B=[4,2,2];
// B=[6,4,2];
B=[10,8,6];
print(A);
print(B);
return hgpoly(A,B);
}
/* C111, h-mle/A-hg/pf-hg4.tex */
def check21() {
if (type(getopt(b)) < 0) {
// B=[4,4,4];
B=[1,1,1];
}else{
B=getopt(b);
}
At=[[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],[1,1,1]];
A=matrix_matrix_to_list(matrix_transpose(matrix_list_to_matrix(At)));
P_A=A; P_B=B;
print(A);
print(B);
P2=hgpoly(A,B); P=P2[0];
print(P);
MyRule=[]; Vx=[];
for (I=1; I<=7; I++) {
MyRule=cons([util_v(x,[I]),eval_str("x"+rtostr(I))],MyRule);
Vx=cons(eval_str("x"+rtostr(I)),Vx);
}
Vx=reverse(Vx);
Brule=[[b_1,B[0]],[b_2,B[1]],[b_3,B[2]]];
PP=base_replace(P,MyRule);
T=mytoric(A,[[dx1,1]]);
printf("Check if it satisfies A-hg system.\n");
for (I=0; I<length(T[0]); I++) {
printf("%a -> %a\n",I,odiff_act(base_replace(T[0][I],Brule),PP,Vx));
}
for (I=0; I<length(T[1]); I++) {
printf("%a -> %a\n",I,odiff_act(T[1][I],PP,Vx));
}
return [append(P2,T),A,B];
}
def check21b() {
F111=check21(|b=[1,1,1]);
F222=check21(|b=[2,2,2]);
F444=check21(|b=[4,4,4]);
return [F111[0][0],F222[0][0],F444[0][0]];
}
/* 2014.07.11 */
#define IDX(I,J,K) (I*Q*R+J*R+K)
def esum(L) {
return base_sum(L,0,0,length(L)-1);
}
def contingency3(P,Q,R) {
A=[];
N=P*Q*R;
for (I=0; I<P; I++) for (J=0; J<Q; J++) {
L = newvect(N);
for (K=0; K<R; K++) L[IDX(I,J,K)] = 1;
if (esum(L) > 1) A=cons(vtol(L),A);
}
for (I=0; I<P; I++) for (K=0; K<R; K++) {
L = newvect(N);
for (J=0; J<Q; J++) L[IDX(I,J,K)] = 1;
if (esum(L)>1) A=cons(vtol(L),A);
}
for (K=0; K<R; K++) for (J=0; J<Q; J++) {
L = newvect(N);
for (I=0; I<P; I++) L[IDX(I,J,K)] = 1;
if (esum(L)>1) A=cons(vtol(L),A);
}
return reverse(A);
}
def cprune(A) {
if (type(A) == 4) A=matrix_list_to_matrix(A);
Size=size(A);
D=Size[0];
Rank = matrix_rank(A);
if(Verbose) printf("Rank=%a\n",Rank);
while (Rank < D) {
Set = newvect(D);
for (I=0; I<D; I++) Set[I]=I;
Done=0;
for (I=0; I<D; I++) {
TrySet = base_set_minus(vtol(Set),[I]);
A2 = matrix_submatrix(A,TrySet);
if (matrix_rank(A2) == Rank) {
A=A2;
Done=1; D--;
break;
}
}
if (!Done) error("cprune failed.");
}
return(A);
}
/* Find a feasible point of AU=B
Calls mytoric in ot_hgm_ahg.rr
Returns 0 iff (A,B) is not feasible.
*/
def feasible(A,B) {
/* A copy from nk_toric.rr */
if (type(A) == 6) {
M = size(A)[0];
N = size(A)[1];
} else if (type(A) == 4) {
M = length(A);
N = length(A[0]);
}
if (M != length(B)) {
if (Verbose) print("A and B do not match.");
return 0;
}
/* List of variables VT=[t0, ..., tM],VX=[x1, ..., xN] */
VT = nk_toric.var_list("t", 0, M);
VX = nk_toric.var_list("x", 1, N);
L = [];
T = 1;
for (I = 0; I < M + 1; I++)
T *= VT[I];
T = T - 1;
L = cons(T, L);
for (J = 0; J < N; J++) {
Plus = 1;
Minus = 1;
for (I = 0; I < M; I++)
if (A[I][J] > 0)
Plus *= VT[I]^A[I][J];
else if (A[I][J] < 0)
Minus *= VT[I]^(-A[I][J]);
L = cons(VX[J]*Minus - Plus, L);
}
if(Verbose) printf("ideal : \n%a\n", L);
/* Computing GB with the elimination order t0, ..., tM > x1, ..., xM */
Ord = [[0, M + 1], [0, N]];
G = nd_gr(L, append(VT, VX), 0, Ord);
if(Verbose) printf("gb : \n%a\n", G);
/* Find a feasible point FPv */
FF=1;
for (I=0; I<M; I++) FF *= VT[I]^B[I];
FP = p_true_nf(FF,G,append(VT,VX),Ord)[0];
if(Verbose) printf("FP=%a\n",FP);
for (I=0; I<length(VT); I++) {
if (base_memberq(VT[I],vars(FP))) {
if(Verbose) print("Normal form contains t.");
return 0;
}
}
FPv=dp_ptod(FP,VX);
FPv=vtol(dp_etov(FPv));
/* intersection(G,K[VX]) */
Id = [];
for (I = 0; I < length(G); I++) {
Vars = vars(G[I]);
for (J = 0; J < length(Vars); J++) {
T = rtostr(Vars[J]);
if (sub_str(T, 0, 0) == "t")
break;
}
if (J == length(Vars))
Id = cons(G[I], Id);
}
return [FPv,FP,reverse(Id),VX];
}
/*
Find AU=B such that inner_product(U,W) is minimized.
W must have non-negative entries.
*/
def optip(A,B,W) {
F=feasible(A,B);
Start = F[1];
Id = F[2];
V = F[3];
Ord = poly_weight_to_omatrix(W,V);
G = nd_gr(Id,V,0,Ord);
OP = p_true_nf(Start,G,V,Ord)[0];
OPv=dp_ptod(OP,V);
OPv=vtol(dp_etov(OPv));
return(OPv);
}
def testip1() {
A = [[1,1,1,1],[0,1,2,3]];
B = [20,46];
print(feasible(A,B));
printf("A=%a\nB=%a\n",A,B);
return(optip(A,B,[1,1,0,0]));
}
/* 2014.12.12 */
def testip2() {
A =
[[1,0,0,1,1,0,1],
[0,1,0,1,0,1,1],
[0,0,1,0,1,1,1]];
B = [30,20,22];
print(feasible(A,B));
printf("A=%a\nB=%a\n",A,B);
return(optip(A,B,[1,1,1,1,1,1,1]));
}
/* 2014.12.20
X=[x_1, x_2, ..., x_n] is used to express hgpoly and random variable.
Expectation E[RandomVarible] at X=Pt for A-HG distribution for (A,B).
If fd=1, A is assumed to be FD.
*/
def expectation(A,B,Pt,RandomVariable)
"expectation(A,B,Pt,RandomVariable)\
Ex. expectation([[1,1,1],[0,1,2]],[5,6],[1,1,1/3],[x_1,x_2,x_3,x_1*x_2]);"
{
FD = getopt(fd);
if (type(FD)<0) FD=0;
GM = getopt(gm);
if (type(GM)<0) GM=0;
N = length(A[0]);
V=[]; DV=[];
for (I=1; I<=N; I++) {
V=cons(util_v(x,[I]),V);
DV=cons(util_v(dx,[I]),DV);
}
V = reverse(V); DV = reverse(DV);
Rule = assoc(V,Pt);
if (!FD) Z = hgpoly(A,B)[0];
else {
Mm = idiv(length(V),2);
X = newmat(2,Mm);
for (I=0; I<Mm; I++) {X[0][I]=V[I]; X[1][I]=V[Mm+I];}
R=cdr(B);
C=[base_sum(R,0,0,length(R)-1)-B[0],B[0]];
ABC = tk_fd.marginal2abc(C,R);
if (ABC[2] <= 0) error("C is negative for FD");
Z = tk_fd.fdah(ABC[0],ABC[1],ABC[2],X | approx=-ABC[0]+1)[0];
}
/* print(Z); */
if (type(RandomVariable)!=4) RandomVariable=[RandomVariable];
if (!GM) {
L = map(taka_dr1.euler2diff,RandomVariable,V,V,DV);
}else {
L=RandomVariable;
}
E = map(odiff_act,L,Z,V);
E = newvect(length(E),base_replace(E,Rule));
if (!GM) {
return(vtol(E/base_replace(Z,Rule)));
}else {
return(vtol(E));
}
}
def gmvector(A,B,P,Diff)
"gmvector(A,B,P,Diff)\
Ex. gmvector([[1,1,1],[0,1,2]],[5,6],[1,1,1/3],[1,dx_3,dx_1,dx_2]);"
{
return expectation(A,B,P,Diff | gm=1);
}
def test_exp1() {
X=[[1,1/2,1/3],[1,1,1]];
AB=tk_fd.abc2ahg(-3,[-2,-4],3);
E=tk_fd.expectation_abc(-3,[-2,-4],3,X);
printf("E by tk_fd.rr = %a\n",E);
A=AB[0]; B=AB[1];
E10 = expectation(A,B,base_flatten(X),[x_4,x_5,x_6]);
printf("E10=%a\n",E10);
}
/* construct hg polynomial by a basis of Ker A */
def toric_in(A) {
W=[[dx1, 1]]; /* dummy */
T0=mytoric(A,W);
E=T0[0];
T=T0[1];
V=T0[2]; Vb=T0[3]; Vd=add_d(V);
Gt=nd_gr(T,Vd,0,0);
dp_ord(0); /* rev lex */
T2=map(dp_ptod,Gt,Vd); T2=map(dp_ht,T2); T2=map(dp_dtop,T2,Vd);
return([T2,E,V,Vb]);
}
/*
In = toric_in([[1,1,1,1],[0,1,2,3]]);
*/
def distraction(F) {
if (type(F) > 3) return(map(distraction,F));
V = vars(F);
for (I=0; I<length(V); I++) F = distraction0(F,V[I]);
return(F);
}
def lfac(V,I) {
if (I == 0) return(1);
Ans=1;
for (J=0; J<I; J++) {
Ans = Ans*(V-J);
}
return(Ans);
}
def distraction0(F,V) {
Ans = 0;
N = deg(F,V);
for (I=0; I<=N; I++) {
Ans = coef(F,I,V)*lfac(V,I);
}
}
/*
In = toric_in([[1,1,1,1],[0,1,2,3]]);
V = In[2]; DV=tk_sm1emu.add_d(V);
L = base_replace(In[1],assoc(In[3],[4,7]));
L = base_replace(L,assoc(In[2],[1,1,1,1]));
Tdecomp = primadec(In[0],DV);
D=Tdecomp[0][0]; Base=Tdecomp[0][1];
Es = primadec(append(D,L),DV);
poly_solve_linear(Es[1][1],DV); --> all cooridnates are in N_0
construct series solution.
*/
def teststartc111() {
A = [[1,1,1,1,1,1,1,1],
[0,1,0,0,1,1,0,1],
[0,0,1,0,1,0,1,1],
[0,0,0,1,0,1,1,1]];
In = toric_in(A);
One = newvect(length(A[0]));
for (I=0; I<length(A[0]); I++) One[I] = 1;
One = vtol(One);
V = In[2]; DV=tk_sm1emu.add_d(V);
L = base_replace(In[1],assoc(In[3],[40,7,10,15]));
/* [4,7,10,15] is not feasible. Interesting case. Study later. */
L = base_replace(L,assoc(In[2],One));
Tdecomp = primadec(In[0],DV);
printf("Tdecomp=%a\n",Tdecomp);
for (I=0; I<length(Tdecomp); I++) {
D=Tdecomp[I][0]; Base=Tdecomp[I][1];
Es = primadec(append(D,L),DV);
for (J=0; J<length(Es); J++) {
E=poly_solve_linear(Es[J][1],DV);
printf("(I,J)=(%a,%a),Base=%a, E=%a\n",I,J,Base,E);
}
}
}
/* 2015.02.22 8:37 */
def isNonNegativeIntegerVector(E) {
// E=[[dx1,25],[dx2,0],[dx3,0],[dx4,5],[dx5,0],[dx6,0],[dx7,3],[dx8,7]]
N = length(E);
for (I=0; I<N; I++) {
if (type(E[I][1]) >=2) return(0);
if ( !number_is_integer(E[I][1])) return(0);
if (E[I][1] < 0) return(0);
}
return(1);
}
def isNonNegativeIntegerVector0(E) {
N = length(E);
for (I=0; I<N; I++) {
if (type(E[I]) >=2) return(0);
if ( !number_is_integer(E[I])) return(0);
if (E[I] < 0) return(0);
}
return(1);
}
def hgpoly_start(A,B) {
if (type(getopt(verbose))>=0) Verbose=getopt(verbose);
else Verbose=0;
In = toric_in(A);
if (Verbose) printf("Initial=%a\n",In);
One = newvect(length(A[0]));
for (I=0; I<length(A[0]); I++) One[I] = 1;
One = vtol(One);
V = In[2]; DV=tk_sm1emu.add_d(V);
L = base_replace(In[1],assoc(In[3],B));
L = base_replace(L,assoc(In[2],One));
Tdecomp = primadec(distraction(In[0]),DV);
if (Verbose) printf("Tdecomp=%a\n",Tdecomp);
for (I=0; I<length(Tdecomp); I++) {
D=Tdecomp[I][0]; Base=Tdecomp[I][1];
Es = primadec(append(D,L),DV);
for (J=0; J<length(Es); J++) {
Found=0;
E=poly_solve_linear(Es[J][1],DV);
if (Verbose) printf("(I,J)=(%a,%a),Base=%a, E=%a\n",I,J,Base,E);
if (isNonNegativeIntegerVector(E)) {
Found=1;
break;
}
}
if (Found) break;
}
if (!Found) {
printf("No integral starting point, but there might exists hg poly (Todo).\n");
return(0);
}
printf("(I,J)=(%a,%a),Base=%a, E=%a\n",I,J,Base,E);
Bindex=[];
for (I=0; I<length(Base); I++) {
Pos=base_position(vars(Base[I])[0],DV);
if (Pos>=0) Bindex=cons(Pos, Bindex);
}
Mstart=base_replace(DV,E);
Bindex=qsort(Bindex);
D = length(A); N = length(A[0]);
Amat=newmat(D,N,A);
Bmat=matrix_transpose(matrix_submatrix(matrix_transpose(Amat),Bindex));
return([Mstart,Amat,Bmat,Bindex,V]);
}
def teststartc111b() {
A = [[1,1,1,1,1,1,1,1],
[0,1,0,0,1,1,0,1],
[0,0,1,0,1,0,1,1],
[0,0,0,1,0,1,1,1]];
B = [40,7,10,15];
In = toric_in(A);
// return(hgpoly_start(A,B));
return(hgpolyk(A,B));
}
def base_subset_by_index(S,Idx) {
A = newvect(length(Idx));
for (I=0; I<length(Idx); I++) {
A[I] = S[Idx[I]];
}
return(vtol(A));
}
/*
cf. fd in tk_fd.rr
*/
def kset(N,Deg) {
for (I=0; I<N; I++) {
M=newvect(N); M[I]=1;
F = F+dp_vtoe(M);
}
E=dp_vtoe(newvect(N));
for (I=0; I<Deg; I++) {
E = E*F;
}
return(E);
}
/* defined in c111c-ser.rr */
def multifac(K) {
F=1;
for (I=0; I<length(K); I++) {
F = F*facByTable(K[I]);
if (F == 0) return(F);
}
return(F);
}
/* hgpolyk by a basis of Ker A. It returns distributed polynomial.
factorials are evaluated from scratch, (then it is less efficient).
The case A=[[1,1,1],[0,1,2]]; B=[3,3]; has not been implemented.
--> Fixed. 2015.02.24.
hgpoly_start(A,B|verbose=1); returns fractional exponents.
B=A*[p,0,q] case is fine.
Example:
A=[[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]];
tk_hgpoly.hgpolyk(A,[9,3,3,4]|verbose=1);
tk_hgpoly.hgpolyk(A,[9,3,3,4]|std=[1,dx5,dx6],x=[1,1/2,1/3,1,1,1]);
*/
def hgpolyk(A,B) {
Approx=-1;
if (type(getopt(approx)) <0) {
Approx=upperBoundOfDegree(A,B);
}else {
Approx=getopt(approx);
}
if (Approx < 0) {
printf("Error: the upper bound of the degree cannot be determined automatically. Give it by the approx=m option\n");
return(0);
}
if (type(getopt(verbose)) >0) Verbose=getopt(verbose); else Verbose=0;
if (type(getopt(x)) >0) X=getopt(x); else X=0;
if (type(getopt(std))>0) Std=getopt(std); else Std=[newvect(length(A[0]))];
S=hgpoly_start(A,B);
if (S == 0) {
return(0);
}
Mstart=S[0]; Amat=S[1]; Bmat=S[2]; Idind=S[3]; /* ind = independent */
Vorig=S[4];
Std = matrix_list_to_matrix(Std);
Dvars = add_d(Vorig);
if (Verbose) printf("Dvars=%a\n",Dvars);
for (I=0; I<length(Std); I++) {
if (type(Std[I])<4) {
Std[I] = dp_etov(dp_ptod(Std[I],Dvars));
}
Std[I] = matrix_list_to_matrix(Std[I]);
}
if (Verbose) printf("Std=%a\n",Std);
N = length(Vorig); D=length(B);
Idall = newvect(N); V=newvect(N);
for (I=0; I<N; I++) {Idall[I]=I; V[I] = util_v(x,[I]);}
Idall = vtol(Idall); V=vtol(V);
Iddep=base_set_minus(Idall,Idind); /* dep = dependent */
Eq=Amat*newvect(N,V)-newvect(D,B);
Eq=vtol(Eq);
Vdep=base_subset_by_index(V,Iddep);
Vind=base_subset_by_index(V,Idind);
Rule = poly_solve_linear(Eq,Vdep);
Zvalue = newvect(length(Std));
for (I=0; I<=Approx; I++) {
K = kset(length(Vind),I);
while (K != 0) {
U0 = vtol(dp_etov(dp_ht(K))); K = dp_rest(K);
U = base_replace_n(base_replace(V,Rule),assoc(Vind,U0));
if (isNonNegativeIntegerVector0(U)) {
if (Verbose) printf("U=%a\n",U);
U = newvect(N,U);
C = multifac(U);
for (J=0; J<length(Std); J++) {
if (C != 0) {
Cj=1/C;
Diff = Std[J];
Cj = Cj*dfac(U,Diff);
if (X != 0) {
Zvalue[J] = Zvalue[J]+Cj*xpower(X,U-Diff);
}else {
Zvalue[J] = Zvalue[J]+Cj*dp_vtoe(U-Diff);
}
}
}
}
}
}
if (length(Zvalue)==1) return(Zvalue[0]);
else return(vtol(Zvalue));
}
def hgpolyk_v1(A,B) {
Approx=-1;
if (type(getopt(approx)) <0) {
Approx=upperBoundOfDegree(A,B);
}else {
Approx=getopt(approx);
}
if (Approx < 0) {
printf("Error: the upper bound of the degree cannot be determined automatically. Give it by the approx=m option\n");
return(0);
}
if (type(getopt(x)) >0) X=getopt(x); else X=0;
if (type(getopt(std))>0) Std=getopt(std); else Std=0;
S=hgpoly_start(A,B);
if (S == 0) {
return(0);
}
Mstart=S[0]; Amat=S[1]; Bmat=S[2]; Idind=S[3]; /* ind = independent */
Vorig=S[4];
N = length(Vorig); D=length(B);
Idall = newvect(N); V=newvect(N);
for (I=0; I<N; I++) {Idall[I]=I; V[I] = util_v(x,[I]);}
Idall = vtol(Idall); V=vtol(V);
Iddep=base_set_minus(Idall,Idind); /* dep = dependent */
Eq=Amat*newvect(N,V)-newvect(D,B);
Eq=vtol(Eq);
Vdep=base_subset_by_index(V,Iddep);
Vind=base_subset_by_index(V,Idind);
Rule = poly_solve_linear(Eq,Vdep);
F=0;
for (I=0; I<=Approx; I++) {
K = kset(length(Vind),I);
while (K != 0) {
U0 = vtol(dp_etov(dp_ht(K))); K = dp_rest(K);
U = base_replace_n(base_replace(V,Rule),assoc(Vind,U0));
if (isNonNegativeIntegerVector0(U)) {
C = multifac(U);
if (X != 0) {
if (C != 0) F = F+(1/C)*xpower(X,U);
}else {
if (C != 0) F = F+(1/C)*dp_vtoe(newvect(N,U));
}
}
}
}
return(F);
}
def testpolyc111() {
A = [[1,1,1,1,1,1,1,1],
[0,1,0,0,1,1,0,1],
[0,0,1,0,1,0,1,1],
[0,0,0,1,0,1,1,1]];
// B = [40,7,10,15];
B=[16,8,8,8];
return(hgpolyk(A,B));
}
/* 2015.02.23, 6:15 */
def upperBoundOfDegree(A,B) {
D=length(A); N=length(A[0]);
/* Two heuristic upper bound */
for (I=0; I<D; I++) {
Found=1;
for (J=0; J<N; J++) {
if (A[I][J] <= 0) { Found=0; break; }
}
if (Found) return(B[I]);
}
/* */
Found = 1;
for (J=0; J<N; J++) {
S=0;
for (I=0; I<D; I++) {
S += A[I][J];
}
if (S<=0) { Found=0; break; }
}
if (Found) return(base_sum(B,0,0,length(B)-1));
return(-1);
}
def facByTable(N) {
if (TableOfFacSize <=0) {
TableOfFacSize=50;
TableOfFac=newvect(TableOfFacSize);
}
if (N<0) return(0);
if (N >= TableOfFacSize) {
T=newvect(N*2);
for (I=0; I<TableOfFacSize; I++) {
T[I] = TableOfFac[I];
}
TableOfFac=T; TableOfFacSize=2*N;
}
if (TableOfFac[N] <= 0) {
TableOfFac[N] = fac(N);
}
return(TableOfFac[N]);
}
/* import("tk_fd.rr"); but it might cause an error. */
def assert_hgpolyk_1() {
A=-5; B=[-6]; C=4;
F=tk_fd.fdah(A,B,C,[[x1,x2],[x3,x4]] | approx=6);
F=F[0]*F[3];
G=hgpolyk([[0,0,1,1],[1,0,1,0],[0,1,0,1]],[C-1-B[0],C-1-A,-B[0]]);
G=dp_dtop(G,[x1,x2,x3,x4]);
printf("2x2, F-G=%a\n",F-G);
A=-7; B=[-4,-5]; C=3;
F=tk_fd.fdah(A,B,C,[[x1,x2,x3],[x4,x5,x6]] | approx=7);
F=F[0]*F[3];
G=hgpolyk([[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
[C-1-B[0]-B[1],C-1-A,-B[0],-B[1]]);
G=dp_dtop(G,[x1,x2,x3,x4,x5,x6]);
printf("2x3, F-G=%a\n",F-G);
}
/* 2015.02.23, 15:41 */
def xpower(X,U) {
N = length(X);
C=1;
for (I=0; I<N; I++) {
C = C*X[I]^U[I];
}
return(C);
}
def dfac0(U,D) {
Dfac=1;
for (I=0; I<D; I++) {
Dfac = Dfac*(U-I);
}
return(Dfac);
}
def dfac(U,Diff) {
Dfac=1;
N=length(U);
for (I=0; I<N; I++) Dfac=Dfac*dfac0(U[I],Diff[I]);
return(Dfac);
}
def assert_hgpolyk_2() {
A=-7; B=[-4,-5]; C=3;
V=[x1,x2,x3,x4,x5,x6];
F=tk_fd.fdah(A,B,C,[[x1,x2,x3],[x4,x5,x6]] | approx=7);
F=F[0]*F[3];
G=hgpolyk([[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
[C-1-B[0]-B[1],C-1-A,-B[0],-B[1]]);
G=dp_dtop(G,V);
printf("2x3, F-G=%a\n",F-G);
G=hgpolyk([[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
[C-1-B[0]-B[1],C-1-A,-B[0],-B[1]] | std=[dx5,dx5^2,dx6], x=V);
H=[diff(F,x5)-G[0],diff(diff(F,x5),x5)-G[1],diff(F,x6)-G[2]];
printf("differential->%a\n",H);
Xval=[1,1/2,1/3,1,1,1];
G=hgpolyk([[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
[C-1-B[0]-B[1],C-1-A,-B[0],-B[1]] | std=[dx5,dx5^2,dx6], x=Xval);
H=[diff(F,x5)-G[0],diff(diff(F,x5),x5)-G[1],diff(F,x6)-G[2]];
H=base_replace(H,assoc(V,Xval));
printf("differential by num->%a\n",H);
}
#if !USE_MODULE
#else
endmodule$ // endmodule of tk_hgpoly
#endif
end$
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tk_hgpoly.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