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);
    term:=Terms(f)[1];

    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)

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

    term:=v[1];
    P:=Parent(term);
    p:=Characteristic(P);
    w:=Remove(v,1);

    g:=&+w;

    return fSpcFct(term,g) + $$(w);

end function;


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


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;


// THIS IS WORKING BUT IS SLOWER!!!  CHECK THE OLD ONE OUT!!
//
// It's slower if not many terms, but get faster with more terms.
//
// psi_2 (\ssfct in comp.tex) = Spc2Fct + Spc3Fct (but simplified)

// This is the new version.  The old one is below!
Spc4Fct := function(f)

    v:=Terms(f);
    n:=#v;

    if n le 1 then
	return 0;
    end if;
	    
    
    P:=Parent(f);
    p:=Characteristic(P);

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

    delete v;
    

    res1:= -&+[ Coef4(p,i)*g1^i*g2^(p^2-i) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    res2:= &+[ Coef3(p,i)*g1^(p*i)*g2^(p*(p-i)) : i in [1..(p-1)] ];
    
    v:=[ -Coef(p,i)*g1^i*g2^(p-i) : i in [1..(p-1)] ];
    
    res3:= SpcFctv(v);

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

    sg1:=SpcFct(g1);
    
    sg2:=SpcFct(g2);
    
    res4:=fSpcFct(sg1,sg2);

    res5:=fSpcFct(sg1+sg2,sf);
    
    return $$(g1) + $$(g2) + res1 + res2 + res3 + res4 + res5;

end function;


// orginal function!
Spc4Fct2 := function(f)

    v:=Terms(f);
    n:=#v;

    if n le 1 then
	return 0;
    end if;
	    
    
    P:=Parent(f);
    p:=Characteristic(P);

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

    delete v;
    

    res1:= -&+[ Coef4(p,i)*g1^i*g2^(p^2-i) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    res2:= &+[ Coef3(p,i)*g1^(p*i)*g2^(p*(p-i)) : i in [1..(p-1)] ];
    
    v:=[ -Coef(p,i)*g1^i*g2^(p-i) : i in [1..(p-1)] ];
    
    res3:= SpcFctv(v);

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

    sg:=SpcFct(g2);
    
    res4:=fSpcFct(SpcFct(g1),sf+sg);

    res5:=fSpcFct(sg,sf);
    
    return $$(g1) + $$(g2) + res1 + res2 + res3 + res4 + res5;

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);

    res1:= -&+[ Coef4(p,i)*x^i*y^(p^2-i) : i in [1..(p^2-1)] | (i mod p) ne 0 ]; 

    res2:= &+[ Coef3(p,i)*x^(p*i)*y^(p*(p-i)) : i in [1..(p-1)] ];
    
    v:=[ -Coef1(1,p,i)*x^i*y^(p-i) : i in [1..(p-1)] ];
    
    res3:= SpcFctv(v);

    return res1 + res2 + res3;

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;

*/

