-- 
-- This is a collection of Macaulay2 functions written by Takehiko Yasuda. May 10, 2012. 
--


isBlowupSmooth = S ->( 
    -- checks if the Rees algebra of an ideal is smooth.
    all(affineCharts(S),isSmooth)  ) ;


isSmooth = R -> (
  -- checks if an affine ring is smooth
  if (isPolynomialRing R) then true else 
    (
      S := ambient R;
  	  I := ideal R;
      (ideal singularLocus I) == ideal(1_S) 
    )
  ) ;

      
isBlowupNormal = S ->( 
   -- checs if the Rees algebra of an ideal is normal.
      all(affineCharts(S),isNormal)) ;



affineCharts = S ->( 
     -- affine charts of a blowup with simplification of coordinate rings
      T := (flattenRing S)_0;
      varsOfS := apply(flatten entries vars S,i->sub(i,T));
      apply(varsOfS, i-> minimalPresentation(T/ideal(i - 1 )))) ;
      
      
affineCharts2 = S ->(
      -- affine charts of a blowup without simplification
      T := (flattenRing S)_0;
      U := ambient T;
      I := ideal T;
      varsOfS := apply(flatten entries vars S,i->sub(i,U));
      apply(varsOfS, i->  U / (I + ideal(i - 1)))) ;      


isFPure = R -> (
    -- checks if a ring is F-pure.
    p := char R;
    assert(p>0); 
    S := ambient R; 
    I := ideal R; 
    K:= ideal(apply(flatten entries vars S, i -> i^p));
    if isCI(R) then  (product flatten entries gens I)^(p-1)%K != 0
    else not isSubset( ideal(apply(flatten entries gens I, i -> i^p)) : I, K));
         
isCI = R -> ( 
     -- checks if the given expression of R is a complete intersection
     S := ambient R;
     I := ideal R;
     n := numgens S;
     r := numgens I;
     dim R == n-r );


frobeniusPushForward = (M,e) -> (  
      -- computes the pushforward F^e_*M of a module M by the e-iterate Frobenius map F^e.
      -- The base field must be a prime field, because the Frobenius map on the field is not the identity.  
      R := ring M;
      p := char R;
      assert(p>0);
      q := p^e; 
      I := ideal R;
      l := numgens I;
      B := gens ideal R;
      S := ambient R;                            
      n := numgens S;  
      qSequence := i -> apply(0..n-1, j -> (i % q ^ (n-j)) // q ^(n-j-1) );
      toNumber := i -> sum(n , j ->  i_j * q ^ (n-j-1)  );
      qQuotient := i -> apply(i, j -> j // q);
      qRemainder := i -> apply(i, j -> j % q);
      monoToMatrix := m ->  (coefficients m)_1_(0,0)* map(S^(q^n),S^(q^n),   
        (i,j) -> (e =( toList qSequence i ) + (exponents m)_0;  
        if (toNumber qRemainder e) == j then S_(toList qQuotient e) else 0));
      polyToMatrix := f -> 
         if f == 0_S then map(S^(q^n),S^(q^n),0_S) else sum(terms f, i -> monoToMatrix i);
      basisToMatrix := b -> 
         fold((i,j)->(i|j), apply(( flatten entries b), polyToMatrix ) );
      matrixToMatrix := m -> 
         fold((i,j)->(i||j), apply(apply( entries m , i ->matrix{i}),basisToMatrix));      
      ROverS := coker map(S^1,S^l, entries B); 
      PresenOverR := presentation minimalPresentation M;
      PresenOverS := presentation minimalPresentation(coker(sub(PresenOverR,S))**ROverS);
      L := matrixToMatrix PresenOverS;
      minimalPresentation coker sub(L,R)
      ) ;      
           

jetAmbient = (R,n) -> (
     -- produces the ambient ring which should become the ambient of the jet ring
     K := coefficientRing R;
     v := flatten entries vars R;
     v = apply(v, toString);
     w := for i to n list apply(v, x -> x | (toString i));
     w = flatten w;
     K[w]); 

jetIdeal = (R,n) -> (
     -- gets the defining ideal of the n-th jet scheme of Spec R
     r := numgens R;
     S := ambient R;
     jetAmb := jetAmbient (R,n);
     jetAmbVars := flatten entries vars jetAmb;
     jetAmbT :=jetAmb[t]/ideal(t^(n+1));
     polys := for i to (r-1) list sum( n+1, j -> jetAmbVars_( i + r*j )*t_jetAmbT^j ); -- polynomials like x0 + x1*t + x2 * t^2 + ... in jetAmbT
     f := map(jetAmbT, S, polys); -- map x0 |-> x0 + x1*t + x2 * t^2 + ...
     fIR := f ideal R; 
     fIRGens := flatten entries gens fIR ;
     jetIdealGens := for g in fIRGens list flatten entries (coefficients(g))_1 ;
     jetIdealGensList := flatten jetIdealGens ;
     jetIdeal := ideal jetIdealGensList;
     sub(jetIdeal, jetAmb)
     );
	       
	       

jetRing = (R,n) -> (
     -- gets the coordinate ring of the n-th jet scheme of Spec R
     I:= jetIdeal(R,n);
     quotient I
     );



villamayorIdeal= M ->(
     -- computes the center ideal of the blowup at a module M, following [Villamayor, On flattening of ...., J. Alg.]
        r := rank M;
        P := presentation M;
        s := rank source P;
        t := rank target P;
        I := {}; for j to s-1 when #I < t-r do (J := append(I,j); 
               if rank coker P_J  == t- #J then I = J );
        fittingIdeal(r,coker P_I)) ;

universalFlattening = M -> (
     -- computes the blowup at a module M
        I := villamayorIdeal M;
        reesAlgebra I);  


affineMonomialCurve = (L,K) -> ( 
     -- get the ring K[t^n_1, ...] of the affine monomial curve associated to the list L={n_1,...}
      l := #L; 
      f := map(K[z],K[vars(0..l-1),Degrees=>L],apply(L,n->z^n));
      coimage f);
 

modulePn = (R,n) ->(  
     -- computes the module Pn for R = affineMonomialCurve(?)
        l:= numgens R;
        L:= flatten degrees R;
        R1:=newRing(R,Variables=>vars(l..(2*l-1)));
        S := tensor(R,R1,Degrees=>(L|L));
        T := ambient S;    
        I = ideal(apply(0..l-1, i-> T_i -T_(i+l)));
        U := T/(ideal S + I^(n+1));
        f := map(U,R,(gens U)_(toList (0..l-1)));
        pushForward(f,U^1)) ;



nashBlowup = (R,n) -> ( 
     -- computes the n-th Nash blowup of affineMonomialCurve(??)
        P := modulePn(R,n);
        I := villamayorIdeal P;
        reesAlgebra I);



fBlowup = (R,e) -> ( 
     -- compute the e-th F-blowup
       reesAlgebra fBlCenter(R,e));



fBlSpecial = (R,e) -> ( 
     -- computes the special fiber of the e-th F-blowup
       specialFiber fBlCenter(R,e));



fBlCenter = (R,e) -> (
   -- computes the center ideal of the e-th F-blowup
     villamayorIdeal frobeniusPushForward(R^1,e));



pCyc3D = p -> (
     -- computes the p-cyclic 3-dimensional quotient singularity V_3 / C_p = Spec k[V_3]^(C_p) 
     S := ZZ/p [x,y,z];
     f := map(S,S,{x,-x+y,x-y+z});
     Ny := y * product(1..(p-1), i -> f^i(y));
     Nz := z * product(1..(p-1), i -> f^i(z));
     d :=  y^2 -2*x*z -x*y;
     T := ZZ/p [X,Y,Z,W];
     F := map(S,T,{x,Ny,Nz,d});
     coimage F);