GAP code for computing the number of monogenic transformation semigroups of a given degree

gistfile1.txt

# Sieve of Eratosthenes based on code from
# https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Pseudocode

MyPrimes := function(n)
  local A, limit, j, i;

  A := BlistList([1 .. n], [2 .. n]);
  limit := Int(Sqrt(Float(n)));
  for i in [1 .. limit] do 
    if A[i] then 
      j := i ^ 2; 
      while j <= n do 
        A[j] := false;
        j := j + i;
      od;
    fi;
  od;
  return ListBlist([1 .. n], A);
end;

PRIMES := MyPrimes(10 ^ 6);;

# Based on Maple code on
# http://oeis.org/A051613 

B := [];

a := function(n, i)
  local p, nr, j;

  if IsBound(B[n + 1]) then 
    if IsBound(B[n + 1][i + 1]) then 
      return B[n + 1][i + 1];
    fi;
  else 
    B[n + 1] := [];
  fi;
  
  if n = 0 then 
    B[n + 1][i + 1] := 1;
    return B[n + 1][i + 1];
  elif i < 1 then 
    B[n + 1][i + 1] := 0;
    return B[n + 1][i + 1];
  fi;
  p := PRIMES[i];
  
  nr := a(n, i -1);

  for j in [1 .. LogInt(n, p)] do 
    nr := nr + a(n - p ^ j, i - 1);
  od;
  
  B[n + 1][i + 1] := nr;
  return nr;
end;

# The b(n) is number of partitions of n into powers of distinct primes

b := n -> a(n, Length(MyPrimes(n)));

# Based on Vladeta Jovovic's formula for the number of cyclic subgroups of the
# symmetric group on n points on:
# http://oeis.org/A009490

LCMS := [];

NrCyclicPermGroups := function(n)
  if not IsBound(LCMS[n]) then 
    LCMS[n] := Sum(List([0 .. n], k -> b(k)));
  fi;
  return LCMS[n];
end;

NrMonogenicTransformationSemigroups := function(n)
  local nr, m;
  nr := 0;
  for m in [1 .. n] do
    nr := nr + NrCyclicPermGroups(m);
  od;
  return nr;
end;