[BACK]Return to tk_fd.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / packages / src

Diff for /OpenXM/src/asir-contrib/packages/src/tk_fd.rr between version 1.3 and 1.12

version 1.3, 2014/06/29 00:53:37 version 1.12, 2014/12/12 08:12:55
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir-contrib/packages/src/tk_fd.rr,v 1.2 2014/06/04 23:44:17 takayama Exp $ */  /* $OpenXM: OpenXM/src/asir-contrib/packages/src/tk_fd.rr,v 1.11 2014/08/27 02:14:40 takayama Exp $ */
 /*  /*
   2015.05.23   h-mle/FD/Prog/fdpf.rr --> tk_fd.rr    2015.05.23   h-mle/FD/Prog/fdpf.rr --> tk_fd.rr
  */   */
Line 145  localf ygdmat_abc$
Line 145  localf ygdmat_abc$
 localf ygfvect_poly$  localf ygfvect_poly$
 localf ygahvec$  localf ygahvec$
 localf ygtry15$  localf ygtry15$
   localf marginal2abc$
   localf abc2marginal$
   localf ygQmat$
   localf fdA$
   localf hgpoly_fd$
   localf hgpoly_fd_beta$
   localf hgpoly_fd0$
   localf hgpoly_fd_beta0$
   localf check21$
   localf checkahg$
   localf abc2alpha$
   localf ygNmat2$
   localf ygCmat2$
   localf ygRmat2$
   localf ygQmat2$
   localf beta2marginal$
   localf beta2abc$
   localf atofd$
   localf atofd_beta$
   localf check22$
   localf checkfd$
   localf ygDk$
   localf ygDk2$
   localf containNegative$
   localf check23$
   localf ygDc2$
   localf ygDc$
   localf check24$
   localf ygDa2$
   localf ygDa$
   localf check25$
   localf poly_fd$
   localf poly_fd_beta$
   localf poly_fdvec$
   localf poly_fdvec_beta$
   localf check26$
   localf ygDca2$
   localf hgpoly_fd_beta1$
   localf check27$
   localf hgpoly_fd1$
   localf fvect_poly2$
   localf check28$
   localf check28b$
   localf ahvec_beta$
   localf ahvec_abc$
   localf expectation_abc$
 #endif  #endif
 /* Matsumoto, page 3 */  /* Matsumoto, page 3 */
   
Line 1408  def fdah_poly(A,B,C,Y) {
Line 1454  def fdah_poly(A,B,C,Y) {
 2014.05.15  2014.05.15
  */   */
 /*&usage begin:  /*&usage begin:
   ahvec(A,B,C,Y | norecc=y, approx=n)    ahvec(A,B,C,Y | norecc=n, approx=n)
   It returns [F,Gamma]. F*Gamma is equal to    It returns R=[F,Gamma]. F*Gamma is equal to
   (dy_21 Z, dy_22 Z, ..., dy_2p Z) at y = Y where p = length(B)+1.    (dy_21 Z, dy_22 Z, ..., dy_2p Z) at y = Y where p = length(B)+1.
   y = [[y_11,y_12, ..., y_1p],[y_21,y_22,...,y_2p]].    y = [[y_11,y_12, ..., y_1p],[y_21,y_22,...,y_2p]].
   y_2* must not be 0.    y_2* must not be 0.
     When the option all=1,  R[2]*R[1] is Z (normalizing constant or hg polynomial).
   example:    example:
   A=-3;    A=-3;
   B=[-2,-5];    B=[-2,-5];
   C=2;    C=2;
   Yval=[[1,1/2,1/3],[1,1,1]];    Yval=[[1,1/2,1/3],[1,1,1]];
   ahvec(A,B,C,Yval);    ahvec(A,B,C,Yval);
     example:
      F=ahvec(-1,[-2],3,[[x11,x12],[x21,x22]] | all=1);
      red(F[2]*F[1]);
      returns hg polynomial  sum(x^u/u!)  (where Au=Beta).
   end: */    end: */
 def ahvec(A,B,C,Y) {  def ahvec(A,B,C,Y) {
   Opt=getopt();    Opt=getopt();
   if (type(getopt(norecc)) >0) UseRecc=0; else UseRecc=1;    if (type(getopt(norecc)) >0) UseRecc=0; else UseRecc=1;
     if (type(getopt(all)) > 0) All=1; else All=0;
   M = length(B);    M = length(B);
   E = newmat(2,M+1);    E = newmat(2,M+1);
   E[0][0] = -A; E[1][0] = C-1;    E[0][0] = -A; E[1][0] = C-1;
Line 1448  def ahvec(A,B,C,Y) {
Line 1500  def ahvec(A,B,C,Y) {
     Y2j1 = Y[1][J]; E2j1 = E[1][J];      Y2j1 = Y[1][J]; E2j1 = E[1][J];
     Ahvec[J] = (-Fd[J]+E2j1*Fd[0])*Ye/Y2j1;      Ahvec[J] = (-Fd[J]+E2j1*Fd[0])*Ye/Y2j1;
   }    }
     if (All) { Zvalue=Ye*F[0];}
   
   NoGamma = getopt(nogamma);    NoGamma = getopt(nogamma);
   if (type(NoGamma)<0) NoGamma=0;    if (type(NoGamma)<0) NoGamma=0;
    /* Gamma factors, todo double check */     /* Gamma factors, todo double check */
   if (NoGamma) {    if (NoGamma) {
      Gamma = "[1/gamma(e+1)]";       Gamma = "[1/gamma(e+1)]";
   }else{    }else if (C > 0) {
     Gamma=tk_number_invgamma(-A+1)*tk_number_invgamma(C-1+1);      Gamma=tk_number_invgamma(-A+1)*tk_number_invgamma(C-1+1);
     for (I=0; I<M; I++) {      for (I=0; I<M; I++) {
       Gamma=Gamma*tk_number_invgamma(-B[I]+1);        Gamma=Gamma*tk_number_invgamma(-B[I]+1);
     }      }
   }    }else { Gamma=1; } /* because fvect_poly2 is called. */
   return [vtol(Ahvec),Gamma];    R=[vtol(Ahvec),Gamma];
     if (All) return append(R,[Zvalue]);
     else return R;
 }  }
 def check17() {  def check17() {
   A=-3;    A=-3;
Line 1485  def check17() {
Line 1540  def check17() {
   return taka_base_equal(Favec,Favec3);    return taka_base_equal(Favec,Favec3);
 }  }
 def fvect_poly(A,B,C,X) {  def fvect_poly(A,B,C,X) {
     if (C <= 0) return(fvect_poly2(A,B,C,X));
   EvalDmat=1;    EvalDmat=1;
   if (type(X) == 4) X=newvect(length(X),X);    if (type(X) == 4) X=newvect(length(X),X);
   if (taka_base_equal(B,FD_BB) &&    if (taka_base_equal(B,FD_BB) &&
Line 1864  def  ah_init_value(A,B,C,Y)  {
Line 1920  def  ah_init_value(A,B,C,Y)  {
   return [Amat,V,BRule,0,Val,0,Base,Fexact];    return [Amat,V,BRule,0,Val,0,Base,Fexact];
 }  }
   
 /* 2014.06.23.  Ref: @s/2014/06/24-intersection-sabun.pdf by Y.Goto , 25-my-note-fdpf-yg-funcs.pdf */  /* 2014.06.23.  Ref: @s/2014/06/24-intersection-sabun.pdf by Y.Goto , 25-my-note-fdpf-yg-funcs.pdf
      2014/07/14-fdpf-memo-ygfuncs.pdf   mynote in goto's note with yg-function names.
   */
 def ygNmat() {  def ygNmat() {
   N=newmat(MM+1,MM+1);    N=newmat(MM+1,MM+1);
   for (I=0; I<MM+1; I++) for (J=0; J<MM+1; J++) N[I][J] = 1;    for (I=0; I<MM+1; I++) for (J=0; J<MM+1; J++) N[I][J] = 1;
Line 1914  def ygdmat_abc(A,B,C) {
Line 1972  def ygdmat_abc(A,B,C) {
 }  }
   
 def ygfvect_poly(A,B,C,X) {  def ygfvect_poly(A,B,C,X) {
     if (C <=0 ) return(fvect_poly2(A,B,C,X));
   YgEvalDmat=1;    YgEvalDmat=1;
   if (type(X) == 4) X=newvect(length(X),X);    if (type(X) == 4) X=newvect(length(X),X);
   if (taka_base_equal(B,YgFD_BB) &&    if (taka_base_equal(B,YgFD_BB) &&
Line 1947  def ygfvect_poly(A,B,C,X) { 
Line 2006  def ygfvect_poly(A,B,C,X) { 
 def ygahvec(A,B,C,Y) {  def ygahvec(A,B,C,Y) {
   Opt=getopt();    Opt=getopt();
   if (type(getopt(norecc)) >0) UseRecc=0; else UseRecc=1;    if (type(getopt(norecc)) >0) UseRecc=0; else UseRecc=1;
     if (type(getopt(all)) > 0) All=1; else All=0;
   M = length(B);    M = length(B);
   E = newmat(2,M+1);    E = newmat(2,M+1);
   E[0][0] = -A; E[1][0] = C-1;    E[0][0] = -A; E[1][0] = C-1;
Line 1958  def ygahvec(A,B,C,Y) {
Line 2018  def ygahvec(A,B,C,Y) {
   Ye = 1;    Ye = 1;
   for (J=0; J<=M; J++) Ye = Ye*(Y[0][J]^E[0][J])*(Y[1][J]^E[1][J]);    for (J=0; J<=M; J++) Ye = Ye*(Y[0][J]^E[0][J])*(Y[1][J]^E[1][J]);
   if (number_is_integer(A) && A < 0 && UseRecc) F=ygfvect_poly(A,B,C,X);    if (number_is_integer(A) && A < 0 && UseRecc) F=ygfvect_poly(A,B,C,X);
   else F=ygfvect(A,B,C,X | option_list=Opt);  /* not implemented.*/    else F=fvect(A,B,C,X | option_list=Opt);
   Fd=newvect(M+1);    Fd=newvect(M+1);
   Fd[0] = F[0];  /* Fd are Euer derivatives of F_D */    Fd[0] = F[0];  /* Fd are Euer derivatives of F_D */
   for (I=1; I<=M; I++) {    for (I=1; I<=M; I++) {
Line 1972  def ygahvec(A,B,C,Y) {
Line 2032  def ygahvec(A,B,C,Y) {
     Y2j1 = Y[1][J]; E2j1 = E[1][J];      Y2j1 = Y[1][J]; E2j1 = E[1][J];
     Ahvec[J] = (-Fd[J]+E2j1*Fd[0])*Ye/Y2j1;      Ahvec[J] = (-Fd[J]+E2j1*Fd[0])*Ye/Y2j1;
   }    }
     if (All) { Zvalue=Ye*F[0];}
   
   NoGamma = getopt(nogamma);    NoGamma = getopt(nogamma);
   if (type(NoGamma)<0) NoGamma=0;    if (type(NoGamma)<0) NoGamma=0;
    /* Gamma factors, todo double check */     /* Gamma factors, todo double check */
   if (NoGamma) {    if (NoGamma) {
      Gamma = "[1/gamma(e+1)]";       Gamma = "[1/gamma(e+1)]";
   }else{    }else if (C>0) {
     Gamma=tk_number_invgamma(-A+1)*tk_number_invgamma(C-1+1);      Gamma=tk_number_invgamma(-A+1)*tk_number_invgamma(C-1+1);
     for (I=0; I<M; I++) {      for (I=0; I<M; I++) {
       Gamma=Gamma*tk_number_invgamma(-B[I]+1);        Gamma=Gamma*tk_number_invgamma(-B[I]+1);
     }      }
   }    }else { Gamma=1; } /* ygfvect_poly2 is called. */
   return [vtol(Ahvec),Gamma];    R= [vtol(Ahvec),Gamma];
     if (All) return append(R,[Zvalue]);
     else return R;
 }  }
   
   
Line 2007  def ygtry15(M) {
Line 2070  def ygtry15(M) {
   T=Fam - base_replace(Dmat,[[a,Aval]])*Fa;    T=Fam - base_replace(Dmat,[[a,Aval]])*Fa;
   printf("If %a is not zero vector , then there is an error.\n",T);    printf("If %a is not zero vector , then there is an error.\n",T);
   return Dmat;    return Dmat;
   }
   
   
   /* cf. fdah()  */
   def marginal2abc(RSum,CSum) {
     if (length(RSum) != 2) error("length of RSum must be 2.");
     Size = length(CSum);
     CSum=newvect(Size,CSum);
     S = RSum[0]+RSum[1];
     if (S != base_sum(CSum,0,0,Size-1)) {
       T=S-base_sum(CSum,0,0,Size-2);
       printf("Warning. RSum != CSum. The last value of CSum is changed from %a to %a.\n",CSum[Size-1],T);
       CSum[Size-1] = T;
       if (T < 0) error("Changed value is negative.");
     }
     CSum = vtol(CSum);
     B=vtol(-newvect(Size-1,cdr(CSum)));
     C=RSum[1] + 1 + base_sum(B,0,0,Size-2);
     A=-RSum[0];
     if ((C < 0) && MyDebug) printf("Advise: C=%a is negative. Exchange columns so that a big value comes to the first and Exchange rows so that a big value come to the second. \n",C);
     return [A,B,C];
   
   }
   def abc2marginal(A,B,C) {
     R2=C-1-base_sum(B,0,0,length(B)-1);
     R1=-A;
     RSum=[R1,R2];
     C1=-A+C-1;
     CSum = -newvect(length(B),B);
     CSum = cons(-A+C-1,vtol(CSum));
     return [RSum,CSum];
   }
   // check20 in A-hg/Prog/hgpoly.rr
   
   def ygQmat(X,K) {
     return ygRmat(X,K); /* Implement without global variables. */
   }
   
   /* 2014.08.03 from hgpoly.rr, should use hgpoly.rr with deg=...
      2 x (M+1)
   */
   def fdA(M) {
     A = newmat(M+1+1,2*(M+1));
     for (I=0; I<M+1;I++) A[0][I] =1;
     for (I=0; I<M+1; I++) {
       A[I+1][I] = 1;
       A[I+1][M+1+I] = 1;
     }
     return A;
   }
   def hgpoly_fd(A,B,C) {
     Marginal = abc2marginal(A,B,C);
     if (containNegative(Marginal)) error("Negative marginal");
     M = length(B);
     Beta =cons(Marginal[0][0],Marginal[1]);
     return hgpoly_fd_beta(M,Beta);
   }
   def hgpoly_fd0(A,B,C) {
     Marginal = abc2marginal(A,B,C);
     if (containNegative(Marginal)) error("Negative marginal");
     M = length(B);
     Beta =cons(Marginal[0][0],Marginal[1]);
     return hgpoly_fd_beta0(M,Beta);
   }
   /*
     Table. B[0] is the first row sum.
     B[0] |
     CB   |
            B[1], B[2], ..., B[M+1]
    */
   def hgpoly_fd_beta(M,B) {
     CB = (base_sum(B,0,0,length(B)-1)-B[0])-B[0];
     if (CB == 1) return hgpoly_fd_beta1(M,B);
     else return hgpoly_fd_beta0(M,B);
   }
   /* _beta0  is the bruce force method. cf _beta1 */
   def hgpoly_fd_beta0(M,B) {
     /* printf("hgpoly_fd_beta0(bruce force), M=%a, Beta=%a\n",M,B); */
     A = fdA(M);
     Deg  = base_sum(B,0,0,length(B)-1)-B[0];
     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,Vx];
   }
   /* from Ogawa data. */
   def check21() {
     B=[2,1,1,1,1]; print(B);
     F=hgpoly_fd_beta(3,B); print(F);
     B=[10,8,6,8]; print(B);
     F=hgpoly_fd_beta(2,B);  print(F);
     B=[4,3,1,3,1,1]; print(B);
     F=hgpoly_fd_beta(4,B);  print(F);
     B=[5,4,2,3,2,1]; print(B);
     F=hgpoly_fd_beta(4,B);  print(F);
     B=[6,5,4,1,1,1]; print(B);
     F=hgpoly_fd_beta(4,B);  print(F); /* takes a few seconds */
   }
   
   /*
      Beta=[3,6,3,4];
      F=hgpoly_fd_beta(2,Beta);
      checkahg(F[0],Beta,F[2]);
    */
   def checkahg(F,Beta,V) {
     N=(length(Beta)-1)*2;  M=length(Beta)-2;
     if (N != length(V)) error("N != length(V)");
     DV = newvect(N) ;
     for (I=0; I<N; I++) DV[I] = eval_str("d"+rtostr(V[I]));
     EV = newvect(N);
     for (I=0; I<N; I++) EV[I] = V[I]*DV[I];
     A = fdA(M);
     E = A*EV;
     for (I=0; I<length(Beta);I++) {
       if (red(tkdiff_act(E[I]-Beta[I],F,V)) != 0) {print("No"); return(0); }
     }
     for (I=0; I<M+1; I++) {
       for (J=I+1; J<M+1; J++) {
         T=DV[I]*DV[M+1+J]-DV[J]*DV[M+1+I];
         if (red(tkdiff_act(T,F,V)) != 0) {print("No"); return(0); }
       }
     }
     return(1);
   }
   
   /* B[0] is not a dummy. cf. in the setparam(), B[0] is a dummy. */
   def abc2alpha(A,B,C) {
     M=length(B);
     Alph = newvect(M+3);
     Alph[0]=-C+base_sum(B,0,0,M-1);
     for (I=1; I<=M; I++) Alph[I] = -B[I-1];
     Alph[M+1] = C-A;
     Alph[M+2] = A;
     return(vtol(Alph));
   }
   
   def ygNmat2(M) {
     N=newmat(M+1,M+1);
     for (I=0; I<M+1; I++) for (J=0; J<M+1; J++) N[I][J] = 1;
     return N;
   }
   
   def ygCmat2(A,B,C) {
     M = length(B);
     Alp = abc2alpha(A,B,C);
     L=newvect(M+1);
     L[0] = 1/Alp[M+2];
     for (I=1; I<=M; I++) L[I] = 1/Alp[I];
     return (ygNmat2(M)/Alp[M+1]) + matrix_diagonal_matrix(L);
   }
   
   def ygRmat2(A,B,C,X,K) {
     M = length(B);
     Alp = abc2alpha(A,B,C);
     X2=newvect(M+2);
     X2[0]=0;
     for (I=1; I<=M; I++) X2[I] = X[I-1];
     X2[M+1]=1;  /* cf. def xx(I) */
     X=X2;
   
     Lx = newvect(M+1);
     for (I=1; I<=M; I++) Lx[I] = 1-X[I];
     Lx = matrix_diagonal_matrix(Lx);
     L1 = newvect(M+1); L1[0]=1;
     L1 = matrix_diagonal_matrix(L1);
     N = ygNmat2(M);
     R=N*(1-X[K])/Alp[M+1] + Lx*N*L1/Alp[M+2] - L1*N*Lx/(1-Alp[M+2]);
   
     L = newvect(M+1);
     AX=0; for (I=0; I<=M+1; I++) AX += Alp[I]*X[I];
     L[0] = (1-X[K])/Alp[M+2] - (1/(1-Alp[M+2]))*(AX/Alp[M+2] + 1);
     for (I=1; I<=M; I++) {
       L[I] = (X[I]-X[K])/Alp[I];
     }
     return (R + matrix_diagonal_matrix(L));
   }
   
   def ygQmat2(A,B,C,X,K) {
     return ygRmat2(A,B,C,X,K);
   }
   
   def beta2marginal(Beta) {
     R=cdr(Beta);
     C=[Beta[0],base_sum(R,0,0,length(R)-1)-Beta[0]];
     return([C,R]);
   }
   def beta2abc(Beta) {
     Marginal=beta2marginal(Beta);
     return(marginal2abc(Marginal[0],Marginal[1]));
   }
   
   /* From A-hg series  to FD. */
   def atofd(F,A,B,C,Z) {
     M = length(B);
     E = newmat(2,M+1);
     Y = newmat(2,M+1);
     Y[0][0] = util_v(y,[0,0]);
     for (J=0; J<=M; J++) Y[1][J] = util_v(y,[1,J]);
     E[0][0] = -A; E[1][0] = C-1;
     for (J=1; J<=M; J++) E[1][J] = -B[J-1];
     X=newvect(M);
     for (I=0; I<M; I++) {
       X[I] =util_v(x,[I+1]);
     }
     for (I=0; I<M; I++) {
       /*     X[I] = Y[0][I+1]*Y[1][0]/(Y[0][0]*Y[1][I+1]); */
       Y[0][I+1] = X[I]/(Y[1][0]/(Y[0][0]*Y[1][I+1]));
     }
     Y2 = newvect(2*(M+1));
     for (J=0; J<=M; J++) {Y2[J] = Y[0][J]; Y2[M+1+J]=Y[1][J];}
     Ye = 1;
     for (J=0; J<=M; J++) Ye = Ye*(Y[0][J]^E[0][J])*(Y[1][J]^E[1][J]);
     Rule = assoc(Z,vtol(Y2));
     if (MyDebug) printf("F=%a,\nYe=%a,\nRule=%a\n",F,Ye,Rule);
     Fd = red(base_replace(F/Ye,Rule));
     return([Fd,vtol(X)]);
   }
   /*
   Beta=[2,1,1,1,1];
   FdA=hgpoly_fd_beta(3,Beta);
   atofd_beta(FdA[0],Beta,FdA[2]);
     Beta is the marginal sum.
     FdA[0] is the polynomial.
     FdA[1] is dp_ptod(FdA[0]);
     FdA[2] is the list of variables.
   */
   def  atofd_beta(F,Beta,Z) {
     Marginal = beta2marginal(Beta);
     CSum=Marginal[1];
     RSum=Marginal[0];
     ABC=marginal2abc(RSum,CSum);
     return(atofd(F,ABC[0],ABC[1],ABC[2],Z));
   }
   /*
   check22(-2,[-1,-2,-1],3);
    */
   def check22(A,B,C) {
     M = length(B);
     Beta=newvect(M+2);
     Beta[0]=-A;
     Beta[1]=C-1-A;
     for (I=0; I<M; I++) Beta[2+I]=-B[I];
     /* Beta=[-A,C-1-A,-B[0],-B[1],-B[2]]; */
     F=hgpoly_fd_beta(M,Beta); print(F);
     G=atofd(F[0],A,B,C,F[2]);
     R=checkfd(G[0],A,B,C,G[1]); printf("R=[0,...,0]? ===> R=%a\n",R);
     return(G);
   }
   
   /* 2014.08.05. Does F satisfy the system of equations for F_D? */
   def checkfd(F,A,B,C,V) {
     M = length(B);
     DV = [];
     for (I=0; I<M; I++) DV=cons(eval_str("d"+rtostr(V[I])),DV);
     DV = reverse(DV);
     T=0; for (I=0; I<M; I++) T = T+V[I]*DV[I];
     Ans=[];
     for (I=0; I<M; I++) {
       G=odiff_act(V[I]*DV[I],odiff_act(T+C-1,F,V),V)
         -V[I]*odiff_act(T+A,odiff_act(V[I]*DV[I]+B[I],F,V),V);
       Ans = cons(G,Ans);
     }
     return(reverse(Ans));
   }
   
   /* F(b-e_k) = Dk F(b)
      e_1, e_2, ..., but b_1=B[0]
    */
   def ygDk(A,B,C,X,K) {
     Qk = ygQmat2(A+1,B,C+1,X,K);
     Q0 = ygQmat2(A+1,B,C+1,X,0);
     return map(red,Qk*matrix_inverse(Q0));
   }
   def ygDk2(A,B,C,X,K) {
     return (1/(-B[K-1]+1))*ygDk(A,B,C,X,K);
   }
   
   def containNegative(L) {
     if (type(L) <= 1)  return(L<0);
     if (type(L) >=4) {
       for (I=0; I<length(L); I++) {
         if (containNegative(L[I])) return(1);
       }
       return(0);
     }
     error("Not an ordered object.");
   }
   /* 2014.08.09,  K=1,2,... */
   def check23(K) {
      A=-3; C=-1; B=[-3,-2,-3];
      /* A=-5; C=-6; B=[-1,-2,-3];   negative marginal */
      /* A=-5; C=-2; B=[-1,-2,-3];   */
     /* A=-5; C=2; B=[-1,-2,-3];   */
     /* A=-2; C=2; B=[-4]; */
      if (containNegative(abc2marginal(A,B,C))) error("Negative marginal");
     M = length(B);
     Alp = abc2alpha(A,B,C);
     F=hgpoly_fd(A,B,C);
     F=atofd(F[0],A,B,C,F[2]);
     V=F[1]; F=F[0];
     Fvec=newvect(M+1);
     Fvec[0] = F;
     for (I=1; I<=M; I++)  Fvec[I] =  ((V[I-1]-1)/Alp[I])*diff(F,V[I-1]);
     B2 = newvect(M,B);
     B2[K-1] = B[K-1]-1;
     B2 = vtol(B2);
     Alp2=abc2alpha(A,B2,C);
     F2=hgpoly_fd(A,B2,C);
     F2=atofd(F2[0],A,B2,C,F2[2])[0];
     printf("V=%a, B=%a, B2=%a\n",V,B,B2);
     printf("F=%a\n",F);
     printf("F2=%a\n",F2);
     /* Todo, how to modify the constant term? --> (-B[0]) in case of C>0*/
     Fvec2=newvect(M+1);
     Fvec2[0] = F2;
     for (I=1; I<=M; I++)  Fvec2[I] =  ((V[I-1]-1)/Alp2[I])*diff(F2,V[I-1]);
     /*  R=map(red,Fvec2-(1/(-B[K-1]+1))*ygDk(A,B,C,V,K)*Fvec); */
     R=map(red,Fvec2-ygDk2(A,B,C,V,K)*Fvec);
     return R;
   }
   
   /* 2014.08.10 */
   /* F(c-1) = Dc F(c)
    */
   def ygDc2(A,B,C,X) {
     M=length(B);
     Q0 = ygQmat2(A+1,B,C,X,0);
     Qm1 = ygQmat2(A+1,B,C,X,M+1);
     return map(red,(C-A-1)*Q0*matrix_inverse(Qm1));
   }
   def ygDc(A,B,C,X) {
     return ((1/(C-1))*ygDc2(A,B,C,X));
   }
   def check24() {
      A=-3; C=-1; B=[-3,-2,-3];
      /* A=-5; C=-6; B=[-1,-2,-3];   negative marginal */
      /* A=-5; C=-2; B=[-1,-2,-3];  */
      /* A=-5; C=2; B=[-1,-2,-3];   */
     /* A=-5; C=2; B=[-4];  */
      if (containNegative(abc2marginal(A,B,C))) error("Negative marginal");
     M = length(B);
     Alp = abc2alpha(A,B,C);
     F=hgpoly_fd(A,B,C);
     F=atofd(F[0],A,B,C,F[2]);
     V=F[1]; F=F[0];
     Fvec=newvect(M+1);
     Fvec[0] = F;
     for (I=1; I<=M; I++)  Fvec[I] =  ((V[I-1]-1)/Alp[I])*diff(F,V[I-1]);
     C2 = C-1;
     Alp2=abc2alpha(A,B,C2);
     F2=hgpoly_fd(A,B,C2);
     F2=atofd(F2[0],A,B,C2,F2[2])[0];
     printf("V=%a, C=%a, C2=%a\n",V,C,C2);
     printf("F=%a\n",F);
     printf("F2=%a\n",F2);
     /* Todo, how to modify the constant term? --> (-B[0]) in case of C>0*/
     Fvec2=newvect(M+1);
     Fvec2[0] = F2;
     for (I=1; I<=M; I++)  Fvec2[I] =  ((V[I-1]-1)/Alp2[I])*diff(F2,V[I-1]);
     R=map(red,Fvec2-ygDc2(A,B,C,V)*Fvec);
     return R;
   }
   
   /* F(a-1) = Da F(a)
    */
   def ygDa2(A,B,C,X) {
     M=length(B);
     Cmat = ygCmat2(A,B,C);
     Qm1 = ygQmat2(A,B,C,X,M+1);
     return map(red,-(1/(C-A))*Qm1*matrix_inverse(Cmat));
   }
   def ygDa(A,B,C,X) {
     Factor=-(A-1); /* Todo */
     return (Factor*ygDa2(A,B,C,X));
   }
   def check25() {
      A=-3; C=-1; B=[-3,-2,-3];
      /* A=-5; C=-6; B=[-1,-2,-3];   negative marginal */
     /* A=-5; C=-2; B=[-1,-2,-3]; */
      /* A=-5; C=2; B=[-1,-2,-3]; */
      /* A=-5; C=2; B=[-4];   */
      if (containNegative(abc2marginal(A,B,C))) error("Negative marginal");
     M = length(B);
     Alp = abc2alpha(A,B,C);
     F=hgpoly_fd(A,B,C);
     F=atofd(F[0],A,B,C,F[2]);
     V=F[1]; F=F[0];
     Fvec=newvect(M+1);
     Fvec[0] = F;
     for (I=1; I<=M; I++)  Fvec[I] =  ((V[I-1]-1)/Alp[I])*diff(F,V[I-1]);
     A2 = A-1;
     Alp2=abc2alpha(A2,B,C);
     F2=hgpoly_fd(A2,B,C);
     F2=atofd(F2[0],A2,B,C,F2[2])[0];
     printf("V=%a, A=%a, A2=%a\n",V,A,A2);
     printf("F=%a\n",F);
     printf("F2=%a\n",F2);
     Fvec2=newvect(M+1);
     Fvec2[0] = F2;
     for (I=1; I<=M; I++)  Fvec2[I] =  ((V[I-1]-1)/Alp2[I])*diff(F2,V[I-1]);
     R=map(red,Fvec2-ygDa2(A,B,C,V)*Fvec);
     return R;
   }
   
   /* ygDc(-4,[-1,-1],-2,[1/2,1/3]); --> det(Dm1)=(c-a)/(c-b1-b2)=0
       if b=[-2,..] then OK.
      a,[b1,b2],c1 --fctr--> c-a-1, c-b1-b2-1
   */
   
   /* 2014.08.19 */
   /* Return value: [F_D like poly, V]
     Note that the constant term is not 1 */
   def poly_fd(A,B,C) {
     F=hgpoly_fd(A,B,C);
     return(atofd(F[0],A,B,C,F[2]));
   }
   def poly_fd_beta(Beta) {
     F=hgpoly_fd_beta(length(Beta)-2,Beta);
     return(atofd_beta(F[0],Beta,F[2]));
   }
   
   def poly_fdvec(A,B,C) {
     M = length(B);
     Alp = abc2alpha(A,B,C);
     F=poly_fd(A,B,C);
     V=F[1]; F=F[0];
     Fvec=newvect(M+1);
     Fvec[0] = F;
     for (I=1; I<=M; I++)  Fvec[I] =  ((V[I-1]-1)/Alp[I])*diff(F,V[I-1]);
     return [Fvec,V];
   }
   def poly_fdvec_beta(Beta) {
     Marginal = beta2marginal(Beta);
     CSum=Marginal[1];
     RSum=Marginal[0];
     ABC=marginal2abc(RSum,CSum);
     return(poly_fdvec(ABC[0],ABC[1],ABC[2]));
   }
   
   def check26() {
      A=-3; C=-1; B=[-3,-2,-3];
      /* A=-5; C=-6; B=[-1,-2,-3];   negative marginal */
     /* A=-5; C=-2; B=[-1,-2,-3]; */
      /* A=-5; C=2; B=[-1,-2,-3]; */
      /* A=-5; C=2; B=[-4];   */
      if (containNegative(abc2marginal(A,B,C))) error("Negative marginal");
     Fvec = poly_fdvec(A,B,C);
     V=Fvec[1]; Fvec=Fvec[0];
     Fvec2 = poly_fdvec(A-1,B,C)[0];
     R=map(red,Fvec2-ygDa2(A,B,C,V)*Fvec);
     return R;
   }
   
   /* Dc2*Da2 c->c-1, a->a-1.
     C-A is canceled.  Not yet tested.
   */
   def ygDca2(A,B,C,X) {
     Q0 = ygQmat2(A,B,C,X,0);
     Cmat = ygCmat2(A,B,C);
     return map(red,-Q0*matrix_inverse(Cmat));
   }
   
   /*
     Table. B[0] is the first row sum.
     B[0] |
     CB   |
            B[1], B[2], ..., B[M+1]
     hg poly in the case of CB==1.
     2014.08.20
    */
   def hgpoly_fd_beta1(M,B) {
     CB = (base_sum(B,0,0,length(B)-1)-B[0])-B[0];
     if (CB != 1) error("The second row sum CB is not equal to 1.");
     /* Assume that CB == 1 */
     M = length(B)-2;
     N = (M+1)*2;
     Vx=[];
     for (I=0; I<N; I++) Vx = cons(util_v(x,[I+1]),Vx);
     Vx=reverse(Vx);
     F = 0;
     for (K=0; K<=M; K++) {
       Mon=1;
       for (J=0; J<=M; J++) {
         if (J != K) {
            Mon = Mon*(Vx[J]^(B[J+1]))/fac(B[J+1]);
         }
       }
       if (B[K+1] == 0) error("B[K+1]-1 is negative.");
       Mon = Mon*((Vx[K]^(B[K+1]-1))/fac(B[K+1]-1))*Vx[M+1+K];
       F += Mon;
     }
     return [F,dp_ptod(F,Vx),Vx];
   }
   def check27() {
      Beta = [7,1,5,2];
      Beta = [6,1,3,2,1]; /* A=-6, B=[-3,-2,-1],C=-4 */
      Beta = [8,2,3,2,2]; /* abc2marginal(-8,[-3,-2,-2],-5); [[8,1],[2,3,2,2]] */
   
      F1=hgpoly_fd_beta1(length(Beta)-2,Beta);
      F2=hgpoly_fd_beta0(length(Beta)-2,Beta);
      return F1[0]-F2[0];
   }
   def hgpoly_fd1(A,B,C) {
     Marginal = abc2marginal(A,B,C);
     if (containNegative(Marginal)) error("Negative marginal");
     M = length(B);
     Beta =cons(Marginal[0][0],Marginal[1]);
     return hgpoly_fd_beta1(M,Beta);
   }
   
   /* fvect_poly for C <= 0 */
   def fvect_poly2(A,B,C,X) {
     if (C > 0) error("C is postive.");
     M=length(B);
     if (type(X) == 4) X=newvect(length(X),X);
     Xrule = assoc(xvars(M),vtol(X));
     CB = C-base_sum(B,0,0,length(B)-1)-1;
     if (CB == 0) {
       F = poly_fdvec(A,B,C)[0];
       return(base_replace(F,Xrule));
     }
     if (CB < 0) error("c-sum(b)-1 is negative.");
     Steps = CB-1;
     printf("Computing Dmat(ca) for parameters B=%a,X=%a\n",B,X);
     Dmat = ygDca2(a,B,c,X);
     Umat = map(red,matrix_inverse(Dmat));
   
     F = poly_fdvec(A-Steps,B,C-Steps)[0];
     F = base_replace(F,Xrule);
     for (I=0; I<Steps; I++) {
       Up = base_replace(Umat,[[a,A-Steps+1+I],[c,C-Steps+1+I]]);
       F = Up*F;
     }
     return(F);
   }
   def check28() {
     A=-5; B=[-4,-3];C=-2; X=[1/2,1/3];
     A=-5; B=[-4,-3,-2];C=-2; X=[1/2,1/3,1/5];
     printf("marginal=%a\n",abc2marginal(A,B,C));
     F1=fvect_poly2(A,B,C,X);
     printf("F1=%a, computing F2 by a bruce force method.\n",F1);
     F2=base_replace(poly_fdvec(A,B,C)[0],assoc(xvars(length(X)),X));
     return(F1-F2);
   }
   def check28b(T) {
     Beta=[6,5,4,1,1,1]; /* ogawa data in check21, it takes a few seconds */
     X = [1/7,1/3,1/2,1/11];
     Ma = beta2marginal(Beta);
     ABC=marginal2abc(Ma[0],Ma[1]);
     A=ABC[0]; B=ABC[1]; C=ABC[2];
     printf("marginal=%a\n",abc2marginal(A,B,C));
     F1=fvect_poly2(A,B,C,X);
     if (T == 0) return(F1);
     printf("F1=%a\n",F1);
     F=poly_fdvec_beta(Beta)[0]$
     F2=base_replace(F,assoc(xvars(4),X));
     return(F1-F2);
   }
   
   /*
    works in case of row sum is zero or column sum is 0.
   */
   def ahvec_beta(Beta,Y) {
     Opt = getopt();
     Ma=beta2marginal(Beta);
     M=length(Beta)-2;
     RS=newvect(2,Ma[0]);
     CS=newvect(M+1,Ma[1]);
     Mnew=M;
     Fvec=newvect(M+1);
     for (I=0, J=0; I<=M; I++) {
       if (CS[I]==0) {
         Mnew--; Fvec[I] = -1;
       }else{ Fvec[I] = J; J++; }
     }
     if (Mnew < 0) error("zero matrix");
     Ynew=newmat(2,Mnew+1);
     CSnew=newvect(Mnew+1);
     for (I=0; I<=M; I++) {
       if (Fvec[I] >= 0) {
         J = Fvec[I];
         Ynew[0][J] = Y[0][I];
         Ynew[1][J] = Y[1][I];
         CSnew[J] = CS[I];
       }
     }
     if (RS[1] == 0) {
       T=RS[0]; RS[0] = RS[1]; RS[1] = T;
       for (J=0; J<=Mnew; J++) {
         T=Ynew[0][J];
         Ynew[0][J] = Ynew[1][J];
         Ynew[1][J] = T;
       }
     }
     ABCnew=marginal2abc(vtol(RS),vtol(CSnew));
     printf("RS=%a, CSnew=%a, Ynew=%a\n",RS,CSnew,Ynew);
     Tnew = ygahvec(ABCnew[0],ABCnew[1],ABCnew[2],Ynew | option_list=Opt);
     Fnew =Tnew[0];
     for (I=0; I<=M; I++) {
       if (Fvec[I] == -1) Fvec[I] = 0;
       else {
         J=Fvec[I];
         Fvec[I] = Fnew[J];
       }
     }
     return cons(Fvec,cdr(Tnew));
   }
   /*
   ahvec_beta([3,1,2,0,3],[[1,1/2,1/3,1/5],[1,1,1,1]]|all=1);
   
   [2402] abc2marginal(-5,[-4,-3],-2);
   [[5,4],[2,4,3]]
   [2403] marginal2abc([4,5],[2,4,3]);
   [-4,[-4,-3],-1]   c is again negative.
   
   2014.08.26
   ahvec_beta([10,1,2,3,4],[[3,1,1,1],[1,1/2,1/3,1/5]]|all=1);
   returns "devision by 0".  Reason?
   
   ahvec_abc(-6,[-2,-3],-4,[[1,1/2,1/3],[1,1,1]] | all=1);
    -> 2014.08.27 it works.
   cgi/r-fd.rr  r_ahvec(-6,[-2,-3],-4,[[1,1/2,1/3],[1,1,1]] | all=1);
   */
   
   def ahvec_abc(A,B,C,X) {
     Opt=getopt();
     Ma=abc2marginal(A,B,C);
     Beta=cons(Ma[0][0],Ma[1]);
     return ahvec_beta(Beta,X | option_list=Opt);
   }
   /*
     It returns [E[U_10], E[U_11], ..., E[U_1m]]
    */
   def expectation_abc(A,B,C,X) {
     R=ahvec_abc(A,B,C,X|all=1);
     Gamma=R[1];
     Der=R[0];
     /* printf("Der=%a\n",Der);  */
     Z=R[2]*Gamma;
     Der2 = newvect(length(Der));
     for (I=0; I<length(Der); I++) Der2[I] = Der[I]*Gamma;
     Der2 = vtol(Der2);
     E=newvect(length(Der));
     for (I=0; I<length(Der); I++) E[I] = X[1][I]*Der2[I]/Z;
     return(vtol(E));
 }  }
   
 #if  !USE_MODULE  #if  !USE_MODULE
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.12
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>