Compute OEIS A005434
Rust, n ≈ 660
use std::collections::HashMap;
use std::iter::once;
use std::rc::Rc;
type Memo = HashMap<(u32, u32, Rc<Vec<u32>>), u64>;
fn f(memo: &mut Memo, mut n: u32, p: u32, mut s: Rc<Vec<u32>>) -> u64 {
debug_assert!(p != 0);
let d = n / p;
debug_assert!(d >= 1);
let r = n - p * if d >= 2 { d - 1 } else { 1 };
let k = s.binary_search(&(n - r + 1)).unwrap_or_else(|i| i);
for &i in &s[..k] {
if i % p != 0 {
return 0;
}
}
if d >= 3 {
let o = n - (p + r);
n = p + r;
s = Rc::new(s[k..].iter().map(|i| i - o).collect());
} else if n == p {
return 1;
} else if k != 0 {
s = Rc::new(s[k..].to_vec());
}
let query = (n, p, s);
if let Some(&c) = memo.get(&query) {
return c;
}
let (n, p, s) = query;
let t = Rc::new(s.iter().map(|i| i - p).collect::<Vec<_>>());
let c = if d < 2 {
(1..r + 1).map(|q| f(memo, r, q, t.clone())).sum()
} else if r == p {
(1..p + 1)
.filter(|&q| p % q != 0 || q == p)
.map(|q| f(memo, r, q, t.clone()))
.sum()
} else {
let t = match t.binary_search(&p) {
Ok(_) => t,
Err(k) => {
Rc::new(t[..k]
.iter()
.cloned()
.chain(once(p))
.chain(t[k..].iter().cloned())
.collect::<Vec<_>>())
}
};
(1..t.first().unwrap() + 1)
.filter(|&q| p % q != 0 || q == p)
.map(|q| f(memo, r, q, t.clone()))
.sum()
};
memo.insert((n, p, s), c);
c
}
fn main() {
let mut memo = HashMap::new();
let s = Rc::new(Vec::new());
for n in 1.. {
println!("{} {}",
n,
(1..n + 1)
.map(|p| f(&mut memo, n, p, s.clone()))
.sum::<u64>());
}
}
Try it online!
How it works
This is a memoized implementation of the recursive predicate Ξ given in Leo Guibas, “Periods in strings” (1981). The function f(memo, n, p, s)
finds the number of correlations of length n
with smallest period p
and also each of the periods in the set s
.
Just a simple brute-force search, to get the challenge started:
#include <stdio.h>
#include <stdint.h>
#include <string.h>
typedef uint16_t u16;
typedef uint64_t u64;
static u64 map[1<<16];
int main(void)
{
for (u64 n = 1;; ++n) {
u64 result = 1;
u64 mask = (1ul << n) - 1;
memset(map, 0, sizeof(map));
#pragma omp parallel
#pragma omp for
for (u64 x = 1ul << (n - 1); x < 1ul << n; ++x) {
u64 r = 0;
for (u64 i = 1; i < n; ++i)
r |= (u64) (x >> i == (x & (mask >> i))) << i;
if (!r)
continue;
u16 h = (u16) (r ^ r >> 13 ^ r >> 27);
while (map[h] && map[h] != r)
++h;
if (!map[h]) {
#pragma omp critical
if (!map[h]) {
map[h] = r;
++result;
}
}
}
printf("%ld\n", result);
}
}
Compile with clang -fopenmp -Weverything -O3 -march=native
. On my machine it reaches n=34 in 2 minutes.
EDIT: sprinkled some OMP directives for easy parallelism.