Heisenberg Group Distance Histogram

HeisenbergGroupDistanceHistogram.cs

Console.OutputEncoding = System.Text.Encoding.UTF8;
const int Order = 6;
int basisSize = Order * 2;

var reachabilities = GetReachabilities(Enumerable.Range(0, basisSize).Select(i => 1ul << (i)).ToArray());
Console.WriteLine("Reachabilities dump: ");
WriteArray(reachabilities);
Console.WriteLine(); 
Console.WriteLine("Histogram:"); 
PrintHistogram(reachabilities);
 
void WriteArray(int[] array)
{
    for (int i = 0; i < array.Length; i++)
    {
        Console.Write(array[i] == -1 ? "*" : array[i]);
    }
}
void PrintHistogram(int[] points)
{
    var distinct = points.Distinct().OrderBy(x => x).ToArray();
    for (int i = 0; i < distinct.Length; i++)
    {
        Console.WriteLine($"{distinct[i]}: {points.Count(x => x == distinct[i])}");
    }
} 
void WriteMatrix(int[,] matrix)
{
    for (int i = 0; i < matrix.GetLength(0); i++)
    {
        for (int j = 0; j < matrix.GetLength(1); j++)
        {
            Console.Write(matrix[i, j] == -1 ? "*" : matrix[i, j]);
        }
        Console.WriteLine();
    }
}
int[] GetReachabilities(ulong[] basis)
{
    ulong size = 1ul << (2 * Order + 1);
    int[] jumpCounts = new int[size];
    HashSet<ulong> currentVertices = new(basis);
    for (ulong i = 0; i < size; i++)
    {
        for (ulong j = 0; j < size; j++)
        {
            jumpCounts[i] = (currentVertices.Contains(i) ? 0 : -1);
        }
    }
    int currentJump = 0;
    bool changed;

    do
    {
        currentJump++;
        changed = false;
        HashSet<ulong> newVertices = new(currentVertices);
        foreach (var a in currentVertices)
            foreach (var b in currentVertices)
                newVertices.Add(Times(a, b));


        foreach (var vertex in newVertices)
        {
            if (jumpCounts[vertex] == -1 || jumpCounts[vertex] > currentJump)
            {
                jumpCounts[vertex] = currentJump;
                changed = true;
            }
        }
        currentVertices = newVertices;
    } while (changed);

    return jumpCounts;
}
ulong Times(ulong a, ulong b)
{
    var xored = a ^ b;

    var zMask = (1UL << (2 * Order));

    var xDotY = a & (b << Order); // in place of y lies the dot product
    for (int i = 0; i < Order; i++)
    {
        xDotY <<= 1;
        xored ^= (xDotY & zMask);
    }
    return xored;
}
ulong[] GetOrbit(ulong element)
{
    HashSet<ulong> result = new();
    ulong size = 1ul << (2 * Order + 1);
    for (ulong i = 0; i < size; i++)
    {
        result.Add(Times(element, i));
    }
    return result.ToArray();
}
void PrintElement(ulong element)
{
    Console.Write("1");
    for (int i = 0; i < Order; i++)
    {
        Console.Write((element & (1UL << i)) == 0 ? "0" : "1");
    }
    Console.WriteLine((element & (1UL << 2 * Order)) == 0 ? "0" : 1);
    for (int i = 0; i < Order; i++)
    {
        Console.WriteLine("1".PadLeft(i + 2).PadRight(Order - i - 1 + i + 2, '0') + ((element & (1UL << (i + Order))) == 0 ? "0" : "1"));
    }
    Console.WriteLine("1".PadLeft(2 * Order - 1));
}