/* I have this file saved as .magmarc in my Unix home
It provides Q = rationals, R[t] = polynomial ring, K[t] = rational
function field, S[[s]] = power series ring, together with the
functions Qorb, Denom, IsgorensteinSymmetric. */

SetIgnorePrompt(true);
Q:=Rationals();
 // K<om> := CyclotomicField(3);
R<t>:=PolynomialRing(Q);
K:=FieldOfFractions(R);
S<s>:=PowerSeriesRing(Q,50);

function Denom(L)
 return &*[1-t^i : i in L];
end function;

function Porb(r,LL,k)
 L := [Integers() | i : i in LL]; // this allows empty list
 if (k + &+L) mod r ne 0
    then error "Error: Canonical weight not compatible";
 end if;
n := #L;
S := [GCD(a,r) : a in L];
D := (1-t^r) * &*[1-t^s : s in S]; // Denom
A := &*[(1-t^(L[i])) div (1-t^(S[i])) : i in [1..n]];
F := (1-t^r) div GCD(1-t^r, &*[1-t^s: s in S]);
dF := Degree(F);
shift0 := Ceiling((k + 1 + &+[s : s in S] + r - dF)/2);
de := Maximum(0,Ceiling(-shift0/r));
shift := shift0+de*r;
 G, al, be := XGCD(t^shift*A, F);
 return t^shift*al/(D*t^(de*r));
end function;

/* I hope that this function can be applied to a polynomial f = N/D
or a rational function. You win if f has Numerator and Denominator
palindromic or antipalindromic polynomials. */
 
function IsgorensteinSymmetric(f)
 K<t> := FunctionField(Q);
 fie := hom<K->K|1/t>;
 g:=f/fie(K!f);
 if #Terms(Numerator(g)) eq 1 then
  return true, Degree(Numerator(g))-Degree(Denominator(g));
 else return false;
 end if;
end function;

function IsGev(f)
 K<t> := FunctionField(Q);
 fie := hom<K->K|1/t>;
 g:=f/fie(K!f);
 if #Terms(Numerator(g)) eq 1 then
  if IsEven(Degree(Numerator(g))-Degree(Denominator(g))) then
   return true, Degree(Numerator(g))-Degree(Denominator(g));
   else return false; end if;
  else return false;
 end if;
end function;

function IsGod(f)
 K<t> := FunctionField(Q);
 fie := hom<K->K|1/t>;
 g:=f/fie(K!f);
 if #Terms(Numerator(g)) eq 1 then
  if IsEven(Degree(Numerator(g))-Degree(Denominator(g))) then
   return false;
   else return true, Degree(Numerator(g))-Degree(Denominator(g)); end if;
  else return false;
 end if;
end function;

function Ph(n)
 return(R!CyclotomicPolynomial(n));
end function;