
// Local power series expansions on the Jacobian

intrinsic LocalPowerSeries(K::Fld : Precision := 7) -> SeqEnum
    {}
    F<f0,f1,f2,f3,f4,f5,f6> := FunctionField(K,7 : Global); 
    return LocalPowerSeries([f0,f1,f2,f3,f4,f5,f6] : Precision := Precision);
end intrinsic;


// RngFrmElt is the dummy type for RngMSerElt
intrinsic LocalPowerSeries(ff::[RngElt] : Precision := 7) -> RngFrmElt 
    {}
    _ , JacRels := JacobianRelations(ff);
    
    require #ff eq 7 : "Argument must have length 7.";
    F := Universe(ff);
    // We pick up the global object:
    P<X00,X01,X02,X03,X04,X05,X06,X07,
      X08,X09,X10,X11,X12,X13,X14,X15> := PolynomialRing(F,16 : Global);

    // We first need to find the 13 relations that express the
    // s_i in terms of other s_j. These are of the form X00Xij - ...
    
    // New polynomial ring "near the origin":
    R<S01, S02, S03, S04, S05, S06, S07, 
        S08, S09, S10, S11, S12, S13, S14, S15> := PolynomialRing(F, 15);

    // The dehomogenisation homomorphism:
    dehomog := hom < P -> R | 1, S01, S02, S03, S04, S05, S06, S07, 
        S08, S09, S10, S11, S12, S13, S14, S15>;
    
    // We already know S01 and S02. We want to find expressions for the others:
    sRel := []; sRel[1] := S01; sRel[2] := S02;
    CurRel := 1;
    for sIndex in [3..15] do
        errChk := CurRel;
        while HomogeneousComponent(dehomog(JacRels[CurRel]), 1) ne R.sIndex do
            CurRel := CurRel ne #JacRels select (CurRel + 1) else 1;
            if CurRel eq errChk then // We've cycled through without finding our Si
                require true: "Error: couldn't find an expression for Si in the relations.";
            end if;
        end while;
        sRel[sIndex] := -HomogeneousComponent(dehomog(JacRels[CurRel]),2);
    end for;
    
    // We now need to recursively generate the power series.
    // We know that we should reach the required precision in Precision steps.
    // errChk is to make sure that we're not in an infinite loop.
    
    Coef<S1,S2> := PolynomialRing(F, 2);
    R<X> := PowerSeriesRing(Coef);
    s := [BigO(X^2) : i in [1..15]];
    s[1] := X*S1; s[2] := X*S2;
    
    errChk := 1;
    repeat
        s := [Evaluate(sRel[i], s) : i in [1..15]];
        errChk +:= 1;
        if errChk gt Precision then
            require true: "Error: couln't reach required precision.";
        end if;
    until Minimum([AbsolutePrecision(s[i]) : i in [3..15]]) ge Precision;
    
    return s;
    
end intrinsic;


intrinsic FindUpperCoordinates(K::Fld : Precision := 8) -> SeqEnum
    {Returns the last three projective coordinates (ie u13, u14 and u15)
        when (1:s1:s2:..) is added to (1:t1:t2:..) on the Jacobian
        corresponding to the curve y^2 = f0 + f1*X + .. + f6*X^6 over K.}
    F<f0,f1,f2,f3,f4,f5,f6> := FunctionField(K,7 : Global);
    a, b := FindUpperCoordinates([f0,f1,f2,f3,f4,f5,f6] : Prec := Precision);
    return a, b;
end intrinsic;

// FIND UPPER COORDINATES

// If c = a + b on the Jacobian this intrinsic calculates approximations 
// for c13, c14 and c15 as power series in K[[s1,s2,t1,t2]];
intrinsic FindUpperCoordinates(ff::[RngElt] : Prec := 8) -> SeqEnum
    {}
    
    // Precision less than 6 gives errors:
    Prec := Maximum(Prec, 6);
    
    // Almost nothing works unless we're over a field!:
    require IsField(Universe(ff)): "The elements of Argument 1 must belong to a field";
    
    require #ff eq 7: "Argument 1 must have length 7";
    
    // The field that we are working over:
    F := Universe(ff);
    
    f0, f1, f2, f3, f4, f5, f6 := Explode(ff);
    
    T<s01, s02, s03, s04, s05, s06, s07, s08, s09, s10, s11, s12, s13, s14, s15, 
    t01, t02, t03, t04, t05, t06, t07, t08, t09, t10, t11, t12, t13, t14, t15> := PolynomialRing(F, 30) ;
    
    // These functions give the cubic that meets the curve at the two sets of points
    // given corresponding to (a01,..,a15) and (b01,..,b15). (See Flynn's article)
    a := [];
    // mu:
    mu := s10*t15 + t10*s15 - s14*t11 - t14*s11 + s13*(t12 + t13) + t13*(s12+ s13);
    // alpha:
    alpha := s07*t15 + t07*s15 + s13*t09 + t13*s09 - s14*t08 - t14*s08;
    // beta:
    beta := -s06*t15 - t06*s15 + s08*(t12 + t13) + t08*(s12 + s13) - s11*t09 - t11*s09;
    // gamma
    gamma := s06*t14 + t06*s14 - s07*(t12 + t13) - t07*(s12 + s13) + s09*t10 + t09*s10;
    // delta
    delta := -s06*t13 - t06*s13 + s07*t11 + t07*s11 - s08*t10 - t08*s10;
    
    // These formulas are lifted straight from Flynn's article:
    c13 := s15^2*t15^2*(delta^2 - f0*mu^2);
    c14 := s13*t13*(s15*t15*(f5*mu^2 - 2*alpha*beta) - (s15*t14 + s14*t15)*(alpha^2 - f6*mu^2));
    c15 := s13*t13*s15*t15*(alpha^2 - f6*mu^2);
    
    PowSer := LocalPowerSeries(ff : Precision:=Prec);
    
    Coeffs<S1,S2,T1,T2> := PolynomialRing(F,4);
    P<X> := PowerSeriesRing(Coeffs);
    
    Q := Universe(PowSer);
    B := BaseRing(Q);
    f := hom< B -> Coeffs | [S1, S2] >;
    g := hom< B -> Coeffs | [T1, T2] >;
    
    s1 := S1*X; s2 := S2*X; t1 := T1*X; t2 := T2*X;
    
    aLocalSeries := [];
    for i in [1..#PowSer] do
        seq, val := Eltseq(PowSer[i]);
        p := RelativePrecision(PowSer[i]);
        aLocalSeries[i] := elt< P | val, [f(seq[j]) : j in [1..#seq]], p>;
    end for;
    
    bLocalSeries := [];
    for i in [1..#PowSer] do
        seq, val := Eltseq(PowSer[i]);
        p := RelativePrecision(PowSer[i]);
        bLocalSeries[i] := elt< P | val, [g(seq[j]) : j in [1..#seq]], p>;
    end for;

    LocalSeries := aLocalSeries cat bLocalSeries;
    
    // We now evaluate our exact projective expressions for c13, c14 and c15
    // at the series approximations for s1,..,s15 and t1,..,t15:
    cp13 := Evaluate(c13, LocalSeries);
    cp14 := Evaluate(c14, LocalSeries);
    cp15 := Evaluate(c15, LocalSeries);
    
    // We take out any common factors to make things nice:
    seq13 := Eltseq(cp13); seq14 := Eltseq(cp14); seq15 := Eltseq(cp15);
    DivRenorm := GreatestCommonDivisor(seq13 cat seq14 cat seq15);
    DivRenorm := DivRenorm*X^(TotalDegree(DivRenorm));
    cp13 := P!(cp13/DivRenorm); cp14 := P!(cp14/DivRenorm); cp15 := P!(cp15/DivRenorm);

    // We know that, after some rescalling we should have:
    // cp13 ~= (s1+t1)^2*(s2+t2)^2
    // cp14 ~= 2*(s1+t1)*(s2+t2)^3
    // cp14 ~= (s2+t2)^4
    // Here we determine what rescalling is necessary:
    
    LeadingTerm := Eltseq(cp13)[1];
    
    ProjectiveFactor := Numerator(LeadingTerm / ((S1+T1)^2*(S2+T2)^2));
    MultRenorm := Denominator(LeadingTerm / ((S1+T1)^2*(S2+T2)^2));
    
    MultRenorm := MultRenorm*X^TotalDegree(MultRenorm);
    ProjectiveFactor := ProjectiveFactor*X^TotalDegree(ProjectiveFactor);
    
    cp13 *:= MultRenorm; cp14 *:= MultRenorm; cp15 *:= MultRenorm;
    
    return ProjectiveFactor, [cp13, cp14, cp15];
    
end intrinsic;

// PROCEDURE: LinearSolve()
// Given a relation, an index, i to solve for (ie if Index = 5 then we
// are solving for X05) and the power series expansions so far this 
// procedure solves for Xi and then clears any denominators.

procedure LinearSolve(Rel, Index, ~Exp)
    // The polynomial ring in X01, X02, X03 etc:
    M := Parent(Rel);
    P<X> := Universe(Exp);
    // We need the power series ring over the function field:
    Coeffs := BaseRing(P);
    BiggerCoeffs := FieldOfFractions(Coeffs);
    Q := PowerSeriesRing(BiggerCoeffs);
    
    top, bottom := Explode(Terms(Rel, M.(Index+1)));
    bottom := ExactQuotient(bottom, M.(Index+1));
    print top, bottom;
    top := Evaluate(top, Exp); 
    bottom := Evaluate(bottom, Exp);
    
    print top, bottom;
    
    printf "Valuation of top: %o, Valuation of bottom: %o\n", Valuation(top), Valuation(bottom);
    printf "RelPrec of top: %o ..of bottom: %o\n", RelativePrecision(top), RelativePrecision(bottom);
    
    result := Q!(-top / bottom);
    
    // We must clear any denominators:
    seq := Eltseq(result);
    denom := &*[Denominator(seq[i]) : i in [1..#seq]];
    denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]]));
    denom := denom*X^TotalDegree(denom);
    // Update all the other coordinates:
    for i in [1..16] do
        Exp[i] *:= denom;
    end for;
        
    Exp[Index+1] := P!(denom*result);
end procedure;

// INTRINSIC: UpdateCoords
intrinsic UpdateCoords(ff::[RngElt], AlreadyKnown::[RngIntElt], Wanted::[RngIntElt], ~SeriesExpansions::SeqEnum)
    {}

    SolvingSequence := JacobianSolvingSequence(ff, AlreadyKnown, Wanted);
    
    // The polynomial ring in X01, X02, X03 etc:
    M := Parent(SolvingSequence[1][2]);
    
    P<X> := Universe(SeriesExpansions);
    
    // We need the power series ring over the function field:
    Coeffs := BaseRing(P);
    BiggerCoeffs := FieldOfFractions(Coeffs);
    Q := PowerSeriesRing(BiggerCoeffs);
    
    for Step in SolvingSequence do
    
        CurrentIndex, Relation, IsSquare := Explode(Step);

        printf "Attempting to calculate s%o using the following relation:\n%o\n", CurrentIndex, Relation;

        if not IsSquare then
            top, bottom := Explode(Terms(Relation, M.(CurrentIndex+1)));
            bottom := ExactQuotient(bottom, M.(CurrentIndex+1));

            top := Evaluate(top, SeriesExpansions); bottom := Evaluate(bottom, SeriesExpansions);

            result := Q!(-top / bottom);
        
        else
            printf "(Note: We must evaluate a square root)\n";
            c, b, a := Explode(Terms(Relation, M.(CurrentIndex+1)));
            require b eq 0: "Unexpected quadratic! (rather than just a square root)";
            
            a := ExactQuotient(a, (M.(CurrentIndex+1))^2);
            a := Evaluate(a, SeriesExpansions); c := Evaluate(c, SeriesExpansions);
            
            // We need to write a SquareRoot function as this generic
            // intrinsic doesn't working in characteristic not equal to 0.

            result := SquareRoot(Q!(-c/a));
        
        end if;
        // We must clear any denominators:
        seq := Eltseq(result);
        denom := &*[Denominator(seq[i]) : i in [1..#seq]];
        // We want the lcm:
        if #seq ne 1 then 
            denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]]));
        end if;
        denom := denom*X^TotalDegree(denom);
        // Update all the other coordinates:
        for i in {1..16} do
            SeriesExpansions[i] *:= denom;
        end for;
            
        SeriesExpansions[CurrentIndex+1] := P!(denom*result);

        printf "Finished (Relative precision: %o. Denominator: %o).\n\n", 
            RelativePrecision(SeriesExpansions[CurrentIndex+1]), Eltseq(denom)[1];
            
    end for;

end intrinsic;

intrinsic FormalGroup(k::Fld : Precision:=7) -> SeqEnum
    {}
    k<f0,f1,f2,f3,f4,f5,f6> := FunctionField(k, 7);
    return FormalGroup([f0,f1,f2,f3,f4,f5,f6] : Precision:= Precision);
end intrinsic;

intrinsic FormalGroup(ff::[RngElt] : Precision := 7) -> SeqEnum
    {}
    
    _, Jac := JacobianRelations(ff);
    M<X00,X01,X02,X03,X04,X05,X06,X07,X08,X09,X10,X11,X12,X13,X14,X15> := Universe(Jac);
    
    f0, f1, f2, f3, f4, f5, f6 := Explode(ff);
    
    Expansion := [];
    
    // First find x13, x14 and x15:
    ProjectiveFactor, b := FindUpperCoordinates(ff : Prec := Precision);
    
    P<X> := Parent(b[1]);
    
    Coeffs<S1,S2,T1,T2> := BaseRing(P);
    s1 := S1*X; s2 := S2*X; t1 := T1*X; t2 := T2*X;
    BiggerCoeffs<u1, u2, v1, v2> := FieldOfFractions(Coeffs);
    Q<Y> := PowerSeriesRing(BiggerCoeffs);
    
    Expansion := [ P!1 : i in [1..16]];
    
    Expansion[13 + 1] := b[1];
    Expansion[14 + 1] := b[2];
    Expansion[15 + 1] := b[3];
    
    
    // We know that x15/x0 ~= (s2 + t2)^4.
    // Since x15 ~= (s2 + t2)^4*ProjectiveFactor we know that
    // x0 ~= ProjectiveFactor.   
    Expansion[0 + 1] := ProjectiveFactor + BigO(X^(Valuation(ProjectiveFactor) + 2));
    
    // We can find x10, x11 and x12 without any need of recursion:
    UpdateCoords( ff, [13, 14, 15], [12, 11, 10], ~Expansion);

    printf "Attempting to calculate s5 using a very ugly relation: (This may take some time!)\n";

    // The horrible relation that we use to get x5:
    rel :=   X05^2*X10*X13 - 1/4*X05^2*X11^2 + f6*X05*X10^2*X13 + 1/2*f5*X05*X10*X11*X13 
        + f4*X05*X10*X13^2 + 1/2*f3*X05*X11*X13^2 + 1/2*f1*X05*X11*X13*X15 + 
        f2*X05*X13^3 + f0*X05*X13^2*X15 + (f4*f6 - 1/4*f5^2)*X10^2*X13^2 + 
        f3*f6*X10*X11*X13^2 + 1/2*f3*f5*X10*X13^3 - 1/2*f1*f5*X10*X13^2*X15 + 
        f2*f6*X11^2*X13^2 + f1*f6*X11*X12*X13^2 + (f1*f6 + f2*f5)*X11*X13^3 + 
        f0*f6*X12^2*X13^2 + (2*f0*f6 + f1*f5)*X12*X13^3 + f0*f5*X12*X13^2*X14 + 
        (f0*f6 + 2*f1*f5 + f2*f4 - 1/4*f3^2)*X13^4 + (f0*f5 + f1*f4)*X13^3*X14 +
        1/2*f1*f3*X13^3*X15 + f0*f4*X13^2*X14^2 + f0*f3*X13^2*X14*X15 + (f0*f2 -
        1/4*f1^2)*X13^2*X15^2;
    
    c, b, a := Explode(Coefficients(rel, M.6));
    
    printf "Expanding a...";
    A := Evaluate(a, Expansion); printf "finished.\nExpanding b...";
    B := Evaluate(b, Expansion); printf "finished.\nEvaluation c...";
    C := Evaluate(c, Expansion); printf "finished.\nRemoving common factors...";
    
    aseq := Eltseq(A); bseq := Eltseq(B); cseq := Eltseq(C); printf "\n --finished eltseq\n";
    gcd := GreatestCommonDivisor(aseq cat bseq cat cseq); printf " --finished gcd\n";
    gcd := gcd*X^(TotalDegree(gcd));
    A := P!(A/gcd); B := P!(B/gcd); C := P!(C/gcd); printf " --finished quotients\n";
    
    printf "gcd: %o; ...finished\nEvaluating square...", gcd;

    sq := B^2 - 4*A*C; printf "finished.\nEvaluating root...";
    
    printf "Testing our square root: \n";
    
    return sq, (s2 + t2)^2*Expansion[0+1];
    
    root := OurSquareRoot(Q!sq, (s2 + t2)^2); printf "finished.\n";
    
    printf "Valuation of top: %o, valuation of bottom: %o\n", Valuation(-B + root), Valuation(A);
    
    result := Q!((-B - root)/(2*A));
    
    // We must clear any denominators:
    seq := Eltseq(result);
    denom := &*[Denominator(seq[i]) : i in [1..#seq]];
    // We want the lcm:
    if #seq ne 1 then 
        denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]]));
    end if;
    denom := denom*X^TotalDegree(denom);
    // Update all the other coordinates:
    for i in {1..16} do
        Expansion[i] *:= denom;
    end for;
    
    printf "result*denom = %o\n", result*denom;
    
    Expansion[5 + 1] := P!(result*denom);
    
    printf "Finished (Relative Precision: %o. Denominator: %o)", RelativePrecision(Expansion[5+1]), denom;
    
    UpdateCoords( ff, [15, 14, 13, 12, 11, 10, 5], [0, 9, 8, 7, 6, 4, 3, 2, 1], ~Expansion);
    
    return Expansion;
    
end intrinsic;

