// PreProcessingTM_GivenConductor.m

/*
INPUT:
    C:= [c_1, c_2, c_3, c_4], the coefficients of the Thue-Mahler form, F, where
        F:= c_1*x^3 + c_2*x^2*y + c_3*x*y^2 + c_4*y^3, having discriminant DF:= +/- 2^{alpha_0} * 3^{beta_0} * N_1
    N:= the conductor of the elliptic curves we seek
    RSolnFile:= text file storing the solutions of the form F which may lead to elliptic curves of conductor N
    
OUTPUT:
     p:= the set of rational primes [2,p_2, ..., p_s] such that F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3*p_4*...*p_s
        NB. the set of primes p is determined by the set of primes which divide N_0 of the conductor (ie. the factor of N which is coprime to 6), as well as alpha_1, beta_1
    
    [alpha_1, beta_1]:= [ [ [alpha_{1,1},A_1], .., [alpha_{1,r}, A_r] ], [ [beta_{1,1},B_1], .., [beta_{1,r}, B_r] ] ] where
        alpha_{1,i}:= the possible value of alpha_1 (ie. the exponent on 2) based on alpha_0 and alpha
        A_1:= 0 or 1, such that 
            if A_1 == 0 then
                alpha_1 is exactly equal to alpha_{1,i}
            if A_1 == 1 then
                alpha_1 is greater or equal to alpha_{1,i}
        beta_{1,i}:= the possible value of beta_1 (ie. the exponent on 3) based on beta_0 and beta
        B_1:= 0 or 1, such that 
            if B_1 == 0 then
                beta_1 is exactly equal to beta_{1,i}
            if B_1 == 1 then
                beta_1 is greater or equal to beta_{1,i}
     
     EnterTM:= a Boolean value determining whether or not to enter the UBC_TM solver
        NB. EnterTM:= false if there is a local obstruction, as in Theorem 1 of the reference
    
    if all A_i == 0:  
        Bounda:= the bound on alpha_1; ie. the largest value of alpha_{1,i}
    else:
        Bounda:= 2^15, a placeholder value, where
            Bounda:= the bound on alpha_1
                NB. here, not enough information is known to bound alpha_1
    
    if all B_i == 0:
        Boundb:= the bound on beta_1; ie. the largest value of beta_{1,i}
            Boundb:= the bound on beta_1
    else:
        Boundb:= 2^15, a placeholder value, where
            Bounda:= the bound on beta_1
                NB. here, not enough information is known to bound beta_1
                
COMMENTS:
    Given a Thue-Mahler form F, F:= c_1*x^3 + c_2*x^2*y + c_3*x*y^2 + c_4*y^3, having discriminant DF:= +/- 2^{alpha_0} * 3^{beta_0} * N_1, and conductor 
        N = 2^{alpha} * 3^{beta} * N_0, 
    one can compute solutions to 
        F(u,v) = 2^{alpha_1} * 3^{beta_1} * p_3^{k_3} * p_4^{k_4} *...* p_s^{k_s}, where p_3, ..., p_s | N_0
    leading to elliptic curves having good reduction outside p, as per Theorem 1 of the reference. 
    
    This algorithm determines a potential upper bound on alpha_1, beta_1 as per Theorem 1 of the reference based on alpha_1, beta_1, alpha, beta.
    
REFERENCE:
    M. Bennett, A. Gherga, A. Rechnitzer. Computing Elliptic Curves over Q. Submitted, University of British Columbia, 2017.

EXAMPLE:
    > C:= [1,0,30,2];
    > N:= 10^6-1;
    > RSolnFile:= "RelevantSolutions.txt";
    > p, alpha_1beta_1, EnterTM, Bounda, Boundb:= PreProcessingTM_GivenConductor(C,N,RSolnFile);    
    > p;
    [ 2, 3, 7, 11, 13, 37 ]
    > alpha_1:= alpha_1beta_1[1];
    > alpha_1;
    [
        [ 0, 0 ],
        [ 3, 0 ]
    ]
    > beta_1:=alpha_1beta_1[2];
    > beta_1;
    [
        [ 0, 0 ],
        [ 1, 0 ]
    ]
    > EnterTM;
    true
    > Bounda;
    3
    > Boundb;
    1

*/

function PreProcessingTM_GivenConductor(C,N,RSolnFile)
    p:= [f[1] : f in Factorization(N) | f[1] gt 3];     // computes the prime factors of the conductor, N; excludes primes 2,3 since those are determined by alpha_1, beta_1

    R<x,y>:= PolynomialRing(Integers(),2); 
    F:= C[1]*x^3 + C[2]*x^2*y + C[3]*x*y^2 + C[4]*y^3;  // generates the relevant Thue-Mahler polynomial, F
    DF:= -27*C[1]^2*C[4]^2 + C[2]^2*C[3]^2 + 18*C[1]*C[2]*C[3]*C[4] - 4*C[1]*C[3]^3 - 4*C[2]^3*C[4];    // computes the discriminant of F
    DF_fac:= [f[1] : f in Factorization(DF)];   // computes the prime factors of the discriminant of F
    
    alpha:= Valuation(N,2);     // determines alpha; the exponent on 2 in the conductor
    beta:= Valuation(N,3);      // determines beta; the exponent on 3 in the conductor
    
    alpha_0:= Valuation(DF,2);  // determines alpha_0; the exponent on 2 in disc(F)
    beta_0:= Valuation(DF,3);   // determines beta_0; the exponent on 3 in disk(F)
    
    if (alpha eq 0) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [0,0], [3,0] ];         
        end if;
    elif (alpha eq 1) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [4,1] ];        
        elif (alpha_0 eq 3) then
            alpha_1:= [ [3,1] ];        
        end if;
    elif (alpha eq 2) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [1,0] ];        
        elif (alpha_0 eq 4) then
            alpha_1:= [ [0,0], [1,0] ];         
        end if;
    elif (alpha eq 3) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [1,0], [2,0] ];
        elif (alpha_0 eq 3) then
            alpha_1:= [ [2,0] ];
        elif (alpha_0 eq 4) then
            alpha_1:= [ [0,0], [1,0] ];
        end if;
    elif (alpha eq 4) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [0,1] ]; 
        elif (alpha_0 eq 3) then
            alpha_1:= [ [2,1] ];
        elif (alpha_0 eq 4) then
            alpha_1:= [ [0,0],[1,0] ];
        end if;
    elif (alpha eq 5) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [2,0] ];
        elif (alpha_0 eq 3) then
            alpha_1:= [ [1,0] ];
        end if;
    elif (alpha eq 6) then
        if (alpha_0 eq 2) then
            alpha_1:= [ [0,1] ];
        elif (alpha_0 eq 3) then
            alpha_1:= [ [1,1] ];
        elif (alpha_0 eq 4) then
            alpha_1:= [ [0,0],[1,0] ];
        end if;
    elif (alpha eq 7) then
        if (alpha_0 eq 3) then
            alpha_1:= [ [0,0] ];
        elif (alpha_0 eq 4) then
            alpha_1:= [ [0,0] ];
        end if;
    elif (alpha eq 8) then
        if (alpha_0 eq 3) then
            alpha_1:= [ [1,0] ];
        end if;
    end if;
    
    if (beta eq 0) then
        if (beta_0 eq 0) then
            beta_1:= [ [0,0] ];
        end if;
    elif (beta eq 1) then
        if (beta_0 eq 0) then
            beta_1:= [ [1,1] ];
        elif (beta_0 eq 1) then
            beta_1:= [ [0,1] ];
        end if;
    elif (beta eq 2) then
        if (beta_0 eq 0) then
            beta_1:= [ [0,1] ];
        elif (beta_0 eq 1) then
            beta_1:= [ [0,1] ];
        elif (beta_0 eq 3) then
            beta_1:= [ [0,0] ];
        end if;
    elif (beta ge 3) then
        if (beta_0 eq beta) then
            beta_1:= [ [0,0],[1,0] ];
        else    // it may occur that the imprimitive form F_1 has all coefficients divisible by 3; in this case, since also beta_1 = {0,1} and 3|LHS, we must have that 3|RHS, hence beta_1 = 1
                    // in this case, may divide 3 from both sides, yielding the form F = F_1/3, whose discriminant has Valuation(DF,3) != beta_0 = beta
            beta_1:= [[0,0]];   // since beta_1 = 1 was divided out, necessarily, beta_1 = 0
        end if;
    end if;
    
    if &and[a[2] eq 0: a in alpha_1 ] then      // if all A_{1,i} == 0
        Bounda:= Max([a[1] : a in alpha_1]);
    else
        Bounda:= 2^15;  // placeholder, no upper bound exists on alpha
    end if;
    
    if &and[b[2] eq 0: b in beta_1 ] then       // if all B_{1,i} == 0
        Boundb:= Max([b[1] : b in beta_1]);
    else
        Boundb:= 2^15;  // placeholder, no upper bound exists on beta
    end if;
    
    for a0 in alpha_1 do        // runs through values of alpha_1
        if (a0[1] gt 0) and (2 notin p) then
            Append(~p,2);       // appends p = 2 if alpha_1 > 0, where F(u,v) = 2^{alpha_1} * 3^{beta_1} * {divisors of N_0}
        elif (a0[2] eq 1) and (2 notin p) then
            Append(~p,2);
        end if;
    end for;
    
    for b0 in beta_1 do         // runs through values of beta_1
        if (b0[1] gt 0) and (3 notin p) then        
            Append(~p,3);       // appends p = 3 if beta_1 > 0, where F(u,v) = 2^{alpha_1} * 3^{beta_1} * {divisors of N_0}
        elif (b0[2] eq 1) and (3 notin p) then
            Append(~p,3);
        end if;
    end for;
    
    Sort(~p);
    
    LocalObstruction:= [];      // stores primes q, where q is a local obstruction
    for q in [i : i in p | i ge 3] do   // runs through potential local obstructions, q; if v_q(N) = 1, for q >= 3, q|( D_F*F(u,v) )
        if (Valuation(DF,q) eq 0) and (Valuation(N,q) eq 1) then        // NB. if Valuation(D_F,q) != 0, then q|D_F, hence q|( D_F*F(u,v) )
            check:= [];         // stores [u,v] such that D_F*F(u,v) = 0 mod q 
            x:= 0;
           
            while x le (q-1) do
                y:= 0;
                while y le (q-1) do
                    if GCD(x,y) eq 1 then
                        F_q:= Integers()! ( Evaluate(F,[x,y]) );
                        if (( (DF*F_q) mod q ) eq 0) then
                            Append(~check, [x,y]);
                            y:= 2*(q-1);
                            x:= 2*(q-1); 
                        end if;
                    end if;
                    y:= y + 1;
                end while;
                x:= x + 1;
            end while;
            
            if (check eq [[0,0]]) or (check eq []) then
                Append(~LocalObstruction, q);
            end if;

        end if;
    end for;
    
    if (LocalObstruction ne []) then
        EnterTM:= false;        // local obstructions present; F(x,y) does not enter TM equation solver
        fprintf RSolnFile, "Relevant solutions corresponding to conductor N = %o, primes p:= %o, and form C:=%o \n ", N, p, C;
        fprintf RSolnFile, "=========================================================================== \n";
        
        fprintf RSolnFile, "Local Obstruction at %o. No solutions possible for conductor N = %o", LocalObstruction, N; 
    else
        EnterTM:= true;         // no local obstructions, solutions are possible; F(x,y) enters TM equation solver
    end if;
    
    return p, [alpha_1, beta_1], EnterTM, Bounda, Boundb;
end function;