Return to pfphom.rr CVS log

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

Fixed wrong node dependencies.

/* -*- mode: C -*- */
/* $OpenXM: OpenXM/src/asir-contrib/packages/src/pfphom.rr,v 1.5 2003/07/27 13:18:48 takayama Exp $ */

/* cycles
     C1 = [t,1] \otimes q_1
     C2 = [0,t] \otimes q_1
     C3 = [0,t] \otimes q_2
     intersection := <lf(C),reg(dual(C))>
*/

/*&usage 
begin: pfphom_intersection(P)
 intersection matrix of homology cycles.
description:  
   Computing intersection matrix of cycles associated to p F_(p-1). 
   As to the meaning of parameters c1, c2, c3, ..., see the paper
   Ohara, Kyushu J. Math. Vol. 51 PP.123.

example:
         SS = pfphom_intersection(3)$
example_description:
   You get the intersection matrix of homologies for 3 F 2.


algorithm: Ohara, Sugiki, Takayama, Quadratic Relations for Generalized
     Hypergeometric Functions p F p-1

author: K.Ohara
end:
*/


extern Pfphom_RuleDual$
extern Pfphom_RuleNormal$
extern Pfphom_Parameters$

Pfphom_Parameters = [c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20]$

Pfphom_RuleDual  = [[c1, 1/c1], [c2, 1/c2], [c3, 1/c3], [c4, 1/c4], [c5, 1/c5], [c6, 1/c6], [c7, 1/c7], [c8, 1/c8], [c9, 1/c9], [c10,1/c10], [c11,1/c11], [c12,1/c12], [c13, 1/c13], [c14, 1/c14], [c16, 1/c16], [c17, 1/c17]]$

Pfphom_RuleNormal = [[c1,1/c2],[c3,1],[c6,1/(c5*c4)],[c9,1/(c8*c7)],[c12,1/(c10*c11)],[c15,1/(c13*c14)],[c18,1/(c16*c17)],[d2,c2-1],[d4,c4-1],[d5,c5-1],[d8,c8-1]]$

def pfphom_dual(F) {
    return base_replace(F, Pfphom_RuleDual);
}

/* Monodromy Functions for pFp-1 (P<=6) */
def pfphom_monodromy_triple(P) {
    C = Pfphom_Parameters;
    L2 = pfphom_monodromy_pair_kyushu(P);
    L3 = [C[3*P+1]*L2[0],L2[1],C[3*P+2]];
    return base_replace(L3,Pfphom_RuleNormal);
}


/*&usage 
begin: pfphom_monodromy_pair_kyushu(P)
description:  
It returns the pair of monodromy matrices.

example:
	     MP = pfphom_monodromy_pair_kyushu(3)$
example_description:
 You get a pair of monodromy matricies for 3F2 standing for
 two paths encircling 0 and 1.

algorithm: Ohara, Kyushu J. Math. Vol.51 PP.123 (1997)

end:
*/


/* Kyushu J. Math. Vol.51 PP.123 (1997) */
def pfphom_monodromy_pair_kyushu(P) {
    C = Pfphom_Parameters;
    M1 = pfphom_monodromy_M1(P);
    M2 = pfphom_monodromy_M2(P);
    MPair = [M1,M2];
    return base_replace(MPair,Pfphom_RuleNormal);
}

def pfphom_monodromy_M1(P) {
    C = Pfphom_Parameters;
    M1 = newmat(P,P);
    for(I=0;I<P;I++) {
        M1[I][I] = C[3*(P-I)];
        for(J=0;J<I;J++) {
            M1[I][J] = (-1)^(I+J+1) *(C[3*(P-I)+3]*C[3*(P-I)+2] - C[3*(P-I)]);
        }
    }
    return M1;
}

def pfphom_monodromy_M2(P) {
    C = Pfphom_Parameters;
    M2 = matrix_identity_matrix(P);
    CC = 1;
    for(I=0;I<P;I++) {
        CC = CC*C[3*I+2];
        M2[0][P-1-I] += (-1)^(P-1-I) * (CC-1);
    }
    return M2;  
}

/* Check Functions */
/* trans(M^-1)* S *Dual(M^-1) == S */
def pfphom_check(M,S) {
    R = (matrix_transpose(matrix_inverse(M)) * S) - (S * pfphom_dual(M));
    return base_cancel(R);
}

def pfphom_checkNumeric(M,S,Rule) {
    D = base_replace(pfphom_dual(M), Rule);
    M = base_replace(M, Rule);
    S = base_replace(S, Rule);
    R = (matrix_transpose(matrix_inverse(M)) * S) - (S * D);
    return base_cancel(R);
}

def pfphom_intersection_checkNumeric(ML,S) {
    I=0;
    /* fixed rule */
    RuleNumeric = [[c2, 1/3],[c4,-1/5],[c5,1/7],[c7,-1/11],[c8,1/13],[c10,-1/17],[c11,1/19],[c13,-1/23],[c14,1/29],[c16,-1/31]];
    print([I,"check(M1)=",pfphom_checkNumeric(ML[0],S,RuleNumeric)]);
    print([I,"check(M2)=",pfphom_checkNumeric(ML[1],S,RuleNumeric)]);
    /* random rule */
    for(I=1;I<=5;I++) {
        RuleNumeric = pfphom_makeRuleNumeric(I);
        print([I,"check(M1)=",pfphom_checkNumeric(ML[0],S,RuleNumeric)]);
        print([I,"check(M2)=",pfphom_checkNumeric(ML[1],S,RuleNumeric)]);
    }
}

def pfphom_makeRuleNumeric(Seed) {
    random(Seed);
    Length = length(Pfphom_Parameters);
    R = newmat(Length,2);
    for(I=0;I<Length;I++) {
        R[I][0]=Pfphom_Parameters[I];
        R[I][1]=random();
    }
    return matrix_matrix_to_list(R);
}

/* Matrix Functions */
def pfphom_matrix_to_identity(M) {

    return matrix_identity_matrix(car(size(M)));
}

def pfphom_di(M) {
    return base_cancel(matrix_inverse(M - pfphom_matrix_to_identity(M)));
}

def pfphom_getN1(M, L) {
    return pfphom_di(M) + (L/(L-1)) * pfphom_matrix_to_identity(M);
}

def pfphom_getN2(M, L) {
    return (L*M-1)/((L-1)*(M-1));
}

/* ML = [M1,M2,Mt]: List of monodormies, 
   M1,M2: monodromy matricies, Mt: a scalar, S: bilinear form */
def pfphom_intersection_ML(ML,S) {
    PP = car(size(S));
    SS = newmat(PP+1,PP+1);
    M1 = ML[0];
    M2 = ML[1];
    Mt = ML[2];

    if (type(Mt) == 6) {
        Mt = Mt[0][0];
    }
    if (type(M2) == 6) {
        M2 = M2[0][0];
    }

    SS[0][0] = red(S[0][0]*pfphom_dual(pfphom_getN2(M2, Mt)));

    CC1 = pfphom_dual(-1/(Mt-1));
    CC2 = pfphom_dual(-Mt/(Mt-1));
    for(I=0;I<PP;I++) {
        SS[0][I+1] = red(CC1 * S[0][I]);
        SS[I+1][0] = red(CC2 * S[I][0]);
    }

    SN1  = S * pfphom_dual(pfphom_getN1(M1, Mt));
    for(I=0;I<PP;I++) {
        for(J=0;J<PP;J++) {
            SS[I+1][J+1] = red(SN1[I][J]);
        }
    }
    return base_replace(SS,Pfphom_RuleNormal);
}

def pfphom_intersection(P) {
    if (P==1) {
        return matrix_identity_matrix(1);
    }
    return pfphom_intersection_ML(pfphom_monodromy_triple(P-1),pfphom_intersection(P-1));
}

/* Remark Ih for P is a PxP matrix,
   F == 1: Numerical check by random numbers */
def pfphom_intersection_check(Ih, F) {
    P = car(size(P));
    ML = pfphom_monodromy_triple(P);
    if (F == 1) {
        intersection_check_numeric(ML,Ih);
    }else {
        print(["check(M1)=",pfphom_check(ML[0],S)]);
        print(["check(M2)=",pfphom_check(ML[1],S)]);
    }
}

/**** Examples ****/

/* 1F0 */
/* monodromy$
XM1 = c4*matrix_list_to_matrix([[c3]])$
XM2 = matrix_list_to_matrix([[c2]])$
XMt = c5$
XM1 = base_replace(XM1,Pfphom_RuleNormal)$
XM2 = base_replace(XM2,Pfphom_RuleNormal)$
XMt = base_replace(XMt,Pfphom_RuleNormal)$
*/
/*
ML = pfphom_monodromy_triple(1)$
XM1 = ML[0]$ XM2 = ML[1]$ XMt = ML[2]$
Ih21 = pfphom_intersection(2);
*/
/* Ih21$
Ih21 = -matrix_list_to_matrix(
    [[(c2*c5-1)/(d5*d2), -c5/d5],
     [-1/d5, (c4*c5-1)/(d4*d5)]])$
Ih21 = base_replace(Ih21,Pfphom_RuleNormal)$
*/

/* 2F1 */
/* monodromy$
M1 = c7*matrix_list_to_matrix([[c6,0],[c5*c6-1,1]])$
M2 = matrix_list_to_matrix([[c2*c5,1-c2],[0,1]])$
Mt = c8$

M1 = base_replace(M1,Pfphom_RuleNormal)$
M2 = base_replace(M2,Pfphom_RuleNormal)$
Mt = base_replace(Mt,Pfphom_RuleNormal)$
*/
/*
ML = pfphom_monodromy_triple(2)$
M1 = ML[0]$ M2 = ML[1]$ Mt = ML[2]$
Ih32 = pfphom_intersection(3)$
*/
/* 3F2 */
/* monodromy$
MM1 = c10*matrix_list_to_matrix(
    [[c9,0,0],
     [c8*c9-c6,c6,0],
     [1-c5*c6,c5*c6-1,1]])$
MM2 = matrix_list_to_matrix(
    [[c2*c5*c8,1-c2*c5,c2-1],
    [0,1,0],
    [0,0,1]])$
MMt = c11$

MM1 = base_replace(MM1,Pfphom_RuleNormal)$
MM2 = base_replace(MM2,Pfphom_RuleNormal)$
MMt = base_replace(MMt,Pfphom_RuleNormal)$
*/
/*
ML = pfphom_monodromy_triple(3)$
MM1 = ML[0]$ MM2 = ML[1]$ MMt = ML[2]$

Ih43 = pfphom_intersection(4)$
print_xdvi_form(poly_factor(Ih43))$

ML = pfphom_monodromy_triple(3)$
MMM1 = ML[0]$ MMM2 = ML[1]$ MMMt = ML[2]$
print(["check(MMM1,Ih43)=",pfphom_checkNumeric(MMM1,Ih43,RuleNumeric1)])$
print(["check(MMM2,Ih43)=",pfphom_checkNumeric(MMM2,Ih43,RuleNumeric1)])$
*/


/* 4F3 */
/* monodromy$
MMM1 = c13*matrix_list_to_matrix(
    [[c12,        0,        0,       0],
     [c11*c12-c9, c9,       0,       0],
     [-c8*c9+c6,  c8*c9-c6, c6,      0],
     [c5*c6-1,    -c5*c6+1, c5*c6-1, 1]])$
MMM2 = matrix_list_to_matrix(
    [[c2*c5*c8*c11, 1-c2*c5*c8, c2*c5-1, 1-c2],
     [0,            1,          0,       0],
     [0,            0,          1,       0],
     [0,            0,          0,       1]])$
MMMt = c14$

MMM1 = base_replace(MMM1,Pfphom_RuleNormal)$
MMM2 = base_replace(MMM2,Pfphom_RuleNormal)$
MMMt = base_replace(MMMt,Pfphom_RuleNormal)$
*/
/*
ML = pfphom_monodromy_triple(3)$
MMM1 = ML[0]$ MMM2 = ML[1]$ MMMt = ML[2]$
Ih54 = pfphom_intersection(5)$
*/


/* Ih32 == Ti0 */
/* Takayama's 3F2 intersection matrix */
/*
Ti0 = base_cancel(matrix_list_to_matrix(
[[ (c8*c5*c2-1)/(((c8-1)*c5-c8+1)*c2+(-c8+1)*c5+c8-1),(-c8*c5*c2+c8)/(((c8-1)*c5-c8+1)*c2+(-c8+1)*c5+c8-1), (c8*c5)/((c8-1)*c5-c8+1) ],
[ (-c5*c2+1)/(((c8-1)*c5-c8+1)*c2+(-c8+1)*c5+c8-1), (((-c7+1)*c4*c5^2+((-c7*c8+c7)*c4+c7^2*c8-c7)*c5)*c2+((c7*c8-1)*c4-c7*c8+c7)*c5+(-c7^2+c7)*c8)/((((-c7+1)*c8+c7-1)*c4*c5^2+(((c7-1)*c8-c7+1)*c4+(c7^2-c7)*c8-c7^2+c7)*c5+(-c7^2+c7)*c8+c7^2-c7)*c2+((c7-1)*c8-c7+1)*c4*c5^2+(((-c7+1)*c8+c7-1)*c4+(-c7^2+c7)*c8+c7^2-c7)*c5+(c7^2-c7)*c8-c7^2+c7), ((c7*c8-1)*c5)/(((-c7+1)*c8+c7-1)*c5+(c7-1)*c8-c7+1) ],
[ (1)/((c8-1)*c5-c8+1), (c7*c8-1)/(((-c7+1)*c8+c7-1)*c5+(c7-1)*c8-c7+1), ((c7*c8-1)*c4*c5-c7*c8+1)/((((c7-1)*c8-c7+1)*c4+(-c7+1)*c8+c7-1)*c5+((-c7+1)*c8+c7-1)*c4+(c7-1)*c8-c7+1) ]]))$
*/


/* ---------------------------------------------------- pfp functions ---*/
/* taka_pfp from misc/200205/runge-kuttan.rr */
/* (th{x}+cc[1])(th{x}+cc[2])(th{x}+cc[3])(th{x}+cc[4])(th{x}+cc[5])  
     - x (th{x}+aa[1])(th{x}+aa[2])(th{x}+aa[3])(th{x}+aa[4])(th{x}+aa[5]) */
/* 
Taka_pfp_rule = [[cc[1],0],[cc[2],c1-1],[cc[3],c2-1],[cc[4],c3-1],[cc[5],c4-1],
                 [aa[1],a1],[aa[2],a2],[aa[3],a3],[aa[4],a4],[aa[5],a5]]$
*/
def taka_pfp_rule(CC,Size) {
  Rule = [[taka_pfp_iv(CC,1),0]];
  for (I=2; I<=Size; I++) {
    Rule = cons([taka_pfp_iv(CC,I),taka_pfp_iv(CC,I-1)-1],Rule);
  }
  return Rule;
}

def taka_pfp_iv(V,I) {
  return strtov(rtostr(V)+rtostr(I));
}

/*&usage begin: pfp_omega(P)
 It returns the Gauss-Manin connection Omega 
 for the generalized hypergeometric function 
  P F P-1 (aa1,aa2, ...; cc1, cc2, ...;x) .
 description:
 Define a vector valued function Y of which elements are 
 generalized hypergeometric function f_1=F and
 f_2=xdf_1/dx, f3=xd f_2/dx, ...
 It satisfies  dY/dx= Omega Y.
 Generalized hypergeometric function is defined by the series
 p F p-1(aa1,aa2, ...; cc1, cc2, ...;x) =
 sum(k=0,infty; (aa1)_k (aa2)_k .../( (1)_k (cc1)_k (cc2)_k ... ) x^k)
 example: pfp_omega(3);
end: */
def pfp_omega(P) { return taka_pfp_omega(P); }
/* Differential equation for p F p-1 is  dY/dx=omega[] Y       */
def taka_pfp_omega(Size) {
     Aaa = taka_pfp_fun(aa,Size); Ccc = taka_pfp_fun(cc,Size); 
     Tmp = newvect(Size);
     for (K=0; K<Size-1; K++) Tmp[K] = taka_pfp_nvec(K+1,Size);
     Tmp[Size-1] = taka_pfp_nnvec(Aaa,Ccc,Size);
     Tmp = map(vtol,Tmp);
     Tmp = vtol(Tmp);
     Tmp = newmat(Size,Size,Tmp);
     return base_replace(Tmp,taka_pfp_rule(cc,Size));
}

def taka_pfp_fun(Bb,Size) {
  Tmp = 1;
  for (K=1; K<=Size; K++) {
    Tmp = Tmp*(x+taka_pfp_iv(Bb,K));
  }
  Ans = newvect(Size);
  for (K=0; K<Size; K++) {
    Ans[K] = coef(Tmp,K,x);
  }
  return Ans;
}

def taka_pfp_nvec(K,Size) {
  Tmp = newvect(Size);
  Tmp[K] = 1/x;
  return Tmp;
}

def taka_pfp_nnvec(Aaa,Ccc,Size) {
  Tmp = newvect(Size);
  for (K=0; K<Size; K++) {
    Tmp[K] = (Aaa[K]*x-Ccc[K])/(x*(1-x));
  }
  return Tmp;
}

/*
 Taka_pfp_rule1 
=[[aa1,1/5], [aa2,1/5], [aa3,1/5], [aa4,1/5],[aa5,1],
  [cc1,2/5], [cc2, 3/5],[cc3,4/5], [cc4,1]]
*/
def taka_pfp_1_omega() {
   Taka_pfp_rule1 
   =[[aa1,1/5], [aa2,1/5], [aa3,1/5], [aa4,1/5],
     [cc1,2/5], [cc2,3/5], [cc3,4/5]]$
   Omega = taka_pfp_omega(4);
   Omega = base_replace(Omega, Taka_pfp_rule1);
   return Omega;
}

end$