def keccak_f_1600(state):
    """The inner permutation for the Keccak sponge function.

    The Keccak-f permutation is an iterated construction consisting of a
    sequence of almost identical rounds. It operates on a state array with
    each of the twenty-four rounds performing five steps, described below
    with detail.

    The loops above and below the core of the permutation are used to save and
    restore the state array to a stream of bytes, used outside the permutation.
    The original state array has three dimensions, whereas this characteristic
    can be cleverly optimized to a 5x5 matrix with 64-bit words. As such, this
    implementation makes use of this trick and stores an entire lane (a z-axis
    set of bits within the state) as a single word.

    The θ step diffuses the bits alongside the state array by calculating the
    parity of nearby columns relative to a lane.

    The ρ and π steps are merged; together, they move more bits around
    according to two alike recurrence relations.

    The χ step is similar to an S-box permutation; it makes the whole round
    non-linear with a few logic operations on bits inside a line.

    The ι step is a simple LFSR that breaks the symmetry of the rounds. It
    generates constants by doing computations according to the round number
    and its previous output, modulo polynomials over GF(2)[x].

    Args:
        state:  square matrix of order 5 that holds the input bytes.

    Returns:
        state:  bytes permuted by Keccak-f[1600].
    """

    def load64(b):
        """
        Saves each byte on its respective position within a 64-bit word.

        Args:
            b:  partial list of bytes from input.

        Returns:
            Sum of list with numbers shifted.
        """
        return sum((b[i] << (8 * i)) for i in range(8))

    def store64(a):
        """
        Transforms a 64-bit word into a list of bytes.

        Args:
            a:  64-bit word.

        Returns:
            List of bytes separated by position on the word.
        """
        return list((a >> (8 * i)) % 256 for i in range(8))

    def rotate(a, n):
        """
        Denotes the bitwise cyclic shift operation, moving bit at position
        `i` into position `i + n` (modulo the lane size).

        Args:
            a:  lane with a 64-bit word, or elements from the state array.
            n:  offset for rotation.

        Returns:
            The rotated lane.
        """
        return ((a >> (64 - (n % 64))) + (a << (n % 64))) % (1 << 64)

    A = [[0 for _ in range(5)] for _ in range(5)]
    for x in range(5):
        for y in range(5):
            i = 8 * (x + 5 * y)
            A[x][y] = load64(state[i : i + 8])

    R = 1
    for _ in range(24):
        C = [A[x][0] ^ A[x][1] ^ A[x][2] ^ A[x][3] ^ A[x][4] for x in range(5)]
        D = [C[(x - 1) % 5] ^ rotate(C[(x + 1) % 5], 1) for x in range(5)]
        A = [[A[x][y] ^ D[x] for y in range(5)] for x in range(5)]

        x, y, current = 1, 0, A[1][0]
        for t in range(24):
            x, y = y, (2 * x + 3 * y) % 5
            offset = ((t + 1) * (t + 2)) // 2
            current, A[x][y] = A[x][y], rotate(current, offset)

        for y in range(5):
            T = [A[x][y] for x in range(5)]
            for x in range(5):
                A[x][y] = T[x] ^ ((~T[(x + 1) % 5]) & T[(x + 2) % 5])

        for j in range(7):
            R = ((R << 1) ^ ((R >> 7) * 0x71)) % 256
            if R & 2:
                A[0][0] ^= 1 << ((1 << j) - 1)

    for x in range(5):
        for y in range(5):
            i = 8 * (x + 5 * y)
            state[i : i + 8] = store64(A[x][y])

    return state