Return to dsolv.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

XM_debug is set to 0 for this package.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/dsolv.rr,v 1.2 2019/08/09 02:19:55 takayama Exp $ */
/* Old: dsolv */
/* dsolv: Documents are at dsolv.texi */
load("ok_diff.rr")$
load("ok_dmodule.rr")$

Dsolv_debug_dual_by_table=0$
Dsolv_debug_normal_vec=0$
Dsolv_debug_dual=0$
Dsolv_message_starting_term=1$
Dsolv_message_initial=0$

XM_debug=0$ print("XM_debug is set to 0 for this package.")$

def dsolv_start() {
  print("Starting ox_sm1 ...",0)$
  sm1.auto_reduce(1)$
  print("Done.")$
  print("All Grobner bases to be computed will be reduced Grobner bases.")$
}
print("Type in dsolv_start(); before using functions in this module.")$
def dsolv_consw(V,W) {
  N = length(V);
  Ans = [ ];
  for (I=0; I<N; I++) {
    Ans = append(Ans,[V[I],-W[I]]);
  }
  for (I=0; I<N; I++) {
    Ans = append(Ans,[strtov("d"+rtostr(V[I])),W[I]]);
  }
  return([Ans]);
}

def sm1_d(X) {
  return(strtov("d"+rtostr(X)));
}

def dsolv_initial(F1,V1,W1) {
  extern Dsolv_message_initial;
  S=[F1,V1,W1];
  F=S[0];
  V=S[1];  
  DV=map(sm1_d,V);
  W=S[2];
  N = length(V);
  G = sm1.gb([F,V,dsolv_consw(V,W)]);
  In = G[1]; 
  In = map(subst,In,h,1);
  In = sm1.gb([In,V])[0];  /* Computing the reduced basis. */
  if (Dsolv_message_initial) {
    print("Initial ideal is ",0)$ print(In)$
    print([In,V,V,DV,V]);
  }
  Ans = [ ];
  for (I=0; I<length(In); I++) {
    D = sm1.distraction([In[I],V,V,DV,V]);
    Ans = append(Ans,[D]);
  }
  return(Ans);
}

def dsolv_test_initial() {
  Mu = 2; Beta=1/2;
  F = sm1.appell1([Mu+Beta,Mu+1,Beta,Beta]);
  A = dsolv_initial(F[0],F[1],[1,3]);
  return(A);
}

def dsolv_normal_vec(F,G,Std,V) {
  Dp_f = dp_ptod(F,V);
  Dp_g =map(dp_ptod,G,V);
  Dp_std = map(dp_ptod,Std,V);
  return(dsolv_dp_normal_vec(Dp_f,Dp_g,Dp_std));
}

def dsolv_dp_normal_vec(F,G,Std) {
  extern Dsolv_debug_normal_vec;
  Index = [ ];
  for (I=0; I<length(G); I++) {
    Index = append(Index,[I]);
  }
  Vec_g = newvect(length(G),G);
  FF = dp_true_nf(Index,F,Vec_g,1);
  FN = FF[0]; FD = FF[1]; 
  A = newvect(length(Std));
  if (Dsolv_debug_normal_vec) print(FN);
  while (FN != 0) {
    H = dp_hm(FN);
    FN = FN - H;
    I = dsolv_where(dp_ht(H),Std);
    if (I < 0) error("dsolv_dp_normal_vec: index is out of bound for Std");
    A[I] = dp_hc(H)/FD;
  }
  return([FF[0]/FF[1],A]);
}

def dsolv_where(M,Std) {
  M2 = dp_ht(M);
  for (I=0; I<length(Std); I++) {
    if (M2 == Std[I]) return(I);
  }
  return(-1);
}

def dsolv_companion_matrix(F,V,K) {
  G = gr(F,V,0);
  Dp_g = map(dp_ptod,G,V);
  Dp_std = dp_mbase(Dp_g);
  Std = map(dp_dtop,Dp_std,V);
  Ans = [ ];
  for (I=0; I<length(Std); I++) {
    T = dsolv_normal_vec(Std[I]*V[K],G,Std,V);
    Vec = T[1];
    print(Vec);
    Ans = append(Ans, [ vtol(Vec) ]);
  }
  return(newmat(length(Ans),length(Ans),Ans));
}

def dsolv_test_companion_matrix() {
  return(dsolv_companion_matrix([x^2+y^2-4,x*y-1],[x,y],0));
}

def dsolv_dual_by_table(G) {
  extern Dsolv_debug_dual_by_table;
  Std = dp_mbase(G);
  N = size(dp_etov(dp_ht(G[0])))[0];  /* N is the number of variables. */
  L = [ ];
  NF = [ ];  /* NF[I] is the vector of coeffients of the normal form of L[I].*/
  I = [ ];   /* Finally, it is equal to G \cap Monomials */

  M = dp_vtoe(newvect(N));
  while(1) {
    if (Dsolv_debug_dual_by_table) {
      print("Reducing the monomial ",0);
      print(M);
    }
    if (dsolv_member(M,I)) {
       M = dsolv_inc(M);
       if (M == 0) return([L,NF]);
    }else{
       T = dsolv_dp_normal_vec(M,G,Std);
       NM = T[0];
       if (NM == 0) {
          I = append(I,[M]);
          M = dsolv_inc(M);
          if (M == 0) return([L,NF]);
       }else{
          L = append(L,[ M ]);
          NF = append(NF, [ T[1] ]);
          M = dsolv_inc_by_one(M);
       }
     }
   }
   return([L,NF]);
}

/* Increase the last exponent by one. */
def dsolv_inc_by_one(M) {
  MV = dp_etov(M);
  N = size(MV)[0];
  MV[N-1] = MV[N-1]+1;
  return(dp_vtoe(MV));
}

/* Get the next non-divisible monomial by the lex order
   when the given M belongs to the ideal. */
def dsolv_inc(M) {
  MV = dp_etov(M);
  N = size(MV)[0];
  for (I=N-1; I>=0; I--) {
    if (MV[I] != 0) {
      MV[I] = 0;
      if (I == 0) return(0);
      else {
         MV[I-1] = MV[I-1]+1;  
         return(dp_vtoe(MV));
      }
    }
  }
}

def dsolv_member(M,I) {
  if (I == [ ]) return(0);
  Index = newvect(length(I));
  for (J=0; J<length(I); J++) {Index[J] = J;}
  Index = vtol(Index);

  T = dp_nf(Index,M,newvect(length(I),I),0);
  if (T == 0) return(1);
  else return(0);
}

def dsolv_test_dual_by_table() {
  /* G = [x^3,x^2*y,y^3,x*y*z,z^2]; V=[x,y,z]; */
  G = [x1+x2+x3+x4+x5,x1+x2-x4,x2+x3-x4,x1*x3,x2*x4];
  V = [x1,x2,x3,x4,x5];
  G = gr(G,V,0);
  Dp_g = map(dp_ptod,G,V);
  return(dsolv_dual_by_table(Dp_g));
}

/* A is alpha and B is beta in the book [SST; page 73].
  it returns B!/A!
*/
def dsolv_weight(FA,FB) {
  A = dp_etov(FA);
  B = dp_etov(FB);
  N = size(A)[0];
  Ans = 1;
  for (I=0; I<N; I++) {
    Ans = Ans*(fac(B[I])/fac(A[I]));
  }
  return(Ans);
}

/* F must be primary to the maximal ideal */
def dsolv_dual(F,V) {
  extern Dsolv_debug_dual;
  G = gr(F,V,0);
  Dp_g = map(dp_ptod,G,V);
  Std = dp_mbase(Dp_g);
  T = dsolv_dual_by_table(Dp_g);
  if (Dsolv_debug_dual) {
    print("------------- dual table --------------------");
    print(T);
  }
  MI = T[0]; Table=T[1];
  Ans = newvect(size(Table[0])[0]);
  for (I=0; I<size(Ans)[0]; I++) {
    H = 0;
    for (J=0; J<length(MI); J++) {
      H = H+MI[J]*Table[J][I]*dsolv_weight(MI[J],Std[I]);
    }
    Ans[I] = H;
  }
  Ans = vtol(Ans);
  Ans = map(dp_dtop,Ans,V);
  return(Ans);
}

def dsolv_act(L,F,V) {
  LL = dmodule_d_op_fromasir([L],V)[0];
  A = diff_act(LL,F,V);
  A = red(A);
  return(A);
}

def dsolv_test_dual() {
  print("test1 -------------------------");
  F = [x^3,x^2*y,y^3,x*y*z,z^2]; V=[x,y,z]; 
  print(F);
  print(dsolv_dual(F,V));
  print("test2 ----------------------------");
  F = [x1+x2+x3+x4+x5,x1+x2-x4,x2+x3-x4,x1*x3,x2*x4];
  /* In this case, just replacement is fine. In general, use sm1.mul */
  Ftheta = map(subst,F,
             x1,x1*dx1, x2,x2*dx2, x3,x3*dx3, x4,x4*dx4, x5,x5*dx5);
  V = [x1,x2,x3,x4,x5];
  /* V = [x5,x4,x3,x2,x1];*/
  print(F);
  print("Grobner basis is -------------");
  print(gr(F,V,0));
  A=dsolv_dual(F,V);
  print(map(fctr,A));
  print("Are they solutions?");
  for (I=0; I<length(A); I++) {
    H = subst(A[I],
        x1,log(x1),x2,log(x2),x3,log(x3),x4,log(x4),x5,log(x5));
    print(H,0); print(" ===> ",0);
    print(map(dsolv_act,Ftheta,H,V));
  }
  /* Solutions at the page 75. */
  print(map(dsolv_act,Ftheta, log(x1)+log(x3)-log(x2)-log(x5),V));
  print(map(dsolv_act,Ftheta, log(x2)+log(x4)-2*log(x5),V));
  return(A);
}


def dsolv_solv1a(F) {
  V = var(F);
  if (deg(F,V) != 1) return([]);
  return([V,red(-coef(F,0)/coef(F,1))]);
}


def dsolv_solv_linear(LL,V) {
  L = gr(LL,V,2);
  N = length(L);
  Ans = newvect(length(V));
  for(J=0; J<length(V); J++) {
     Ans[J] = "?";
  }
  for (I=0; I<N; I++) {
    S = dsolv_solv1a(L[I]);
    if (S == []) return([]);
    for (J=0; J<length(V); J++) {
      if (V[J] == S[0]) {
        Ans[J] = S[1];
      }
    }
  }
  return(Ans);
}

def dsolv_subst(H,V,P) {
  A = H;
  N = length(V);
  for (I=0; I<N; I++) {
    A = subst(A,V[I],V[I]+P[I]);
  }
  return(A);
}
def dsolv_subst_log(H,V) {
  A = H;
  N = length(V);
  for (I=0; I<N; I++) {
    A = subst(A,V[I],log(V[I]));
  }
  return(A);
}
def dsolv_starting_term(F,V,W) {
  extern Dsolv_message_starting_term;
  Ans = [ ];
  Exp = [ ];
  if (Dsolv_message_starting_term) {
    print("Computing the initial ideal.");
  }
  S = dsolv_initial(F,V,W);
  if (Dsolv_message_starting_term) {
    print("Done.");
    print("Computing a primary ideal decomposition.");
  }
  S = primadec(S,V);
  if (Dsolv_message_starting_term) {
    print("Primary ideal decomposition of the initial Frobenius ideal to the direction ",0);
    print(W,0);
    print(" is ");
    print(S);
    print(" ");
  }
  N = length(S);
  for (I=0; I<N; I++) {
    /* Check if S[I][1] is linear or not. If it is not, output an error.
       Not yet implemented. */
    P = dsolv_solv_linear(S[I][1],V);
    Exp = append(Exp, [ P ]);
    if (Dsolv_message_starting_term) {
       print("----------- root is ",0); print(P);
    }
    Q = S[I][0];  /* Primary component associated at X=P */
    Q = map(dsolv_subst,Q,V,P);
    Dual = dsolv_dual(Q,V);
    if (Dsolv_message_starting_term) {
       print("----------- dual system is ",0); print(Dual);
    }
    Dual = map(dsolv_subst_log,Dual,V);

    /* 
      P = dp_dtop(dp_vtoe(P),V); buggy for rational exponents.
    */
    P1 = 1;
    for (I1 = 0; I1 <length(V); I1++) {
      P1 = P1 * V[I1]^P[I1];
    }
    P = P1;

    Ans = append(Ans, [ vtol(P*newvect(length(Dual),Dual)) ]);
  }
  print("  ");
  return([Exp,Ans]);
}

def dsolv_test_starting_term() {
  /* Example 2.6.4. p.99 of [SST]. */
  F = sm1.gkz( [ [[1,1,1,1,1],[1,1,0,-1,0],[0,1,1,-1,0]], [1,0,0]]);
  return(dsolv_starting_term(F[0],F[1],[1,1,1,1,0]));
}

def dsolv_ttest_starting_term() {
  /* Example 2.6.4. p.99 of [SST] with beta = [-1,0,0]  */
  F = sm1.gkz( [ [[1,1,1,1,1],[1,1,0,-1,0],[0,1,1,-1,0]], [-1,0,0]]);
  A = dsolv_starting_term(F[0],F[1],[1,1,1,1,0]);
  A = A[1][0];
  print(A);
  B=map(subst,A,x1,x,x2,x,x3,x,x4,x,x5,1); /* Do not use t */
  print("Drawing a graph of ");
  print(B);
  gnuplot.plot_function(B);
  return(0);
}

def dsolv_test_starting_term_1() {
  F = sm1.gkz( [ [[2,1,0,0,1,2,1],[1,2,2,1,0,0,1],[0,0,1,2,2,1,1]], [1,1,1]]);
  W = [1,1,1,1,1,1,0];  /* We need to get a reduced GB. 
                           This part has not been implemented.==> Done. 2/26.*/
  A = dsolv_starting_term(F[0],F[1],W);
  A = A[1][0];
  print(A);
  B=map(subst,A,x1,x,x2,x,x3,x,x4,x,x5,x,x6,x,x7,1); /* Do not use t */
  print("Drawing a graph of ");
  print(B);
  gnuplot.plot_function(B);
  return(0);
}

/* Functions to consruct higher order terms. */
def dsolv_euler_diff(F,G,VD,VV) {
  H=getopt(usage);
  if (type(H) != -1) {
    /* help message for debugging. */
    print("dsolv_euler_diff(diff op in terms of euler op, poly, euler op, variables).");
    print("Example 1: dsolv_euler_diff(s1*s2^2+y1*s1^2,y1^2*y2+y2^3,[s1,s2],[y1,y2])");
    return 0;
  }
  
  FF = dp_ptod(F,VD);
  GG = dp_ptod(G,VV);
  Ans = 0;
  while (FF != 0) {
    HT = dp_ht(FF); HC = dp_hc(FF);
    FF = dp_rest(FF);
    HT = dp_etov(HT);
    TT = dsolv_euler_diff_1(HT,GG);
    Ans += dp_dtop(HC*TT,VV);
  }
  return Ans;
}

def dsolv_pfac_n(HE,E) {
  if (type(HE) == 4) {
    N = length(HE);
  }else{
    N = size(HE)[0];
  }
  A = 1;
  for (I=0; I<N; I++) {
    for (J=0; J<E[I]; J++) {
      A *= HE[I]-J;
    }
  }
  return A;
}

def dsolv_euler_diff_1(E,F) {
  /* E is a vector (degree of euler operators) */
  /* F is a distributive polynomial. */
  Ans = 0;
  while (F != 0) {
    HE = dp_etov(F);
    Ans += dsolv_pfac_n(HE,E)*dp_hm(F);
    F = dp_rest(F);
  }
  return Ans;
}


def dsolv_d_mult(P,Q,DV,VV) {
  H = getopt(usage);
  if (type(H) != -1) {
    print("dsolv_d_mult(P,Q, d-variables, x-variables)");
    print("Example: dsolv_d_mult(x*dx^2,x^2,[dx],[x])");
    return 0;
  }
  N = length(DV);
  Kmax = newvect(N);
  for (I=0; I<N; I++) {
    KK = deg(P,DV[I]);
    KK2 = deg(Q,VV[I]);
    if (KK > KK2) {
      KK = KK2;
    }
    Kmax[I] = KK+1;
  }
  K = newvect(N);
  Ans = 0;
  while (K[0] < Kmax[0]) {
    PP = P;
    QQ = Q;
    for (I=0; I<N; I++) {
       PP = dsolv_diff(PP,DV[I],K[I]);
       QQ = dsolv_diff(QQ,VV[I],K[I]);
    }
    Ans += PP*QQ*(1/dsolv_pfac_n(K,K));
    dsolv_increase_vect(K,Kmax,N);
  }  
  return Ans;
}

def dsolv_diff(F,X,M) {
  for (I=0; I<M; I++) {   
     F = diff(F,X);
  }
  return F;
}
def dsolv_increase_vect(K,Kmax,N) {
  for (I=N-1; I>=0; I--) {
     K[I]++;
     if (K[I] >= Kmax[I] && I != 0) {
       K[I] = 0;
       dsolv_increase_vect(K,Kmax,I);
       return ;
     }else{
       return;
     }
  }
}

def dsolv_homogeneous_base(W,K,V,B) {
  H = getopt(usage);
  if (type(H) != -1) {
    print("dsolv_homogeneous_base(weight W, degree K, variables V, bases B)");
    print("It returns a set of exponents E such that W-degree of B^E is equal");
    print("to K.");
    print("Example:  dsolv_homogeneous_base([1,1,1,1,0],6,[x1,x2,x3,x4,x5],[x1*x3/(x2*x5),(x2*x4)/x2^5])");
  }

  N = length(V);
  M = legnth(B);
  WB = newvect(M);
  for (I=0; I<M; I++) {
    WB[I] = 0;
    F = nm(B[I]);
    for (J=0; J<N; J++) {
      WB[I] += deg(F,V[J])*W[J];
    }
    F = dn(B[I]);
    for (J=0; J<N; J++) {
      WB[I] -= deg(F,V[J])*W[J];
    }
  }
  NOT_YET ;  

}

end$