x1;x2;x3;x4;x5;

qf_find_primitive_zeros_mod_p_general(q, p) = {

  n = p^3+p^2+p^1+p^0;
  Z = matrix(n, 5);
  w = [1, 0, 0, 0, 0];
  i = 1;
  if(qfeval(q, w)%p == 0, Z[i, ] = w; i+= 1);
  for(a=0, p-1,
    w = [a, 1, 0, 0, 0];
    if(qfeval(q, w)%p == 0, Z[i, ] = w; i+= 1);
    for(b = 0, p-1,
      w = [a, b, 1, 0, 0];
      if(qfeval(q, w)%p == 0, Z[i, ] = w; i+= 1);
      for(c = 0, p-1,
        w = [a, b, c, 1, 0];
        if(qfeval(q, w)%p == 0, Z[i, ] = w; i+= 1);
        for(d = 0, p-1,
          w = [a, b, c, d, 1];
          if(qfeval(q, w)%p == 0, Z[i, ] = w; i+= 1);
        );
      );
    );
  );
  Z;
};

divisorssquarefree(n) = divisors(vecprod(factor(n)[, 1]));


primitive(v, n:int) = {

  my(i);
  i = #v;
  while(gcd(v[i], n) > 1, i-= 1);
  (v*v[i]^-1)%n;

};


qf_isotropic_vector_random(q, p:int) = {

  my(v);
  while(1, v = vector(5, i, random(p))~; if(v != 0 && (qfeval(q, v)/2)%p == 0, break));
  Mat(v)%p;

};

qf_hensel_lifting(q, p, X, Z) = {

  my(x1, x2, x1p, x2p, z1, z2, z1p, z2p, z1pp, z2pp, M);
  x1 = X[, 1]; x2 = X[, 2];
  z1 = Z[, 1]; z2 = Z[, 2];
  x1p = (x1 - qfeval(q, x1)/2*z1)%p^2;
  x2p = (x2 - qfeval(q, x2)/2*z2 - qfeval(q, x1, x2)*z1)%p^2;
  z1p = (z1 - qfeval(q, z1)/2*x1)%p^2;
  z2p = (z2 - qfeval(q, z2)/2*x2 - qfeval(q, z1, z2)*x1)%p^2;
  z1pp = (z1p + (1 - qfeval(q, x1p, z1p))*z1p - qfeval(q, x2p, z1p)*z2p)%p^2;
  z2pp = (z2p - qfeval(q, x1p, z2p)*z1p + (1 - qfeval(q, x2p, z2p))*z2p)%p^2;
  M = matrix(5, 4);
  M[, 1] = x1p; M[, 2] = x2p; M[, 3] = z1pp; M[, 4] = z2pp;
  return(M);

};

qf_lifts_with_fixed_complement(q, p, X, Z) = {

  x1 = X[, 1]; x2 = X[, 2];
  z1 = Z[, 1]; z2 = Z[, 2];

};

qf_hensel_lifting_vector(q, p:int, v) = {

  my(w, a:int);
  while(1, w = matrix(5, 1, i, j, random(p)); if(qfeval(q, v[,1], w[,1])%p != 0, break));
  w*= qfeval(q, v[,1], w[,1])^1%p;
  w = (w - qfeval(q, w[,1])/2*v);
  a = -(qfeval(q, v[,1])/2)*qfeval(q, v[,1], w[,1])^-1%p^2;
  (v + a*w)%p^2;

}

qf_isotropic_plane_random(q, p:int) = {

  my(v1, v2);
  v1 = qf_isotropic_vector_random(q, p);
  while(q, v2 = qf_orthogonal_random(q, p, v1); if((qfeval(q, v2)/2)%p == 0 && primitive(v1[,1], p) != primitive(v2[,1], p), break));
  matconcat([v1, v2]);

};

qf_orthogonal_random(q, p:int, X) = {

  my(K, v, n:int);
  K = matkermod((q*Mat(X))~, p);
  n = matsize(K)[2];
  while(1, v = K*matrix(n, 1, i, j, random(p)); if(v%p != 0, break));
  v%p; 

};

qf_hyperbolic_complement(q, p:int, X) = {

  my(v3, v4);
  v1 = X[, 1];
  v2 = X[, 2];
  while(1, v3 = qf_orthogonal_random(q, p, v1)[,1]; if(qfeval(q, v2, v3)%p != 0 && primitive(v2, p) != primitive(v3, p), break));
  v3*= qfeval(q, v2, v3)^-1%p;
  v3 = (v3 - qfeval(q, v3)/2*v2)%p;
  while(1, v4 = qf_orthogonal_random(q, p, matconcat([v2, v3]))[,1]; if(qfeval(q, v1, v4)%p != 0 && primitive(v1, p) != primitive(v4, p), break));
  v4*= qfeval(q, v1, v4)^-1%p;
  v4 = (v4 - qfeval(q, v4)/2*v1)%p;
  matconcat([v4, v3]);
};

qf_basis(q, p:int) = {

  my(X, Y, v5);
  X = qf_isotropic_plane_random(q, p);
  Y = qf_hyperbolic_complement(q, p, X);
  v5 = qf_orthogonal_random(q, p, matconcat([X, Y]));
  matconcat([X, Y[, 2], Y[, 1], v5]);

};

qf_lifts_with_fixed_complement(q, p, X, Z) = {

  my(x1, x2, z1, z2, L, x1p, x2p, M);
  x1 = X[, 1]; x2 = X[, 2];
  z1 = Z[, 1]; z2 = Z[, 2];
  L = listcreate();
  for(a = 0, p - 1, x1p = x1+p*a*z2; x2p = x2 - p*a*z1; 
                    M = matrix(5, 2); M[, 1] = x1p; M[, 2] = x2p;
                    listput(L, M));
  L;

};

qf_p_neighbor(q, p:int, X) = {

 my(Kp);
 Kp = matkermod((q*Mat(X))~, p);
 Kp = matconcat([Kp, p*matid(5)]);
 Kp = mathnf(Kp);
 Kp = matconcat([X, p*Kp]);
 Kp = mathnf(Kp)/p;
 [Kp~*q*Kp, Kp];

};

qf_isotropic_planes(q, p) = {

  my(B, v1, v2, v3, v4, v5, L, qv5p, M, w1, w2);
  B = qf_basis(q, p);
  v1 = B[, 1]; v2 = B[, 2]; v3 = B[, 3]; v4 = B[, 4]; v5 = B[, 5];
  L = listcreate();
  qv5p = (qfeval(q, v5)/2)%p;
  listput(L, B[, 3..4]);
  for(a = 0, p - 1,
    w1 = v2 - a^2*qv5p*v3 + a*v5;
    M = matrix(5, 2);
    M[, 1] = w1%p; M[, 2] = v4;
    listput(L, M);
    for(b = 0, p - 1,
      w1 = v1 + a*v2 - b^2*qv5p*v4 + b*v5;
      w2 = v3 - a*v4;
      M = matrix(5, 2);
      M[, 1] = w1%p; M[, 2] = w2%p;
      listput(L, M);
      for(c = 0, p - 1,
        w1 = v1 - (2*a*c*qv5p + b)*v3 - a^2*qv5p*v4 + a*v5;
        w2 = v2 - c^2*qv5p*v3 + b*v4 + c*v5;
        M = matrix(5, 2);
        M[, 1] = w1%p; M[, 2] = w2%p;
        listput(L, M);
      );
    );
  );
  L;

}



qf_p2_neighbors(q, p:int) = {

  my(M, L, X, Z, H, X1, Z1, X2, lifts, isoplanes);
  L = listcreate();
  isoplanes = qf_isotropic_planes(q, p);
  for(i = 1, #isoplanes,
    X = isoplanes[i];
    Z = qf_hyperbolic_complement(q, p, X);
    H = qf_hensel_lifting(q, p, X, Z);
    X1 = H[, 1..2];
    Z1 = H[, 3..4];
    lifts = qf_lifts_with_fixed_complement(q, p, X1, Z1);
    for(j = 1, #lifts,
      X2 = lifts[j];
      listput(L, qf_p_neighbor(q, p, X2));
    );
  );
  L;

};

qf_p_neighbors(q, p:int) = {

  my(M, L, X);
  M = qf_basis(q, p);
  L = listcreate();
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  listput(L, qf_p_neighbor(q, p, X));
  for(i = 0, p-1, 
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
    listput(L, qf_p_neighbor(q, p, X));
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
    listput(L, qf_p_neighbor(q, p, X));
    for(j = 0, p-1, 
      X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
      listput(L, qf_p_neighbor(q, p, X));
      if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        listput(L, qf_p_neighbor(q, p, X));
        );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        listput(L, qf_p_neighbor(q, p, X));
      );
    );
  );
  L;

};

qf_p_neighbors_singular(q, p) = {

  my(M, L, X);
  L = listcreate();
  X = [1, 0, 0, 0, 0]~;
  if(qfeval(q, X)/2%p == 0 && q*X%p != 0, listput(L, qf_p_neighbor(q, p, qf_hensel_lifting_vector(q, p, Mat(X)))));
  for(a = 0, p-1,
    X = [a, 1, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0, listput(L, qf_p_neighbor(q, p, qf_hensel_lifting_vector(q, p, Mat(X)))));
    for(b = 0, p-1,
      X = [a, b, 1, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0, listput(L, qf_p_neighbor(q, p, qf_hensel_lifting_vector(q, p, Mat(X)))));
      for(c = 0, p-1,
        X = [a, b, c, 1, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0, listput(L, qf_p_neighbor(q, p, qf_hensel_lifting_vector(q, p, Mat(X)))));
        for(d = 0, p-1,
          X = [a, b, c, d, 1]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0, listput(L, qf_p_neighbor(q, p, qf_hensel_lifting_vector(q, p, Mat(X)))));
        );
      );
    );
  );
  L;
  
};


qf_p_neighbors_fun(q, p, f, ac) = {

  my(M, L, X);
  M = qf_basis(q, p);
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  f(ac, qf_p_neighbor(q, p, X));
  for(i = 0, p-1,
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
    f(ac, qf_p_neighbor(q, p, X));
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
    f(ac, qf_p_neighbor(q, p, X));
    for(j = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
      f(ac, qf_p_neighbor(q, p, X));
      if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        f(ac, qf_p_neighbor(q, p, X))
        );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        f(ac, qf_p_neighbor(q, p, X))
      );
    );
  );
  ac;

};


qf_is_in_genus(G, q) = {
  
  my(isin, iso);
  for(i=1, #G, iso = qfisom(G[i], q); if(iso != 0, isin = 1; return(1)));
  return(0);

}

nu(n:int, d:int) = (-1)^(matsize(factor(gcd(n, d)))[1]);


qf_tp(X, q, p, qiso, qtheta) = {

  my(tp, A, B, r);
  X = qf_hensel_lifting_vector(q, p, Mat(X));
  [r, A] = qf_p_neighbor(q, p, X);
  A*= matdet(A);
  tp = 0;
  if(qfrep(r, 2) == qtheta, B = qfisom(qiso, r); B*= matdet(B); tp+= nu(qf_spinor_norm(q, A*B), 167));
  tp;

}

qf_tp2(tp, X, q, p, qiso, qtheta) = {

  my(tp, A, B, r);
  [r, A] = X;
  A*= matdet(A);
  if(qfrep(r, 2) == qtheta, B = qfisom(qiso, r); B*= matdet(B); tp+= nu(qf_spinor_norm(q, A*B), 167));
  tp;

}


qf_167_2(p) = {

  my(q, qiso, qtheta, f);
  q = [4, 1, 1, 2, 1; 1, 4, 2, 1, -1; 1, 2, 4, 2, 0; 2, 1, 2, 4, 1; 1, -1, 0, 1, 4];
  qiso = qfisominit(q);
  qtheta = qfrep(q, 2);
  f(tp, X) = qf_tp2(tp, X, q, p, qiso, qtheta);
  qf_p_neighbors_fun(q, p, f, 0);

}


qf_q_i(q, p, M, i) = {

   my(qiso, qtheta, tp, X);
   qiso = qfisominit(q);
   qtheta = qfrep(q, 10);
   tp = 0;
   X = (i*M[, 1]+M[, 2])%p;
    tp+= qf_tp(X, q, p, qiso, qtheta);
    X = (i*M[, 1]+M[, 3])%p;
    tp+= qf_tp(X, q, p, qiso, qtheta);
    for(j = 0, p-1,
      X = ((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p;
      tp+= qf_tp(X, q, p, qiso, qtheta);
      if(j != 0, 
        X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        tp+= qf_tp(X, q, p, qiso, qtheta);
        );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        tp+= qf_tp(X, q, p, qiso, qtheta);
      );
    );
   tp;

};

qf_q_par(q, p) = {

  my(qiso, qtheta, M, X, tp);
  qiso = qfisominit(q);
  qtheta = qfrep(q, 10);
  tp = 0;
  M = qf_basis(q, p);
  X = (M[, 1]);
  tp+= qf_tp(X, q, p, qiso, qtheta);
  tp+= parsum(i=0, p-1, qf_q_i(q, p, M, i));
  tp;
 
}

qf_167(p) = {

  my(q, qiso, qtheta, M, X, A, B, tp, r);
  q = [4, 1, 1, 2, 1; 1, 4, 2, 1, -1; 1, 2, 4, 2, 0; 2, 1, 2, 4, 1; 1, -1, 0, 1, 4];
  qiso = qfisominit(q);
  qtheta = qfrep(q, 2);
  tp = 0;
  M = qf_basis(q, p);
  X = (M[, 1]);
  tp+= qf_tp(X, q, p, qiso, qtheta);
  for(i = 0, p-1,
    X = (i*M[, 1]+M[, 2])%p;
    tp+= qf_tp(X, q, p, qiso, qtheta);
    X = (i*M[, 1]+M[, 3])%p;
    tp+= qf_tp(X, q, p, qiso, qtheta);
    for(j = 0, p-1,
      X = ((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p;
      tp+= qf_tp(X, q, p, qiso, qtheta);
      if(j != 0, 
        X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        tp+= qf_tp(X, q, p, qiso, qtheta);
        );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        tp+= qf_tp(X, q, p, qiso, qtheta);
      );
    );
  );
  tp;

}

qf_reduce(q) = {

  my(A);
  A = qflll(q);
  A*= matdet(A);
  [A~*q*A, A];

}

qf_spinor_norm(q, A) = {

  my(v, s, n);
  A*= matdet(A);
  s = 1;
  while(A!=s,
    v = matimage(A-s)[, 1];
    n = qfeval(q, v);
    s*= n;
    A*= (n - 2*v*v~*q);
  );
  core(numerator(s))*core(denominator(s));

}



/* HAY UN ERROR EN ESTA FUNCION, NO ENCUENTRA EL GENERO ENTERO (DEPENDE DE N) */

qf_find_genus(q, p = 2, N = 20) = {

  my(G, L, Lqq, qq, r, M);
  G = listcreate();
  L = listcreate();
  listput(G, q);
  listput(L, q);
  M = Map();
  mapput(M, qfrep(q, N), q);
  while(#L > 0,
    qq = L[1];
    Lqq = qf_p_neighbors(qq, p);
    for(i = 1, #Lqq, r = qf_reduce(Lqq[i][1])[1]; if(!mapisdefined(M, qfrep(r, N)), mapput(M, qfrep(r, N), r); listput(L, r);listput(G, r)));
    listpop(L, 1);
  );
  G;

}

omf_init(q, p = 2, N = 20) = {

  my(G, L, M, T, k, qq, Lqq, r, A, rt, new, LT, D);
  D = matdet(q)/2;
  G = List();
  L = List();
  M = Map();
  T = Map();
  k = 2;
  listput(L, q);
  listput(G, [q, matid(5)]);
  mapput(M, q, 1);
  mapput(T, qfrep(q, N), List([q]));
  while(#L > 0,
    qq = L[1];
    if(D%p == 0, Lqq = qf_p_neighbors_singular(qq, p), Lqq = qf_p_neighbors(qq, p));
    for(i = 1, #Lqq, 
      [r, A] = qf_reduce(Lqq[i][1]);
      rt = qfrep(r, N);
      if(!mapisdefined(T, rt),
        mapput(T, rt, List([r]));
        mapput(M, r, k);
        k+= 1;
        listput(G, [r, G[mapget(M, qq)][2]*Lqq[i][2]*A]);
        listput(L, r);
      ,
        new = 1;
        for(j=1, #mapget(T, rt), if(0!=qfisom(r, mapget(T, rt)[j]), new = 0; break));
        if(new,
          mapput(M, r, k);
          k+= 1;
          listput(G, [r, G[mapget(M, qq)][2]*Lqq[i][2]*A]);
          LT = mapget(T, rt);
          listput(LT, r);
          mapput(T, rt, LT);
          listput(L, r);
        );
      );
    );
    listpop(L, 1);
  );
  [G, M, T, N, p];

}

omf_get1(omf, X, q, p) = {

  my(r, A, N, G, M, T, Tr, rr, qq);
  [G, M, T, N] = omf;
  [r, A] = qf_p_neighbor(q, p, X);
  Tr = mapget(T, qfrep(r, N));
  if(#Tr>1,
    for(i=1, #Tr, if(0!=qfisom(r, Tr[i]), qq = Tr[i]; break));
    ,
    qq = Tr[1];
  );
  mapget(M, qq);

}

omf_get(omf, X, q, p) = {

  my(r, A, Tr, rr, qq);
  [r, A] = qf_p_neighbor(q, p, X);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, if(0!=qfisom(r, Tr[i]), qq = Tr[i]; break));
    ,
    qq = Tr[1];
  );
  mapget(omf[2], qq);

}


omf_get2(omf, r) = {

  my(Tr, rr);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, if(0!=qfisom(r, Tr[i]), rr = Tr[i]; break));
    ,
    rr = Tr[1];
  );
  mapget(omf[2], rr);

};

omf_get2d(omf, row, r, A) = {

  my(Tr, iso, col, A1, s);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(0!=iso, rr = Tr[i]; break));
    ,
    rr = Tr[1]; iso = qfisom(r, Tr[1]);
  );
  col = mapget(omf[2], rr);
  A1 = omf[1][row][2]*A*iso^-1*omf[1][col][2]^-1;
  s = qf_spinor_norm(omf[1][1][1], A1^-1);
  return([col, s]);

};

omf_get2nd(omf, row, r, A) = {

  my(Tr, iso, col, A1, s);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(0!=iso, rr = Tr[i]; break));
    ,
    rr = Tr[1]; iso = qfisom(r, Tr[1]);
  );
  col = mapget(omf[2], rr);
  iso*= matdet(iso);
  A1 = omf[1][row][2]*A*iso^-1*omf[1][col][2]^-1;
  s = qf_spinor_norm(omf[1][1][1], A1^-1);
  A1 = A*iso^-1;
  return([col, s, A1]);

};


omf_tp1_singular(omf, p) = {
 
  my(G, M, T, N, n, tp1, row, col, X, q);
  [G, M, T, N] = omf;
  n = #G;
  tp1 = matrix(n);
  for(row = 1, n,
    q = G[row][1];
    X = [1, 0, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0, 
    col = omf_get(omf, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); tp1[row, col]+ =1);
    for(a = 0, p-1,
      X = [a, 1, 0, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      col = omf_get(omf, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); tp1[row, col]+ =1);
      for(b = 0, p-1,
        X = [a, b, 1, 0, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        col = omf_get(omf, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); tp1[row, col]+ =1);
        for(c = 0, p-1,
          X = [a, b, c, 1, 0]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          col = omf_get(omf, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); tp1[row, col]+ =1);
          for(d = 0, p-1,
            X = [a, b, c, d, 1]~;
            if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            col = omf_get(omf, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); tp1[row, col]+ =1);
          );
        );
      );
    );
  );
  tp1;

};

omf_tp1d_init_singular(omf, p) = {

  my(D, spins, tp1, G, M, T, N, n, q, X, col, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp1 = matrix(#omf[1], #omf[1], i, j, vector(#spins));
  [G, M, T, N] = omf;
  n = #G;
  for(row = 1, n,
    q = G[row][1];
    X = [1, 0, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s] = omf_getd(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      tp1[row, col][scolrow]+ =1);
    for(a = 0, p-1,
      X = [a, 1, 0, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s] = omf_getd(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        tp1[row, col][scolrow]+ =1);
      for(b = 0, p-1,
        X = [a, b, 1, 0, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s] = omf_getd(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p); 
          scolrow = vecsearch(spins, s);
          tp1[row, col][scolrow]+ =1);
        for(c = 0, p-1,
          X = [a, b, c, 1, 0]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s] = omf_getd(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            tp1[row, col][scolrow]+ =1);
          for(d = 0, p-1,
            X = [a, b, c, d, 1]~;
            if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
              [col, s] = omf_getd(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
              scolrow = vecsearch(spins, s);
              tp1[row, col][scolrow]+ =1);
          );
        );
      );
    );
  );
  [tp1, spins];


};


omf_tp2(omf, p) = {

  my(G, M, T, N, n, tp2, q, col, p2neighs, r);
  [G, M, T, N] = omf;
  n = #G;
  tp2 = matrix(n);
  for(row = 1, n,
    q = G[row][1];
    p2neighs = qf_p2_neighbors(q, p);
    for(i = 1, #p2neighs,
      r = p2neighs[i][1];
      col = omf_get2(omf, r);
      tp2[row, col]+= 1;
    );
  );  
  tp2;

};

omf_tp2d_init(omf, p) = {

  my(D, spins, tp2, G, M, T, N, n, q, p2neighs, r, A, col, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp2 = matrix(#omf[1], #omf[1], i, j, vector(#spins));
  [G, M, T, N] = omf;
  n = #G;
  for(row = 1, n,
    q = G[row][1];
    p2neighs = qf_p2_neighbors(q, p);
    for(i = 1, #p2neighs,
      [r, A] = p2neighs[i];
      [col, s] = omf_get2d(omf, row, r, A);
      scolrow = vecsearch(spins, s);
      tp2[row, col][scolrow]+= 1;
    );
  );
  [tp2, spins];

};

omf_tp2nd_init(omf, p, n) = {

  my(D, spins, tp2, G, M, T, N, m, q, p2neighs, r, A, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp2n = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  \\tp2 = matrix(#omf[1], #omf[1], i, j, vector(#spins));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    p2neighs = qf_p2_neighbors(q, p);
    for(i = 1, #p2neighs,
      [r, A] = p2neighs[i];
      [col, s, AA] = omf_get2nd(omf, row, r, A);
      scolrow = vecsearch(spins, s);
      tp2n[row, col][scolrow]+= action_n(AA^-1, n);
    );
  );
  [tp2n, spins];


};

omf_tp2abd_init(omf, p, a, b) = {

  my(D, spins, tp2ab, G, M, T, N, m, q, p2neighs, r, A, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp2ab = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(ngl(a, b))));
  \\tp2 = matrix(#omf[1], #omf[1], i, j, vector(#spins));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    p2neighs = qf_p2_neighbors(q, p);
    for(i = 1, #p2neighs,
      [r, A] = p2neighs[i];
      [col, s, AA] = omf_get2nd(omf, row, r, A);
      scolrow = vecsearch(spins, s);
      tp2ab[row, col][scolrow]+= actionab(AA^-1, a, b);
    );
  );
  [tp2ab, spins];

};


omf_tp1(omf, p) = {

  my(G, M, T, n, tp1, row, col, X, r, A, N, q);
  [G, M, T, N] = omf;
  n = #G;
  tp1 = matrix(n);
  for(row = 1, n,
    q = G[row][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    col = omf_get(omf, X, q, p);
    tp1[row, col]+= 1;
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      col = omf_get(omf, X, q, p);
      tp1[row, col]+= 1;
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      col = omf_get(omf, X, q, p);
      tp1[row, col]+= 1;
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        col = omf_get(omf, X, q, p);
        tp1[row, col]+= 1;
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          col = omf_get(omf, X, q, p);
          tp1[row, col]+= 1;
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          col = omf_get(omf, X, q, p);
          tp1[row, col]+= 1;
        );
      );
    );
  );
  tp1;
}

omf_tp1d_init(omf, p) = {

  my(spins, tp1, G, M, T, n, row, col, X, r, A, N, q, s, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp1 = matrix(#omf[1], #omf[1], i, j, vector(#spins));
  [G, M, T, N] = omf;
  n = #G;
  for(row = 1, n,
    q = G[row][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    [col, s] = omf_getd(omf, row,  X, q, p);
    scolrow = vecsearch(spins, s);
    tp1[row, col][scolrow]+= 1;
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      [col, s] = omf_getd(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tp1[row, col][scolrow]+= 1;
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      [col, s] = omf_getd(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tp1[row, col][scolrow]+= 1;
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        [col, s] = omf_getd(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        tp1[row, col][scolrow]+= 1;
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s] = omf_getd(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tp1[row, col][scolrow]+= 1;
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s] = omf_getd(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tp1[row, col][scolrow]+= 1;
        );
      );
    );
  );
  [tp1, spins];
};

omf_getdn(omf, row, X, q, p) = {

  my(r, A, Tr, rr, qq, col, s, A1, iso);
  [r, A] = qf_p_neighbor(q, p, X);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(0!=iso, qq = Tr[i]; break));
    ,
    qq = Tr[1]; iso = qfisom(r, Tr[1]);
  );
  col = mapget(omf[2], qq);
  iso*= matdet(iso);
  A1 = omf[1][row][2]*A*iso^-1*omf[1][col][2]^-1;
  s = qf_spinor_norm(omf[1][1][1], A1^-1);
  A1 = A*iso^-1;
  return([col, s, A1]);

}


omf_tpnd_init(omf, p, n) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    [col, s, AA] = omf_getdn(omf, row,  X, q, p);
    scolrow = vecsearch(spins, s);
    tpn[row, col][scolrow]+= action_n(AA^-1, n);
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= action_n(AA^-1, n);
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= action_n(AA^-1, n);
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        tpn[row, col][scolrow]+= action_n(AA^-1, n);
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= action_n(AA^-1, n);
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= action_n(AA^-1, n);
        );
      );
    );
  );
  [tpn*p^n, spins];

};

omf_tpabd_init(omf, p, a, b) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(ngl(a, b))));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    [col, s, AA] = omf_getdn(omf, row,  X, q, p);
    scolrow = vecsearch(spins, s);
    tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= actionab(AA^-1, a, b);
        );
      );
    );
  );
  [tpn, spins];

};


omf_tpnd_row_init(omf, p, n, row) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  [G, M, T, N] = omf;
  \\m = #G;
  q = G[row][1];
  M = qf_basis(q, p);
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  [col, s, AA] = omf_getdn(omf, row,  X, q, p);
  scolrow = vecsearch(spins, s);
  tpn[col][scolrow]+= action_n(AA^-1, n);
  for(i = 0, p-1,
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    tpn[col][scolrow]+= action_n(AA^-1, n);
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    tpn[col][scolrow]+= action_n(AA^-1, n);
    for(j = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[col][scolrow]+= action_n(AA^-1, n);
      if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        tpn[col][scolrow]+= action_n(AA^-1, n);
      );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        tpn[col][scolrow]+= action_n(AA^-1, n);
      );
    );
  );
  [tpn*p^n, spins];
};

omf_tpnd_row_init_i(omf, p, n, row, i, M) = {

  my(D, spins, tpn, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  \\m = #G;
  q = omf[1][row][1];
  X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
  [col, s, AA] = omf_getdn(omf, row, X, q, p);
  scolrow = vecsearch(spins, s);
  tpn[col][scolrow]+= action_n(AA^-1, n);
  X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
  [col, s, AA] = omf_getdn(omf, row, X, q, p);
  scolrow = vecsearch(spins, s);
  tpn[col][scolrow]+= action_n(AA^-1, n);
  for(j = 0, p-1,
    X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    tpn[col][scolrow]+= action_n(AA^-1, n);
    if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
      X = qf_hensel_lifting_vector(q, p, Mat(X));
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[col][scolrow]+= action_n(AA^-1, n);
    );
    for(k=1, p-1,
      X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
      X = qf_hensel_lifting_vector(q, p, Mat(X));
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      tpn[col][scolrow]+= action_n(AA^-1, n);
    );
  );
  tpn*p^n;

};

omf_tpnd_row_init_par(omf, p, n, row) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  \\[G, M, T, N] = omf;
  \\m = #G;
  q = omf[1][row][1];
  M = qf_basis(q, p);
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  [col, s, AA] = omf_getdn(omf, row,  X, q, p);
  scolrow = vecsearch(spins, s);
  tpn[col][scolrow]+= action_n(AA^-1, n)*p^n;
  \\print(omf_tpnd_row_init_i(omf, p, n, row, 0, M));
  tpn+= parsum(i = 0, p-1, omf_tpnd_row_init_i(omf, p, n, row, i, M));
  [tpn, spins];


};


omf_tpnds_init(omf, p, Bn) = {

  my(D, spins, tpns, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    [col, s, AA] = omf_getdn(omf, row,  X, q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [col, s, AA] = omf_getdn(omf, row, X, q, p);
          scolrow = vecsearch(spins, s);
          for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn));
        );
      );
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*= p^nn);
  [tpns, spins];
};

omf_tpnds_row_init(omf, p, Bn, row) = {

  my(D, spins, tpns, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  q = G[row][1];
  M = qf_basis(q, p);
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  [col, s, AA] = omf_getdn(omf, row,  X, q, p);
  scolrow = vecsearch(spins, s);
  for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
  for(i = 0, p-1,
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
    X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
    for(j = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
      if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
      );
      for(k=1, p-1,
        X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
        X = qf_hensel_lifting_vector(q, p, Mat(X));
        [col, s, AA] = omf_getdn(omf, row, X, q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
      );
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*= p^nn);
  [tpns, spins];
};

omf_tpnds_row_init_i(omf, p, Bn, row, i, M) = {

  my(D, spins, tpns, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  q = omf[1][row][1];
  X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
  [col, s, AA] = omf_getdn(omf, row, X, q, p);
  scolrow = vecsearch(spins, s);
  for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
  X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
  [col, s, AA] = omf_getdn(omf, row, X, q, p);
  scolrow = vecsearch(spins, s);
  for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
  for(j = 0, p-1,
    X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
    [col, s, AA] = omf_getdn(omf, row, X, q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
    if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
      X = qf_hensel_lifting_vector(q, p, Mat(X));
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
    );
    for(k=1, p-1,
      X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
      X = qf_hensel_lifting_vector(q, p, Mat(X));
      [col, s, AA] = omf_getdn(omf, row, X, q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*= p^nn);
  tpns;

};

omf_tpnds_row_init_par(omf, p, Bn, row) = {

  my(D, spins, tpns, G, M, T, N, m, q, X, col, s, AA, scolrow);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  q = omf[1][row][1];
  M = qf_basis(q, p);
  X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
  [col, s, AA] = omf_getdn(omf, row,  X, q, p);
  scolrow = vecsearch(spins, s);
  for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
  for(nn = 0, Bn, tpns[nn + 1]*= p^nn);
  tpns+= parsum(i=0, p-1, omf_tpnds_row_init_i(omf, p, Bn, row, i, M));
  [tpns, spins];
  
};




omf_tpnd_init_singular(omf, p, n) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    X = [1, 0, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= action_n(AA^-1, n));
    for(a = 0, p-1,
      X = [a, 1, 0, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        tpn[row, col][scolrow]+= action_n(AA^-1, n));
      for(b = 0, p-1,
        X = [a, b, 1, 0, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= action_n(AA^-1, n));
        for(c = 0, p-1,
          X = [a, b, c, 1, 0]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            tpn[row, col][scolrow]+= action_n(AA^-1, n));
          for(d = 0, p-1,
            X = [a, b, c, d, 1]~;
            if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
              [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
              scolrow = vecsearch(spins, s);
              tpn[row, col][scolrow]+= action_n(AA^-1, n));
          );
        );
      );
    );
  );
  [tpn*p^n, spins];

};

omf_tpabd_init_singular(omf, p, aa, bb) = {

  my(D, spins, tpn, G, M, T, N, m, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(ngl(aa, bb))));
  [G, M, T, N] = omf;
  m = #G;
  for(row = 1, m,
    q = G[row][1];
    X = [1, 0, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      tpn[row, col][scolrow]+= actionab(AA^-1, aa, bb));
    for(a = 0, p-1,
      X = [a, 1, 0, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        tpn[row, col][scolrow]+= actionab(AA^-1, aa, bb));
      for(b = 0, p-1,
        X = [a, b, 1, 0, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          tpn[row, col][scolrow]+= actionab(AA^-1, aa, bb));
        for(c = 0, p-1,
          X = [a, b, c, 1, 0]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            tpn[row, col][scolrow]+= actionab(AA^-1, aa, bb));
          for(d = 0, p-1,
            X = [a, b, c, d, 1]~;
            if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
              [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
              scolrow = vecsearch(spins, s);
              tpn[row, col][scolrow]+= actionab(AA^-1, aa, bb));
          );
        );
      );
    );
  );
  [tpn*p^aa, spins];

};

omf_tpnd_row_init_singular(omf, p, n, row) = {

  my(D, spins, tpn, G, M, T, N, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpn = vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, n)))));
  [G, M, T, N] = omf;
  q = G[row][1];
  X = [1, 0, 0, 0, 0]~;
  if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
    [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
    scolrow = vecsearch(spins, s);
    tpn[col][scolrow]+= action_n(AA^-1, n));
  for(a = 0, p-1,
    X = [a, 1, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      tpn[col][scolrow]+= action_n(AA^-1, n));
    for(b = 0, p-1,
      X = [a, b, 1, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        tpn[col][scolrow]+= action_n(AA^-1, n));
      for(c = 0, p-1,
        X = [a, b, c, 1, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          tpn[col][scolrow]+= action_n(AA^-1, n));
        for(d = 0, p-1,
          X = [a, b, c, d, 1]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            tpn[col][scolrow]+= action_n(AA^-1, n));
        );
      );
    );
  );
  [tpn*p^n, spins];

};


omf_tpnds_init_singular(omf, p, Bn) = {
  
  my(D, spins, tpns, G, M, T, N, m, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, matrix(#omf[1], #omf[1], i, j, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  m = #G; 
  for(row = 1, m,
    q = G[row][1];
    X = [1, 0, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn)));
    for(a = 0, p-1, 
      X = [a, 1, 0, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn)));
      for(b = 0, p-1, 
        X = [a, b, 1, 0, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn)));
        for(c = 0, p-1, 
          X = [a, b, c, 1, 0]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn)));
          for(d = 0, p-1, 
            X = [a, b, c, d, 1]~;
            if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
              [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
              scolrow = vecsearch(spins, s);
              for(nn=0, Bn, tpns[nn+1][row, col][scolrow]+= action_n(AA^-1, nn)));
          );
        );
      );
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*=p^nn);
  [tpns, spins];

};

omf_tpnds_row_init_singular(omf, p, Bn, row) = {

  my(D, spins, tpns, G, M, T, N, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  q = G[row][1];
  X = [1, 0, 0, 0, 0]~;
  if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
    [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
  for(a = 0, p-1,
    X = [a, 1, 0, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
    for(b = 0, p-1,
      X = [a, b, 1, 0, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
      for(c = 0, p-1,
        X = [a, b, c, 1, 0]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
        for(d = 0, p-1,
          X = [a, b, c, d, 1]~;
          if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
            [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
            scolrow = vecsearch(spins, s);
            for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
        );
      );
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*=p^nn);
  [tpns, spins];

};

omf_tpnds_row_init_singular_a(omf, p, Bn, row, a, M) = {

  my(D, spins, tpns, q, X, col, scolrow, AA, s);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  q = omf[1][row][1];
  X = [a, 1, 0, 0, 0]~;
  if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
    [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
  for(b = 0, p-1,
    X = [a, b, 1, 0, 0]~;
    if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
      [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
      scolrow = vecsearch(spins, s);
      for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
    for(c = 0, p-1,
      X = [a, b, c, 1, 0]~;
      if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
        [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
        scolrow = vecsearch(spins, s);
        for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
      for(d = 0, p-1,
        X = [a, b, c, d, 1]~;
        if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
          [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
          scolrow = vecsearch(spins, s);
          for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
      );
    );
  );
  for(nn = 0, Bn, tpns[nn + 1]*=p^nn);
  tpns;

};

omf_tpnds_row_init_singular_par(omf, p, Bn, row) = {

  my(D, spins, tpns, G, M, T, N, q, X, col, scolrow, AA, s);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tpns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  q = G[row][1];
  X = [1, 0, 0, 0, 0]~;
  if(qfeval(q, X)/2%p == 0 && q*X%p != 0,
    [col, s, AA] = omf_getdn(omf, row, qf_hensel_lifting_vector(q, p, Mat(X)), q, p);
    scolrow = vecsearch(spins, s);
    for(nn=0, Bn, tpns[nn+1][col][scolrow]+= action_n(AA^-1, nn)));
  for(nn = 0, Bn, tpns[nn + 1]*=p^nn);
  tpns+= parsum(a = 0, p - 1, omf_tpnds_row_init_singular_a(omf, p, Bn, row, a, M));
  [tpns, spins];


};




omf_tp2nds_row_init(omf, p, Bn, row) = {

  my(D, spins, tp2ns, G, M, T, N, q, X, col, scolrow, AA);
  D = p*matdet(omf[1][1][1])/2*omf[5];
  spins = divisorssquarefree(D);
  tp2ns = vector(Bn + 1, nn, vector(#omf[1], i, vector(#spins, i, matrix(abs(binomial(-5, nn-1))))));
  [G, M, T, N] = omf;
  q = G[row][1];
  p2neighs = qf_p2_neighbors(q, p);
  for(i = 1, #p2neighs,
    [r, A] = p2neighs[i];
    [col, s, AA] = omf_get2nd(omf, row, r, A);
    scolrow = vecsearch(spins, s);
    for(nn = 0, Bn, tp2ns[nn+1][col][scolrow]+= action_n(AA^-1, nn));
  );
  [tp2ns, spins];

};

omf_tpnd(tpndinit, d, invnd, spins) = {

  my(tpnd, m, Mij, i0, j0);
  m = sum(i=1, #invnd, matsize(invnd[i][2])[2]);
  tpnd = matrix(m);
  i0 = 1;
  for(i = 1, #invnd,
    j0 = 1;
    for(j = 1, #invnd,
      Mij = matsolve(invnd[i][2], sum(s = 1, #spins, nu(d, spins[s])*tpndinit[invnd[i][1], invnd[j][1]][s])*invnd[j][2]);
      \\Mij = matsolve(invnd[i][2], vecsum(tpndinit[i, j])*invnd[j][2]);
      \\print(Mij);
      for(ii = i0, i0 + matsize(Mij)[1]-1,
        for(jj = j0, j0 + matsize(Mij)[2]-1, tpnd[ii, jj] = Mij[ii-i0+1, jj-j0+1]);
      );
      \\tpnd[i0..i0+matsize(Mij)[1], j0..j0+matsize(Mij)[2]] = Mij;
      j0+= matsize(Mij)[2];
      \\print([i, j, matsize(Mij)]);
    );
    i0+= matsize(Mij)[1];
  );
  tpnd;

}

omf_tpnd_row(tpndinit, d, invnd, spins, row) = {

  my(tpnd, m, Mij, j0, row1);
  m = sum(i=1, #invnd, matsize(invnd[i][2])[2]);
  row1 = vecsearch(vector(#invnd, i, invnd[i][1]), row);
  tpnd = matrix(matsize(invnd[row1][2])[2], m);
  j0 = 1;
  for(j = 1, #invnd,
    Mij = matsolve(invnd[row1][2], sum(s = 1, #spins, nu(d, spins[s])*tpndinit[invnd[j][1]][s])*invnd[j][2]);
    for(ii = 1, matsize(Mij)[1],
      for(jj = j0, j0 + matsize(Mij)[2]-1, tpnd[ii, jj] = Mij[ii, jj-j0+1]);
    );
    j0+= matsize(Mij)[2];
  );
  tpnd;

};


omf_tp1d2(tp1dinit, d, invd, spins) = {

  my(tp1d);
  tp1d = matrix(#invd);
  for(row = 1, #invd,
    for(col = 1, #invd,
      tp1d[row, col] = sum(s = 1, #spins, nu(d, spins[s])*tp1dinit[invd[row], invd[col]][s]);
    );
  );
  tp1d

};

omf_invariants1d(omf, d) = {

  my(G, M, T, N, invariants);
  [G, M, T, N] = omf;
  invariants = List();
  for(i = 1, #omf[1], if(matsize(qf_invariant(omf[1][i][1], 0, d))[1] == 1, listput(invariants, i)));
  invariants;

};

omf_invariantsnd(omf, n, d) = {

  my(invariants, invq);
  invariants = List();
  for(i = 1, #omf[1], invq = qf_invariant(omf[1][i][1], n, d); if(matsize(invq)[1] != 0, listput(invariants, [i, invq])));;
  invariants; 

};

omf_invariantsabd(omf, a, b, d) = {

  my(invariants, invq);
  invariants = List();
  for(i = 1, #omf[1], invq = qf_invariantab(omf[1][i][1], a, b, d); if(matsize(invq)[1] != 0, listput(invariants, [i, invq])));;
  invariants;

};



listtomap(L) = {

  my(m);
  m = Map();
  for(i = 1, #L, mapput(m, L[i], i));
  m;

};

nprimedivisors(n) = {

  my(divs);
  divs = divisors(n);
  sum(i = 1, #divs, isprime(divs[i]));

}

nu(d, n) = (-1)^nprimedivisors(gcd(n, d));

omf_getd(omf, row, X, q, p) = {

  my(r, A, Tr, rr, qq, col, s, A1, iso);
  [r, A] = qf_p_neighbor(q, p, X);
  Tr = mapget(omf[3], qfrep(r, omf[4]));
  if(#Tr>1,
    for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(0!=iso, qq = Tr[i]; break));
    ,
    qq = Tr[1]; iso = qfisom(r, Tr[1]);
  );
  col = mapget(omf[2], qq);
  A1 = omf[1][row][2]*A*iso^-1*omf[1][col][2]^-1;
  s = qf_spinor_norm(omf[1][1][1], A1^-1);
  return([col, s]);

}

omf_tp1d(omf, invariants, p, d) = {

  my(G, M, M1, T, N, n, tpd, minvariants, col, col1, q, X, r, A, spin, iso, A1);
  [G, M1, T, N] = omf;
  n = #invariants;
  minvariants = listtomap(invariants);
  tpd = matrix(n);
  for(row = 1, n,
    q = G[invariants[row]][1];
    M = qf_basis(q, p);
    X = qf_hensel_lifting_vector(q, p, Mat(M[, 1]));
    [r, A] = qf_p_neighbor(q, p, X);
    Tr = mapget(T, qfrep(r, N));
    if(#Tr>1,
      for(i=1, #Tr, iso = qfisom(r, Tr[i]);if(0!=iso, qq = Tr[i]; break));
      ,
      qq = Tr[1]; iso = qfisom(r, Tr[1]);
    );
    col = mapget(M1, qq);
    if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
    A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
    spin = qf_spinor_norm(G[1][1], A1^-1);
    tpd[row, col1]+= nu(d, spin));
    for(i = 0, p-1,
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 2])%p);
      [r, A] = qf_p_neighbor(q, p, X);
      Tr = mapget(T, qfrep(r, N));
      if(#Tr>1,
        for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(0!=iso, qq = Tr[i]; break));
        ,
        qq = Tr[1]; iso = qfisom(r, Tr[1]);
      );
      col = mapget(M1, qq);
      if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
      A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
      spin = qf_spinor_norm(G[1][1], A1^-1);
      tpd[row, col1]+= nu(d, spin));
      X = qf_hensel_lifting_vector(q, p, Mat(i*M[, 1]+M[, 3])%p);
      [r, A] = qf_p_neighbor(q, p, X);
      Tr = mapget(T, qfrep(r, N));
      if(#Tr>1,
        for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(iso, qq = Tr[i]; break));
        ,
        qq = Tr[1]; iso = qfisom(r, Tr[1]);
      );
      col = mapget(M1, qq);
      if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
      A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
      spin = qf_spinor_norm(G[1][1], A1^-1);
      tpd[row, col1]+= nu(d, spin));
      for(j = 0, p-1,
        X = qf_hensel_lifting_vector(q, p, Mat((-i*j)*M[, 1]+i*M[, 2]+j*M[, 3]+M[, 4])%p);
        [r, A] = qf_p_neighbor(q, p, X);
        Tr = mapget(T, qfrep(r, N));
        if(#Tr>1,
          for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(iso, qq = Tr[i]; break));
          ,
          qq = Tr[1]; iso = qfisom(r, Tr[1]);
        );
        col = mapget(M1, qq);
        if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
        A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
        spin = qf_spinor_norm(G[1][1], A1^-1);
        tpd[row, col1]+= nu(d, spin));
        if(j != 0, X = (-qfeval(q, M[, 5])/2*M[, 1]+i*j*M[, 2]+j^2*M[, 4]+j*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [r, A] = qf_p_neighbor(q, p, X);
          Tr = mapget(T, qfrep(r, N));
          if(#Tr>1,
            for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(iso, qq = Tr[i]; break));
            ,
            qq = Tr[1]; iso = qfisom(r, Tr[1]);
          );
          col = mapget(M1, qq);
          if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
          A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
          spin = qf_spinor_norm(G[1][1], A1^-1);
          tpd[row, col1]+= nu(d, spin));
        );
        for(k=1, p-1,
          X = (i*k*M[, 1] + (-qfeval(q, M[, 5])/2-i*j)*M[, 2] + k^2*M[, 3] + j*k*M[, 4] + k*M[, 5])%p;
          X = qf_hensel_lifting_vector(q, p, Mat(X));
          [r, A] = qf_p_neighbor(q, p, X);
          Tr = mapget(T, qfrep(r, N));
          if(#Tr>1,
            for(i=1, #Tr, iso = qfisom(r, Tr[i]); if(iso, qq = Tr[i]; break));
            ,
            qq = Tr[1]; iso = qfisom(r, Tr[1]);
          );
          col = mapget(M1, qq);
          if(mapisdefined(minvariants, col), col1 = mapget(minvariants, col);
          A1 = G[invariants[row]][2]*A*iso^-1*G[col][2]^-1;
          spin = qf_spinor_norm(G[1][1], A1^-1);
          tpd[row, col1]+= nu(d, spin));
        );
      );
    );
  );
  tpd;


};

lexi(n) = {

  my(L, L1, i, j);
  if(n == 0, return(List()));
  if(n == 1, return(List([List([1]), List([2]), List([3]), List([4]), List([5])])));
  L = lexi(n-1);
  L1 = List();
  for(i = 1, 5, for(j = 1, #L, if(i <= L[j][1], Lj = L[j]; listinsert(Lj, i, 1); listput(L1, Lj))));
  L1;
}

basis_n(n) = {

  my(nn, v, L, pol);
  if(n == 0, return([1]));
  nn = abs(binomial(-5, n));
  v = vector(nn);
  L = lexi(n);
  L1 = [x1, x2, x3, x4, x5]~;
  for(i = 1, #L, pol = prod(j = 1, n, L1[L[i][j]]); v[i] = pol);
  v;
}

basis_n_A(A, n) = {

  my(nn, v, L, pol);
  if(n == 0, return([1]));
  nn = abs(binomial(-5, n));
  v = vector(nn);
  L = lexi(n);
  L1 = A*[x1, x2, x3, x4, x5]~;
  for(i = 1, #L, pol = prod(j = 1, n, L1[L[i][j]]); v[i] = pol);
  v;

}

pol_aux(p, n, i) = {

  my(v, j);
  if(i == 5, return([polcoef(p, n, x5)]));
  if(i == 4, v = []; for(j = 0, n, v = matconcat([v, pol_aux(polcoef(p, n-j, x4), j, 5)])));
  if(i == 3, v = []; for(j = 0, n, v = matconcat([v, pol_aux(polcoef(p, n-j, x3), j, 4)])));
  if(i == 2, v = []; for(j = 0, n, v = matconcat([v, pol_aux(polcoef(p, n-j, x2), j, 3)])));
  if(i == 1, v = []; for(j = 0, n, v = matconcat([v, pol_aux(polcoef(p, n-j, x1), j, 2)])));
  v;

}

pol_in_basis(p, n) = {

  return(pol_aux(p, n, 1));  

}

action_n(A, n) = {

  my(M, B);
  M = [];
  if(n == 1, return(A~));
  if(n == 2, return(action2(A)));
  if(n == 3, return(action3(A)));
  if(n == 4, return(action4(A)));
  B = basis_n_A(A, n);
  for(i = 1, #B, M = matconcat([M, pol_in_basis(B[i], n)~]));
  M;

}

delta(p) = {

  vars = [x1, x2, x3, x4, x5];
  sum(i = 1, 5, deriv(deriv(p, vars[i]), vars[i]));

}

delta_q(q, p) = {

  my(d, der);
  vars = [x1, x2, x3, x4, x5];
  d = 0;
  for(i = 1, 5, for(j = i, 5, der = deriv(deriv(p, vars[i]), vars[j])*q[i, j]; if(i == j, d+= der, d+=2*der)));
  return(d);

}

delta_basis(q, n) = {

  my(qa, L, v);
  qa = matadjoint(q);
  L = basis_n(n);
  v = [];
  for(i=1, #L, v = matconcat([v, pol_in_basis(delta_q(qa, L[i]), n-2)~]));
  v;

}

qf_invariant(q, n, d) = {

  my(MK, A, auto);
  MK = matker(delta_basis(q, n), 1);
  auto = qfauto(q)[2];
  for(i = 1, #auto, A = auto[i]; A*=matdet(A); MK = matintersect(MK, matker(action_n(A, n)-nu(qf_spinor_norm(q, A), d), 1)));
  MK*mathnf(MK~)~^-1;

}

qf_invariantab(q, a, b, d) = {

  my(auto, A, MK);
  auto = qfauto(q)[2];
  A = auto[1];
  A*= matdet(A);
  MK = matker(actionab(A, a, b)-nu(qf_spinor_norm(q, A), d), 1);
  for(i = 2, #auto, A = auto[i]; A*= matdet(A); MK = matintersect(MK, matker(actionab(A, a, b) - nu(qf_spinor_norm(q, A), d), 1)));
  MK*mathnf(MK~)~^-1;

};

psinverse(A) = (A*(A~*A)^-1)~;

matdecomp(A, B = matid(matsize(A)[1])) = {

  my(chC, V, C);
  C = psinverse(B)*A*B;
  chC = factor(charpoly(C));
  V = matrix(matsize(chC)[1], matsize(chC)[2]);
  for(i = 1, matsize(chC)[1], V[i, 1] = matker(subst(chC[i, 1], x, C)); V[i, 2] = chC[i, 2] == 1);
  [V, chC];

};

omf_decomp(omf) = {

  my(D, p, tp, L, todo, todop, dp, chC);
  D = matdet(omf[1][1][1])/2;
  p = 2;
  print(p);
  if(D%p == 0, tp = omf_tp1_singular(omf, p), tp = omf_tp1(omf, p));
  [dp, chC] = matdecomp(tp);
  L = listcreate();
  todo = listcreate();
  for(i=1, matsize(dp)[1], if(dp[i, 2], listput(L, [dp[i, 1], chC[i, 1]]), listput(todo, dp[i, 1])));
  while(#todo > 0,
   p = nextprime(p+1);
   print(p);
   if(D%p == 0, tp = omf_tp1_singular(omf, p), tp = omf_tp1(omf, p));
   todop = listcreate();
   for(i = 1, #todo,
     [dp, chC] = matdecomp(tp, todo[i]);
     for(j = 1, matsize(dp)[1], if(dp[j, 2], listput(L, [todo[i]*dp[j, 1], chC[j, 1]]), listput(todop, todo[i]*dp[j, 1])));
   );
   todo = todop;
  );
  L;  

};

omf_decompnd(omf, n, d) = {
 
  my(D, p, inv, tpinit, spins, tp, dp, L, todo, todop, chC);
  D = matdet(omf[1][1][1])/2;
  p = 2;
  inv = omf_invariantsnd(omf, n, d);
  if(D%p == 0, [tpinit, spins] = omf_tpnd_init_singular(omf, p, n), [tpinit, spins] = omf_tpnd_init(omf, p, n));
  tp = omf_tpnd(tpinit, d, inv, spins);
  [dp, chC] = matdecomp(tp);
  L = listcreate();
  todo = listcreate();
  for(i=1, matsize(dp)[1], if(dp[i, 2], listput(L, [dp[i, 1], chC[i, 1]]), listput(todo, dp[i, 1])));
  while(#todo > 0,
   p = nextprime(p+1);
   if(D%p == 0, [tpinit, spins] = omf_tpnd_init_singular(omf, p, n), [tpinit, spins] = omf_tpnd_init(omf, p, n));
   tp = omf_tpnd(tpinit, d, inv, spins);
   todop = listcreate();
   for(i = 1, #todo,
     [dp, chC] = matdecomp(tp, todo[i]);
     for(j = 1, matsize(dp)[1], if(dp[j, 2], listput(L, [todo[i]*dp[j, 1], chC[j, 1]]), listput(todop, todo[i]*dp[j, 1])));
   );
   todo = todop;
  );
  L;

};

omf_decompndmod(omf, n, d, l) = {

  my(D, p, inv, tpinit, spins, tp, dp, L, todo, todop, chC);
  D = matdet(omf[1][1][1])/2;
  p = 2;
  inv = omf_invariantsnd(omf, n, d);
  if(D%p == 0, [tpinit, spins] = omf_tpnd_init_singular(omf, p, n), [tpinit, spins] = omf_tpnd_init(omf, p, n));
  tp = Mod(omf_tpnd(tpinit, d, inv, spins), l);
  [dp, chC] = matdecomp(tp);
  L = listcreate();
  todo = listcreate();
  for(i=1, matsize(dp)[1], if(dp[i, 2], listput(L, [dp[i, 1], chC[i, 1]]), listput(todo, dp[i, 1])));
  while(#todo > 0,
   p = nextprime(p+1);
   if(D%p == 0, [tpinit, spins] = omf_tpnd_init_singular(omf, p, n), [tpinit, spins] = omf_tpnd_init(omf, p, n));
   tp = Mod(omf_tpnd(tpinit, d, inv, spins), l);
   todop = listcreate();
   for(i = 1, #todo,
     [dp, chC] = matdecomp(tp, todo[i]);
     for(j = 1, matsize(dp)[1], if(dp[j, 2], listput(L, [todo[i]*dp[j, 1], chC[j, 1]]), listput(todop, todo[i]*dp[j, 1])));
   );
   todo = todop;
  );
  L;

};



omf_decompns(omf, Bn) = {

  my(D, ds, invs, tpninits, spins, tpns, dones, todos, p, dp, todop, tpnss, dps, chCs, chC);
  D = matdet(omf[1][1][1])/2;
  ds = divisorssquarefree(D);
  invs = matrix(Bn+1, #ds, i, j, omf_invariantsnd(omf, i - 1, ds[j]));
  p = 2;
  if(D%p == 0, [tpninits, spins] = omf_tpnds_init_singular(omf, p, Bn), [tpninits, spins] = omf_tpnds_init(omf, p, Bn));
  tpns = matrix(Bn+1, #ds, i, j, omf_tpnd(tpninits[i], ds[j], invs[i, j], spins));
  tpnss = listcreate();
  listput(tpnss, tpns);
  dones = matrix(Bn+1, #ds, i, j, listcreate());
  todos = matrix(Bn+1, #ds, i, j, listcreate());
  dps = matrix(Bn+1, #ds);
  chCs = matrix(Bn+1, #ds);
  for(i=1, Bn+1, for(j=1, #ds, [dps[i, j], chCs[i, j]] = matdecomp(tpns[i, j])));
  \\dps = matrix(Bn+1, #ds, i, j, matdecomp(tpns[i, j]));
  for(i = 1, Bn+1,
    for(j = 1, #ds,
      for(k = 1, matsize(dps[i, j])[1], if(dps[i, j][k, 2], listput(dones[i, j], [dps[i, j][k, 1], chCs[i, j][k, 1]]), listput(todos[i, j], dps[i, j][k, 1])));
    );
  );
  while(sum(i = 1, matsize(todos)[1], sum(j = 1, matsize(todos)[2], #todos[i, j])) > 0,
    p = nextprime(p+1);
    if(D%p == 0, [tpninits, spins] = omf_tpnds_init_singular(omf, p, Bn), [tpninits, spins] = omf_tpnds_init(omf, p, Bn)); 
    tpns = matrix(Bn+1, #ds, i, j, omf_tpnd(tpninits[i], ds[j], invs[i, j], spins));
    listput(tpnss, tpns);
    for(i = 1, Bn+1,
      for(j = 1, #ds,
        todop = listcreate();
        for(k = 1, #todos[i, j],
          [dp, chC] = matdecomp(tpns[i, j], todos[i, j][k]);
          \\dp = matdecomp(tpns[i, j], todos[i, j][k]);
          for(l = 1, matsize(dp)[1], if(dp[l, 2], listput(dones[i, j], [todos[i, j][k]*dp[l, 1], chC[l, 1]]), listput(todop, todos[i, j][k]*dp[l, 1])));
        );
        todos[i, j] = todop;
      );
    );
  );
  [dones, tpnss];
  
};

matdecomprowfirst(tp, dec, row = 1) = {

  my(n, v, proy, m);
  n = matsize(tp)[1];
  m = matsize(dec)[2];
  for(i = row, n,
    v = vector(n, j, i == j);
    proy = matconcat(vector(m, j, (v*tp^(j-1))*dec)~)~;
    if(matdet(proy) != 0, return([i, proy]))
  );
  return(0);

};

matdecompisrow(tp, dec, row) = {

  my(n, v, m, proy);
  n = matsize(tp)[1];
  m = matsize(dec)[2];
  v = vector(n, j, row == j);
  proy = matconcat(vector(m, j, (v*tp^(j-1))*dec)~)~;
  if(matdet(proy) != 0, return([row, proy]), return(0));

};


matdecomprows(tp, dec) = {
 
  my(n, rows, v, proy, m);
  n = matsize(tp)[1];
  m = matsize(dec)[2];
  rows = listcreate();
  for(i = 1, n,
    v = vector(n, j, i == j);
    proy = matconcat(vector(m, j, (v*tp^(j-1))*dec)~)~;
    if(matdet(proy) != 0, listput(rows, [i]))
  );
  if(#rows == 0,
    for(i1 = 1, n,
      for(i2 = i1 + 1, n,
        v = vector(n, j, j == i1 + j == i2);
        proy = matconcat(vector(m, j, (v*tp^(j-1))*dec)~)~;
        if(matdet(proy) != 0, listput(rows, [i1, i2]))
      );
    );
  );
  rows;

};

invariantinverse(invariant) = {
 
  my(inv, i0);
  inv = Map();
  i0 = 0;
  for(i = 1, #invariant,
    for(j = 1, matsize(invariant[i][2])[2],
      mapput(inv, i0+j, i);
    );
    i0+= matsize(invariant[i][2])[2];
  );
  inv;

};

matdecomprowsinv(rows, inverse) = {

  my(rowsinv);
  rowsinv = listcreate();
  for(i = 1, #rows,
    listput(rowsinv, vector(#rows[i], j, mapget(inverse,rows[i][j])));
  );
  List(Set(rowsinv));

};