Computes the number of conjugacy classes of 2023×2023 matrices with coefficients in the finite field 𝔽₂

cargo.toml

[dependencies]
num-bigint = "0.4.4"
num-traits = "0.2.17"

HNY.rs

use num_bigint::{BigUint, ToBigUint};
use num_traits::{One, Zero};

fn divisors(n: usize) -> Vec<usize> {
    let mut divs = Vec::new();
    for i in 1..n + 1 {
        if n % i == 0 {
            divs.push(i);
        }
    }
    divs
}

fn main() {
    const N: usize = 2023;
    let mut c: [BigUint; N + 1] = std::array::from_fn(|_| BigUint::zero());
    for n in 1..N + 1 {
        for d in divisors(n) {
            c[n] += d.to_biguint().unwrap() * (BigUint::one() << (n / d));
        }
    }
    let mut b: [BigUint; N + 1] = std::array::from_fn(|_| BigUint::zero());
    b[0] = BigUint::one();
    for n in 1..N + 1 {
        for k in 1..n {
            b[n] += c[k].clone() * b[n - k].clone();
        }
        b[n] += c[n].clone();
        b[n] /= n;
    }
    print!("{}", b[N]);
}

Runs in about 500ms on my laptop. The output is

3335046873222021964935649536263932166251363070724652188755090828719410225592325865029100255101556534520137476606182718362782153894752880601586793099979765585854303419865707625511654617471452154169628385033298652767392696241343665338383349316096952923944666268608598081258761073649034938897419485527764891333855911298593051397743384659830370155218107973666631111129675410551753616752682250460652644920438563247713847253064157913405521388784520267386239178951932922774580737523194652065934689417098801115075955318904014329469831391039802059119803778333722785092907529948537045754234185294262387375136964882773214