////////////////////////////////////////////////////////////////////////
//                     USES pol_power!!!!!!!!                         //
////////////////////////////////////////////////////////////////////////

/*

Fast functions to preform computation with Witt vectors of
Length 3.

See comp.tex

*/

load 'pol_power.m';


// Coef's are some divisions of binomial coefficients that
// we need to have explicitly as integer (or in GF(p))
Coef1 := function(n,p,i)
    return Binomial(n*p,i) div p;
end function;

Coef := function(p,i)
    return Coef1(1,p,i);
end function;


Coef2 := function(n,p,i)
    return (Binomial(n,i)-Binomial(n*p,i*p)) div p;
end function;


// ////////////////////////////////////////////////////////////

// psi_1 (\sfct from comp.tex)
SpcFct := function(f)

    v:=Terms(f);
    n:=#v;
    
    if n le 1 then
        return 0;
    end if;

    P:=Parent(f);
    p:=Characteristic(P);

    // split in 'half-sized' pols.
    g1:=&+[ v[i] : i in [1..(n div 2)] ];
    g2:=&+[ v[i] : i in [((n div 2)+1)..n] ];

    delete v;
    
    tmp:= &+[ -Coef(p,i)*g1^i*g2^(p-i) : i in [1..(p-1)] ];

        
    return $$(g1) + $$(g2) + tmp;
    
end function;


// eta_1(x,y) (\afctt from comp.tex)
fSpcFct := function(x,y)

    if (x eq 0) or (y eq 0) then
        return 0;
    else
	p:=Characteristic(Parent(x));
	return  &+[ -Coef(p,i)*x^i*y^(p-i) : i in [1..(p-1)] ];
    end if;
	
end function;


/*

TESTED!!!

RR<X,Y,Z,W>:=PolynomialRing(Integers(),4);

p:=17;R<x,y,z,w>:=PolynomialRing(GF(p),4);

p1:=R!(((X^p+Y^p+Z^p+W^p)-(X+Y+Z+W)^p) div p);
q1:=SpcFct(x+y+z+w);

p1 eq q1;


*/


// eta_1 (\afctt from comp.tex)
SpcFctv := function(v)

    r:=#v;
    
    if (r eq 1) or (r eq 0) then
        return 0;
    end if;

    P:=Parent(v[1]);
    p:=Characteristic(P);

    // split in 'half-sized' pols.
    v1:=[ v[i] : i in [1..(r div 2)] ];
    v2:=[ v[i] : i in [((r div 2)+1)..r] ];

    return fSpcFct(&+v1,&+v2) + $$(v1) + $$(v2);

end function;

// ////////////////////////////////////////////////////////////


// second coordinate of an integer as a Witt vector

NextCoord := function(a,p)
    return (a^p - a) div p;
end function;

// more coeffs.

Coef3:= function(p,i)
    return NextCoord(Binomial(p,i) div p,p);
end function;


Coef4 := function(p,i)
    return (Binomial(p^2,i) div p^2);
end function;


// The number of terms is relevant on how to compute this function.
// Picked a random number of terms (namely, 5) to divide.
// ************************************************************
// SHOULD CHECK WHEN ONE METHOD BECOMES FASTER THAN THE OTHER!!!
// ************************************************************


// split in 2 of same size
Spc4Fct1 := function(f)
 
    w:=Terms(f);
    n:=#w;
   
    if n le 1 then
        return 0;
    end if;
    
    
    P:=Parent(f);
    p:=Characteristic(P);

    g1:=&+[ w[i] : i in [1..(n div 2)] ];
    g2:=&+[ w[i] : i in [((n div 2)+1)..n] ];

    delete w;

    // precompute powers of g1 and g2 to use with pol_power2
    pwrs1:=[ 1, g1 ];
    for i in [3..p] do
        pwrs1[i]:=pwrs1[i-1]*g1;
    end for;

    pwrs2:=[ 1, g2 ];
    for i in [3..p] do
        pwrs2[i]:=pwrs2[i-1]*g2;
    end for;
    
    // res1
    res:= -&+[ Coef4(p,i)*pol_power2(g1,i : pwrs:=pwrs1)*pol_power2(g2,(p^2-i) : pwrs:=pwrs2) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    // res2
    res+:= &+[ Coef3(p,i)*pol_p_power(pwrs1[i+1],p,1)*pol_p_power(pwrs2[p-i+1],p,1) : i in [1..(p-1)] ];
    
    v:=[ -Coef(p,i)*pwrs1[i+1]*pwrs2[p-i+1] : i in [1..(p-1)] ];

    // res3
    res+:= SpcFctv(v);

    // sf:=fSpcFct(g1,g2);
    sf := &+v;

    delete v;

    sg1:=SpcFct(g1);
    
    sg2:=SpcFct(g2);

    // res5
    res+:=fSpcFct(sg1+sg2,sf);

    delete sf;
    
    // res4
    res+:=fSpcFct(sg1,sg2);

    delete sg1, sg2;
    
    return $$(g1) + $$(g2) + res;

end function;



// split in monomial + rest
Spc4Fct2 := function(f)
 
    w:=Terms(f);
    n:=#w;
   
    if n le 1 then
        return 0;
    end if;
    
    
    P:=Parent(f);
    p:=Characteristic(P);

    g1:=w[1];
    g2:=f-g1;

    delete w;
    // precompute powers of g1 and g2 to use with pol_power2
    pwrs1:=[ 1, g1 ];
    for i in [3..p] do
        pwrs1[i]:=pwrs1[i-1]*g1;
    end for;

    pwrs2:=[ 1, g2 ];
    for i in [3..p] do
        pwrs2[i]:=pwrs2[i-1]*g2;
    end for;
    
    // res1
    res:= -&+[ Coef4(p,i)*pol_power2(g1,i : pwrs:=pwrs1)*pol_power2(g2,(p^2-i) : pwrs:=pwrs2) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    // res2
    res+:= &+[ Coef3(p,i)*pol_p_power(pwrs1[i+1],p,1)*pol_p_power(pwrs2[p-i+1],p,1) : i in [1..(p-1)] ];
    
    v:=[ -Coef(p,i)*pwrs1[i+1]*pwrs2[p-i+1] : i in [1..(p-1)] ];

    // res3
    res+:= SpcFctv(v);

    // sf:=fSpcFct(g1,g2);
    sf := &+v;

    delete v;

    // zero in this case
    // sg1:=SpcFct(g1);
    
    sg2:=SpcFct(g2);

    // res5
    res+:=fSpcFct(sg2,sf);

    delete sf, sg2;
    
    return $$(g2) + res;

end function;



    
// psi_2 (\ssfct in comp.tex) := Spc2Fct + Spc3Fct (but simplified)
Spc4Fct := function(f)
    // decides which function to call
    
    v:=Terms(f);
    n:=#v;

    if n gt 5 then
        // larger number of terms
        return Spc4Fct1(f);
    else
        // smaller number of terms
        return Spc4Fct2(f);
    end if;
        
end function;


/*

TESTED!!!

Testing Spc4Fct and Spc4Fct2

p:=5;F<a>:=GF(p^3);P<x0,y0>:=PolynomialRing(F,2);
k:=7;deg:=5;
f:=&+[ Random(F)*x0^Random([0..deg])*y0^Random([0..deg]) : i in [1..k]];

time f1:=Spc4Fct(f);
time f2:=Spc4Fct2(f);

f1 eq f2;

*/




/*


Divp := function(f,p)

  ZZ:=Integers();
  res:=0;

  for m in Monomials(f) do
      coef:=MonomialCoefficient(f,m);
      res+:=(ZZ!coef div p)*m;
  end for;

  return res;

end function;

Div2p := function(f,p)

  P:=Parent(f);
  R:=PolynomialRing(Integers(),4);

  return P!(R!f div p);

end function;


p:=11;R<x,y,z,w>:=PolynomialRing(GF(p),4);

time p1:=Spc4Fct(x+y+z+w);
time q1:=Spc2Fct(x+y+z+w) + Spc3Fct(x+y+z+w);

p1 eq q1;

RR<X,Y,Z,W>:=PolynomialRing(Integers(),4);

time q1:=R!(((((X^(p^2)+Y^(p^2)+Z^(p^2)+W^(p^2))-(X+Y+Z+W)^(p^2)) div p)
-(((X^p+Y^p+Z^p+W^p)-(X+Y+Z+W)^p) div p)^p) div p);

p1 eq q1;

Rp3:=ResidueClassRing(p^3);
RRR<XX,YY,ZZ,WW>:=PolynomialRing(Rp3,4);

time q1:=R!(Div2p(Div2p((XX^(p^2)+YY^(p^2)+ZZ^(p^2)+WW^(p^2))-(XX+YY+ZZ+WW)^(p^2),p) -
  (Div2p((XX^(p)+YY^(p)+ZZ^(p)+WW^(p))-(XX+YY+ZZ+WW)^(p),p))^p,p));

p1 eq q1;


*/




// eta_2 (\aafct in comp.tex)
fSpc4Fct := function(x,y)

    P:=Parent(x);
    p:=Characteristic(P);

    // precompute powers of g1 and g2 to use with pol_power2
    pwrs1:=[ 1, x ];
    for i in [3..p] do
        pwrs1[i]:=pwrs1[i-1]*x;
    end for;

    pwrs2:=[ 1, y ];
    for i in [3..p] do
        pwrs2[i]:=pwrs2[i-1]*y;
    end for;
    
    // res1
    res:= -&+[ Coef4(p,i)*pol_power2(x,i : pwrs:=pwrs1)*pol_power2(y,(p^2-i) : pwrs:=pwrs2) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    // res2
    res+:= &+[ Coef3(p,i)*pol_p_power(pwrs1[i+1],p,1)*pol_p_power(pwrs2[p-i+1],p,1) : i in [1..(p-1)] ];
    
    v:=[ -Coef(p,i)*pwrs1[i+1]*pwrs2[p-i+1] : i in [1..(p-1)] ];

    return res + SpcFctv(v);

end function;

// TESTED!!



/*

TESTED!!!!

Now test second coord.

load 'witt_gen.m';

p:=31;R<x0,x1,x2,y0,y1,y2>:=PolynomialRing(GF(p),6);

time p1:=R!WittSumPol(2,p)[3];

time q1:=x2+y2+fSpcFct(x1,y1)+fSpcFct(x1+y1,SpcFct(x0+y0))+Spc3Fct(x0+y0)+Spc2Fct(x0+y0);

p1 eq q1;

*/


/*

Now, test the eq. of the elliptic curve

load 'witt_gen.m';


p:=11;P<a0,a1,a2,b0,b1,b2,x0,x1,x2>:=PolynomialRing(GF(p),9);
PP<A0,A1,A2,B0,B1,B2,X0,X1,X2>:=PolynomialRing(Integers(),9);


x:=[x0,x1,x2];
a:=[a0,a1,a2];
b:=[b0,b1,b2];

f:=x0^3+a0*x0+b0;
sf:=SpcFct(f);

rhs1:=WittPwr(x,3);
rhs1:=WittSum(rhs1,WittProd(a,x));
rhs1:=WittSum(rhs1,b);

rhs1[2] eq (3*x0^(2*p)+a0^p)*x1 + a1*x0^p + b1 + sf;

rhs1:=rhs1[3];


ev:=[ 3*x0^(2*p)*x1 , a0^p*x1, a1*x0^p , b1]; 

term1:= (3*x0^(2*p^2)+a0^(p^2))*x2 + 3*x0^(p^2)*x1^(2*p) +
    (-NextCoord(3,p)*x0^(2*p^2)+a1^p)*x1^p + a2*x0^(p^2) + b2;

term2:= SpcFctv(ev);

term3:= fSpcFct( (3*x0^(2*p)+a0^p)*x1 + a1*x0^p + b1 , sf );

term4:= Spc2Fct(f);

term5:= Spc3Fct(f);

rhs2:=term1+term2+term3+term4+term5;

rhs2 eq rhs1;


Term1:=(3*X0^(2*p^2)+A0^(p^2))*X2 + 3*X0^(p^2)*X1^(2*p) +
    (-NextCoord(3,p)*X0^(2*p^2)+A1^p)*X1^p + A2*X0^(p^2) + B2;

Term2:=(A0^(p^2)*X1^p + A1^p*X0^(p^2) + B1^p + (3*X0^(2*p)*X1)^p -
    (A0^p*X1 + A1*X0^p + B1 + 3*X0^(2*p)*X1 +
        ((X0^(3*p) + A0^p*X0^p + B0^p - (X0^3 + A0*X0 + B0)^p) div p))^p);

Term3:=(X0^(3*p^2) + A0^(p^2)*X0^(p^2) + B0^(p^2) - (X0^3 + A0*X0 + B0)^(p^2));

rhs3:= Term1 + (( p* Term2 + Term3 ) div p^2);

rhs3:= P!rhs3;

rhs3 eq rhs1;

*/

