///////////////////////////////////////////////////////////////////////////
// "Growth of torsion groups of elliptic curves upon base changes"
// Enrique González-Jiménez & Filip Najman
///////////////////////////////////////////////////////////////////////////

// 13/10/2016 - Magma 2.21

// Magma script related to Section 6


function C(p)
// Given a prime p, the function returns the set A with the property that it is the smallest set values a such that p is in R_Q(d) if and only if d is divisible by some a in A.

CM:={3,4,7,8,11,19,43,67,163};
  if p in {2,3,5,7} then return [1]; end if;
  if p in {11} then return [5]; end if;
  if p in {13} then return [3,4]; end if;
  if p in {17} then return [8]; end if;
  if p in {37} then return [12]; end if;
  if p in {19,43,67,163} then return [(p-1)/2]; end if;
  if p mod 3 eq 1 then return [2*(p-1)]; end if;
  if p mod 3 eq 2 and not &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [2*(p-1)]; end if;
  if p mod 9 in {2,5} and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [(p^2-1)/3]; end if;
  if p mod 9 eq 8 and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then 
    return [p^2-1]; 
  end if;
end function;

function RQ(d)
// Given a positive integer d, 
// the function returns the set of primes p for which 
// there exists a number field K with [K:Q]=d 
// and an elliptic curve E/Q , 
// such that p divides the order of the torsion subgroup of E(K)

  bound:=Integers()!Max((163-1)/2,2*d+1);
  R:={2,3,5,7};
  if d le 2 then return R; end if;
  for p in PrimesInInterval(2,bound) do Cp:=C(p); 
    for cp in Cp do
      if IsDivisibleBy(d,Integers()!cp) then R:=R join {p}; end if;
     end for; 
  end for;
  return R;

end function;



function RQstar0(d)
// Given a positive integer d, 
// the function returns the set of primes p 
// that appear in RQ(d) and not in RQ(dd), where dd runs through all proper divisors dd of d.
  bound:=Integers()!Max((163-1)/2,2*d+1);
  Rstar:=RQ(d);
  for D in [dd : dd in Divisors(d) | dd ne d] do
    Rstar:=Rstar diff RQ(D);   
  end for;
  return Rstar;
end function;

/* OUTPUT: Table below Corollary 6.4
> for d in [1..100] do R:=RQstar0(d); if #R ne 0 then print d, R; end if; end for;
1 { 2, 3, 5, 7 }
3 { 13 }
4 { 13 }
5 { 11 }
8 { 17 }
9 { 19 }
12 { 37 }
21 { 43 }
33 { 67 }
44 { 23 }
56 { 29 }
60 { 31 }
80 { 41 }
81 { 163 }
92 { 47 }
*/