12
\$\begingroup\$

The task is to compute OEIS A005434 as quickly as possible.

Consider a binary string S of length n. Indexing from 1, we can determine if S[1..i+1] matches S[n-i..n] exactly for all i in order from 0 to n-1. For example,

S = 01010

gives 

[Y, N, Y, N, Y].

This is because 0 matches 0, 01 does not match 10, 010 matches 010, 0101 does not match 1010 and finally 01010 matches itself.

Define f(n) to be the number of distinct arrays of Ys and Ns one gets when iterating over all 2^n different possible bit strings S of length n.

The observant will notice this question is a simpler variant of another recent question of mine. However, I expect that clever tricks can make this much faster and easier.

Task

For increasing n starting at 1, your code should output n, f(n).

Example answers

For n = 1..24, the correct answers are:

1, 2, 3, 4, 6, 8, 10, 13, 17, 21, 27, 30, 37, 47, 57, 62, 75, 87, 102, 116, 135, 155, 180, 194

Scoring

Your code should iterate up from n = 1 giving the answer for each n in turn. I will time the entire run, killing it after two minutes.

Your score is the highest n you get to in that time.

In the case of a tie, the first answer wins.

Where will my code be tested?

I will run your code under Virtualbox in a Lubuntu guest VM (on my Windows 7 host).

My laptop has 8GB of RAM and an Intel i7 [email protected] GHz (Broadwell) CPU with 2 cores and 4 threads. The instruction set includes SSE4.2, AVX, AVX2, FMA3 and TSX.

Leading entries per language

  • n = 599 in Rust bu Anders Kaseorg.
  • n= 30 in C by Grimy. Parallel version gets to 32 when run natively in cygwin.
Improve this question
\$\endgroup\$
4
  • \$\begingroup\$ math.uni-bielefeld.de/~sillke/SEQUENCES/autocorrelation-range.c (linked from the OEIS page) run with -O3 can calculate up to 100 in <.02 seconds on my machine \$\endgroup\$ Commented Jun 16, 2017 at 15:08
  • \$\begingroup\$ @rogaos Oh dear. I should delete the question but it already has an answer. \$\endgroup\$ Commented Jun 16, 2017 at 15:10
  • \$\begingroup\$ I think it's still a cool problem – but maybe up to 1000 instead? Or ask answers to golf a sufficiently fast program \$\endgroup\$ Commented Jun 16, 2017 at 15:11
  • 1
    \$\begingroup\$ @rogaos I just removed the hard limit completely. \$\endgroup\$ Commented Jun 16, 2017 at 15:19

3 Answers 3

Reset to default
4
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
2
  • \$\begingroup\$ Makes you wonder if there is a faster solution to the other related problem. Very impressive! \$\endgroup\$ Commented Jun 17, 2017 at 13:18
  • \$\begingroup\$ Interestingly this is entirely memory limited. It speeds up to ~500 and then suddenly slows down as it runs out of RAM. \$\endgroup\$ Commented Jun 17, 2017 at 15:57
2
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
4
  • \$\begingroup\$ @Lembik Is existence of a good answer outside of SE grounds for deletion? Shouldn’t you wait for somebody (possibly the commenter) to submit this algorithm as an answer, and accept that answer? \$\endgroup\$ Commented Jun 16, 2017 at 15:16
  • \$\begingroup\$ You make a very good point \$\endgroup\$ Commented Jun 16, 2017 at 15:16
  • \$\begingroup\$ Sadly I can't really test your parallel code in virtualbox as I have a total of two cores on my CPU. \$\endgroup\$ Commented Jun 16, 2017 at 15:54
  • \$\begingroup\$ I ran it in cygwin and it got to 32. \$\endgroup\$ Commented Jun 16, 2017 at 16:07
0
\$\begingroup\$

Python, n≈427

Rewrite @Anders Kaseorg's Rust answer in Python

TIO Link

import bisect from collections import defaultdict import time memo = defaultdict(int) s = [] def f(n, p, s): assert(p != 0) d = n // p assert(d >= 1) r = n - p * (d - 1 if d >= 2 else 1) k = bisect.bisect_left(s, n - r + 1) for i in s[:k]: if i % p != 0: return 0 if d >= 3: o = n - (p + r) n = p + r s = [i - o for i in s[k:]] elif n == p: return 1 elif k != 0: s = s[k:] query = (n, p, tuple(s)) if query in memo: return memo[query] t = [i - p for i in s] if d < 2: c = sum(f(r, q, t) for q in range(1, r + 1)) elif r == p: c = sum(f(r, q, t) for q in range(1, p + 1) if p % q != 0 or q == p) else: k = bisect.bisect_left(t, p) if k < len(t) and t[k] == p: pass else: t.insert(k, p) c = sum(f(r, q, t) for q in range(1, t[0] + 1) if p % q != 0 or q == p) memo[query] = c return c start_time = time.time() n = 1 while time.time() - start_time < 120: # run for at most 2 minutes print(n, sum(f(n, p, s) for p in range(1, n + 1))) n += 1
Improve this answer
\$\endgroup\$