Return to taka_ahg.rr CVS log

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

Fixed a bug of grassman(K,N).

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_ahg.rr,v 1.6 2004/07/04 00:05:00 takayama Exp $ */

module taka_ahg;

/* ------------ list of local functions ---------- */
localf b$
localf bb$
localf toPrimitive$
localf toPrimitive2$
localf gcd_L$
localf a_mc, a_dprism $
localf a_grassman $
localf a_sst0, a_sst26, a_sst40, a_sst99, a_sst107, a_sst125$
localf a_sst136, a_sst141, a_sst146, a_sst156, a_sst163, a_sst166$
localf a_sst173, a_sst177, a_sst188$ 

def gcd_L(L) {
  if (length(L) == 0) return(0);
  return igcd(L[0],gcd_L(cdr(L)));
}

def toPrimitive(A,F) {
/*  A set of points, F is a set of facets. */
  Nfacets = length(F);
  NewF = [ ];
  V=[]; Sign = 1;
  for (I=0; I<Nfacets; I++) {
    for (J=0; J<length(A); J++) {
      Ip = matrix_inner_product(A[J],F[I]);
      if ((Sign == -1) && (Ip > 0)) error("isPrimitive: broken facet data.");
      if (Ip < 0) { Sign = -1; }
      Ip = Sign*Ip;
      if (Ip != 0) V = cons(Ip,V);
    }
    if (length(V) == 0) G = 1; else G = gcd_L(V);
    if (G != 1) {
      FF=vtol(newvect(length(F[I]),F[I])/G);
    }else FF=F[I];
    NewF = cons(FF,NewF);
  }
  return reverse(NewF);
}

def toPrimitive2(A,F) {
/*  A : a point, F is a set of facets */
  Nfacets = length(F);
  V=[]; Sign = 1;
  for (I=0; I<Nfacets; I++) {
    Ip = matrix_inner_product(A,F[I]);
    if ((Sign == -1) && (Ip > 0)) error("isPrimitive2: broken facet data.");
    if (Ip < 0) { Sign = -1; }
    Ip = Sign*Ip;
    if (Ip != 0) V = cons(Ip,V);
  }
  NewF = [ ];
  if (length(V) == 0) G=1; else G = gcd_L(V);
  if (Sign < 0) G = -G;
  if (G != 1) {
    for (I=0; I<Nfacets; I++) {
      FF=vtol(newvect(length(F[I]),F[I])/G);
      NewF = cons(FF,NewF);
    }
  }else return(F);
  return reverse(NewF);
}

def b(A,Idx,V) {
  F = oxshell.facets(A);
  if (A != F[0]) error("b: points are sorted. Not implemented in this case.");
  F = F[1];

  P = A[Idx];

  Nfacets = length(F);
   
  /* to primitive */
  F = toPrimitive2(P,F);
  Bf = 1;
  for (I=0; I<Nfacets; I++) {
    H = matrix_inner_product(P,F[I]);
    if (H != 0) {
      B = matrix_inner_product(F[I],V);
      for (J=0; J<H; J++) {
        Bf = Bf*(B-J);  /* See Paper3/ip/ip2/hg.dvi [SST; compositio] */
      }
    }
  } 
  return(Bf);
}
/* bb is almost idential with b. A is transposed. */
def bb(A,Idx,V) {
  A = matrix_transpose(A);
  A = matrix_matrix_to_list(A);

  F = oxshell.facets(A);
  if (A != F[0]) error("bb: points are sorted. Not implemented in this case.");
  F = F[1];

  P = A[Idx];

  Nfacets = length(F);
   
  /* to primitive */
  F = toPrimitive2(P,F);
  Bf = 1;
  for (I=0; I<Nfacets; I++) {
    H = matrix_inner_product(P,F[I]);
    if (H != 0) {
      B = matrix_inner_product(F[I],V);
      for (J=0; J<H; J++) {
        Bf = Bf*(B-J);  /* See Paper3/ip/ip2/hg.dvi [SST; compositio] */
      }
    }
  } 
  return(Bf);
}

/* matrix A associated to a monomial curve.  Irregular.
   [[[1,2,3]],[b1]]
   Example. sm1.gkz(taka_ahg.a_mc(10));
 */
def a_mc(N) {
  A=[];
  for (I=N; I>=1; I--) {
    A=cons(I,A);
  }
  return [[A],[b1]];
}

/* matrix A associated to a distorted prism.  Irregular.
   cf. a_grassman(2,N)
 */
def a_dprism(N) {
  A=[];
  P = newvect(N+N+1);
  for (I=0; I<N; I++) P[I]=0;
  for (I=N; I<2*N+1; I++) P[I]=1;
  A = cons(vtol(P),A);

  for (J=0; J<N; J++) {
    P = newvect(N+N+1);
    for (I=0; I<N; I++) {
      if (I == J) P[I] = 1;
    }
    for (I=N; I<2*N+1; I++) {
      if (I-N-1 == J) P[I] = 1;
    }
    A = cons(vtol(P),A);
  }
  A = reverse(A);

  B = [];
  for (I=1; I<=N+1; I++) {
    B = cons(eval_str("b"+rtostr(I)),B);
  }
  B = reverse(B);
  return [A,B];
}

/* matrix A associated to the grassman E(k,n). Regular
   cf. E(2,N) <--> prism,  F_D
 */
def a_grassman(K,N) {
  A=[];
  M = N-K;
  for (J=0; J<M; J++) {
    P = newvect(K*M);
    for (I=J*K; I< (J+1)*K; I++) P[I] = 1;
    A = cons(vtol(P),A);
  }

  for (J=0; J<K; J++) {
    P = newvect(N+N+1);
    P = newvect(K*M);
    for (I=J; I< K*M; I = I+K) P[I] = 1;
    A = cons(vtol(P),A);
  }
  A = reverse(A);
  A = cdr(A);

  B = [];
  for (I=1; I<=K+M-1; I++) {
    B = cons(eval_str("b"+rtostr(I)),B);
  }
  B = reverse(B);
  return [A,B];
}

/* A that discussed from page 40 of [SST] (Saito-Sturmfels-Takayama) */
def a_sst40(M) {
  return a_grassman(2,M);
}

def a_sst0() {
  A = [[1,1,1],
       [0,1,2]];
  B=[b1,b2];
  return [A,B];
}

def a_sst26() {
  A = [[1,0,0,-1],
       [0,1,0,1],
       [0,0,1,1]];
  B=[b1,b2,b3];
  return [A,B];
}

def a_sst99() {
  A = [[1,1,1,1,1],
       [1,1,0,-1,0],
       [0,1,1,-1,0]];
  B = [1,0,0];
  return [A,B];
}

def a_sst107() {
  A = [[3,2,1,0],
       [0,1,2,3]];
  B = [b1,b2];
  return [A,B];
}

def a_sst125() {
  A = [[4,3,2,1,0],
       [0,1,2,3,4]];
  B = [b1,b2];
  return [A,B];
}

def a_sst136() {
  A = [[1,1,1,1,1],
       [-1,1,1,-1,0],
       [-1,-1,1,1,0]];
  B = [1,0,0];
  return [A,B];
}

def a_sst141() {
  A = [[1,1,1,1,1,1,1,1,1],
       [0,1,2,0,1,2,0,1,2],
       [0,0,0,1,1,1,2,2,2]];
  B = [b1,b2,b3];
  return [A,B];
}

def a_sst146() {
  A = [[1,1,1],
       [0,1,2]];
  B = [b1,b2];
  return [A,B];
}


def a_sst156() {
  A = [[1,1,1,1,1],
       [0,8,16,21,26]];
  B = [b1,b2];
  return [A,B];
}

def a_sst163() {
  A = [[1,1,1,1,1],
       [0,2,4,7,9]];
  B = [b1,b2];
  return [A,B];
}

def a_sst166() {
  A = [[1,1,1,1,1,1],
       [0,1,1,0,-1,-1],
       [-1,-1,0,1,1,0]];
  B = [1,0,0];
  return [A,B];
}

def a_sst173() {
  A = [[1,1,1,1],
       [0,1,3,4]];
  B = [1,2];
  return [A,B];
}

def a_sst177() {
  A = [[1,1,1,1,1],
       [0,1,2,3,4]];
  B = [b1,b2];
  return [A,B];
}


def a_sst188() {
  A = [[1,1,1,1,1,1],
       [0,1,1,0,-1,-1],
       [-1,-1,0,1,1,0]];
  B = [b1,b2,b3];
  return [A,B];
}

/* sst 222, 223, 224, 225. They have not yet input. */

endmodule$

Loaded_taka_ahg=1$

def test() {
 A=[[1,0,0],[1,1,0],[1,0,1],[1,1,1],[1,2,0]];
 print("A="); print(A); 
 print("P=",0); print(0);
 print(fctr(taka_ahg.b(A,0,[s1,s2,s3])));

 print("P=",0); print(3);
 print(fctr(taka_ahg.b(A,3,[s1,s2,s3])));

 print("P=",0); print(4);
 print(fctr(taka_ahg.b(A,4,[s1,s2,s3])));

}

end$