14
\$\begingroup\$

Consider a binary string S of length n. Indexing from 1, we can compute the Hamming distances between S[1..i+1] and S[n-i..n] for all i in order from 0 to n-1. The Hamming distance between two strings of equal length is the number of positions at which the corresponding symbols are different. For example,

S = 01010

gives

[0, 2, 0, 4, 0].

This is because 0 matches 0, 01 has Hamming distance two to 10, 010 matches 010, 0101 has Hamming distance four to 1010 and finally 01010 matches itself.

We are only interested in outputs where the Hamming distance is at most 1, however. So in this task we will report a Y if the Hamming distance is at most one and an N otherwise. So in our example above we would get

[Y, N, Y, N, Y]

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.

Task

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

Example answers

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

1, 1, 2, 4, 6, 8, 14, 18, 27, 36, 52, 65, 93, 113, 150, 188, 241, 279, 377, 427, 540, 632, 768, 870

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 on my (slightly old) Windows 7 laptop under cygwin. As a result, please give any assistance you can to help make this easy.

My laptop has 8GB of RAM and an Intel i7 5600U@2.6 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 = 40 in Rust using CryptoMiniSat, by Anders Kaseorg. (In Lubuntu guest VM under Vbox.)
  • n = 35 in C++ using the BuDDy library, by Christian Seviers. (In Lubuntu guest VM under Vbox.)
  • n = 34 in Clingo by Anders Kaseorg. (In Lubuntu guest VM under Vbox.)
  • n = 31 in Rust by Anders Kaseorg.
  • n = 29 in Clojure by NikoNyrh.
  • n = 29 in C by bartavelle.
  • n = 27 in Haskell by bartavelle
  • n = 24 in Pari/gp by alephalpha.
  • n = 22 in Python 2 + pypy by me.
  • n = 21 in Mathematica by alephalpha. (Self reported)

Future bounties

I will now give a bounty of 200 points for any answer that gets up to n = 80 on my machine in two minutes.

\$\endgroup\$
  • \$\begingroup\$ Do you know of some trick that will allow someone to find a faster algorithm than a naive brute force? If not this challenge is "please implement this in x86" (or maybe if we know your GPU...). \$\endgroup\$ – Jonathan Allan Jun 4 '17 at 7:15
  • \$\begingroup\$ @JonathanAllan It is certainly possible to speed up a very naive approach. Exactly how fast you can get I am not sure. Interestingly, if we changed the question so that you get a Y if the Hamming distance is at most 0 and an N otherwise, then there is a known closed form formula. \$\endgroup\$ – user9206 Jun 4 '17 at 8:16
  • \$\begingroup\$ @Lembik Do we measure CPU time or real time? \$\endgroup\$ – flawr Jun 11 '17 at 19:32
  • \$\begingroup\$ @flawr I am measuring real time but running it a few times and taking the minimum to eliminate oddities. \$\endgroup\$ – user9206 Jun 11 '17 at 19:57

12 Answers 12

9
+100
\$\begingroup\$

Rust + CryptoMiniSat, n ≈ 41

src/main.rs

extern crate cryptominisat;
extern crate itertools;

use std::iter::once;
use cryptominisat::{Lbool, Lit, Solver};
use itertools::Itertools;

fn make_solver(n: usize) -> (Solver, Vec<Lit>) {
    let mut solver = Solver::new();
    let s: Vec<Lit> = (1..n).map(|_| solver.new_var()).collect();
    let d: Vec<Vec<Lit>> = (1..n - 1)
        .map(|k| {
                 (0..n - k)
                     .map(|i| (if i == 0 { s[k - 1] } else { solver.new_var() }))
                     .collect()
             })
        .collect();
    let a: Vec<Lit> = (1..n - 1).map(|_| solver.new_var()).collect();
    for k in 1..n - 1 {
        for i in 1..n - k {
            solver.add_xor_literal_clause(&[s[i - 1], s[k + i - 1], d[k - 1][i]], true);
        }
        for t in (0..n - k).combinations(2) {
            solver.add_clause(&t.iter()
                                   .map(|&i| d[k - 1][i])
                                   .chain(once(!a[k - 1]))
                                   .collect::<Vec<_>>()
                                   [..]);
        }
        for t in (0..n - k).combinations(n - k - 1) {
            solver.add_clause(&t.iter()
                                   .map(|&i| !d[k - 1][i])
                                   .chain(once(a[k - 1]))
                                   .collect::<Vec<_>>()
                                   [..]);
        }
    }
    (solver, a)
}

fn search(n: usize,
          solver: &mut Solver,
          a: &Vec<Lit>,
          assumptions: &mut Vec<Lit>,
          k: usize)
          -> usize {
    match solver.solve_with_assumptions(assumptions) {
        Lbool::True => search_sat(n, solver, a, assumptions, k),
        Lbool::False => 0,
        Lbool::Undef => panic!(),
    }
}

fn search_sat(n: usize,
              solver: &mut Solver,
              a: &Vec<Lit>,
              assumptions: &mut Vec<Lit>,
              k: usize)
              -> usize {
    if k >= n - 1 {
        1
    } else {
        let s = solver.is_true(a[k - 1]);
        assumptions.push(if s { a[k - 1] } else { !a[k - 1] });
        let c = search_sat(n, solver, a, assumptions, k + 1);
        assumptions.pop();
        assumptions.push(if s { !a[k - 1] } else { a[k - 1] });
        let c1 = search(n, solver, a, assumptions, k + 1);
        assumptions.pop();
        c + c1
    }
}

fn f(n: usize) -> usize {
    let (mut solver, proj) = make_solver(n);
    search(n, &mut solver, &proj, &mut vec![], 1)
}

fn main() {
    for n in 1.. {
        println!("{}: {}", n, f(n));
    }
}

Cargo.toml

[package]
name = "correlations-cms"
version = "0.1.0"
authors = ["Anders Kaseorg <andersk@mit.edu>"]

[dependencies]
cryptominisat = "5.0.1"
itertools = "0.6.0"

How it works

This does a recursive search through the tree of all partial assignments to prefixes of the Y/N array, using a SAT solver to check at each step whether the current partial assignment is consistent and prune the search if not. CryptoMiniSat is the right SAT solver for this job due to its special optimizations for XOR clauses.

The three families of constraints are

SiSk + iDki, for 1 ≤ kn − 2, 0 ≤ i ≤ nk;
Dki1Dki2 ∨ ¬Ak, for 1 ≤ kn − 2, 0 ≤ i1 < i2nk;
¬Dki1 ∨ ⋯ ∨ ¬Dkink − 1Ak, for 1 ≤ kn − 2, 0 ≤ i1 < ⋯ < ink − 1nk;

except that, as an optimization, S0 is forced to false, so that Dk0 is simply equal to Sk.

\$\endgroup\$
  • 2
    \$\begingroup\$ Woohoooooo ! :) \$\endgroup\$ – user9206 Jun 12 '17 at 11:38
  • \$\begingroup\$ I am still trying to compile this in Windows (using cygwin + gcc). I cloned cryptominisat and compiled it. But I still don't know how to compile the rust code. When I do cargo build I get --- stderr CMake Error: Could not create named generator Visual Studio 14 2015 Win64 \$\endgroup\$ – user9206 Jun 12 '17 at 14:17
  • 2
    \$\begingroup\$ @rahnema1 Thanks, but it sounds like the issue is with the CMake build system of the embedded C++ library in the cryptominisat crate, not with Rust itself. \$\endgroup\$ – Anders Kaseorg Jun 12 '17 at 21:26
  • 1
    \$\begingroup\$ @Lembik I'm getting a 404 from that paste. \$\endgroup\$ – Mego Jun 14 '17 at 9:17
  • 1
    \$\begingroup\$ @ChristianSievers Good question. That works but it seems to be a bit slower (2× or so). I’m not sure why it shouldn’t be just as good, so maybe CryptoMiniSat just hasn’t been well optimized for that kind of incremental workload. \$\endgroup\$ – Anders Kaseorg Jun 16 '17 at 10:08
9
\$\begingroup\$

Rust, n ≈ 30 or 31 or 32

On my laptop (two cores, i5-6200U), this gets through n = 1, …, 31 in 53 seconds, using about 2.5 GiB of memory, or through n = 1, …, 32 in 105 seconds, using about 5 GiB of memory. Compile with cargo build --release and run target/release/correlations.

src/main.rs

extern crate rayon;

type S = u32;
const S_BITS: u32 = 32;

fn cat(mut a: Vec<S>, mut b: Vec<S>) -> Vec<S> {
    if a.capacity() >= b.capacity() {
        a.append(&mut b);
        a
    } else {
        b.append(&mut a);
        b
    }
}

fn search(n: u32, i: u32, ss: Vec<S>) -> u32 {
    if ss.is_empty() {
        0
    } else if 2 * i + 1 > n {
        search_end(n, i, ss)
    } else if 2 * i + 1 == n {
        search2(n, i, ss.into_iter().flat_map(|s| vec![s, s | 1 << i]))
    } else {
        search2(n,
                i,
                ss.into_iter()
                    .flat_map(|s| {
                                  vec![s,
                                       s | 1 << i,
                                       s | 1 << n - i - 1,
                                       s | 1 << i | 1 << n - i - 1]
                              }))
    }
}

fn search2<SS: Iterator<Item = S>>(n: u32, i: u32, ss: SS) -> u32 {
    let (shift, mask) = (n - i - 1, !(!(0 as S) << i + 1));
    let close = |s: S| {
        let x = (s ^ s >> shift) & mask;
        x & x.wrapping_sub(1) == 0
    };
    let (ssy, ssn) = ss.partition(|&s| close(s));
    let (cy, cn) = rayon::join(|| search(n, i + 1, ssy), || search(n, i + 1, ssn));
    cy + cn
}

fn search_end(n: u32, i: u32, ss: Vec<S>) -> u32 {
    if i >= n - 1 { 1 } else { search_end2(n, i, ss) }
}

fn search_end2(n: u32, i: u32, mut ss: Vec<S>) -> u32 {
    let (shift, mask) = (n - i - 1, !(!(0 as S) << i + 1));
    let close = |s: S| {
        let x = (s ^ s >> shift) & mask;
        x & x.wrapping_sub(1) == 0
    };
    match ss.iter().position(|&s| close(s)) {
        Some(0) => {
            match ss.iter().position(|&s| !close(s)) {
                Some(p) => {
                    let (ssy, ssn) = ss.drain(p..).partition(|&s| close(s));
                    let (cy, cn) = rayon::join(|| search_end(n, i + 1, cat(ss, ssy)),
                                               || search_end(n, i + 1, ssn));
                    cy + cn
                }
                None => search_end(n, i + 1, ss),
            }
        }
        Some(p) => {
            let (ssy, ssn) = ss.drain(p..).partition(|&s| close(s));
            let (cy, cn) = rayon::join(|| search_end(n, i + 1, ssy),
                                       || search_end(n, i + 1, cat(ss, ssn)));
            cy + cn
        }
        None => search_end(n, i + 1, ss),
    }
}

fn main() {
    for n in 1..S_BITS + 1 {
        println!("{}: {}", n, search(n, 1, vec![0, 1]));
    }
}

Cargo.toml

[package]
name = "correlations"
version = "0.1.0"
authors = ["Anders Kaseorg <andersk@mit.edu>"]

[dependencies]
rayon = "0.7.0"

Try it online!

I also have a slightly slower variant using very much less memory.

\$\endgroup\$
  • \$\begingroup\$ What optimisations have you used? \$\endgroup\$ – user9206 Jun 5 '17 at 5:46
  • 1
    \$\begingroup\$ @Lembik The biggest optimization, besides doing everything with bitwise arithmetic in a compiled language, is to use only as much nondeterminism as needed to nail down a prefix of the Y/N array. I do a recursive search on possible prefixes of the Y/N array, taking along a vector of possible strings achieving that prefix, but only the strings whose unexamined middle is filled with zeros. That said, this is still an exponential search, and these optimizations only speed it up by polynomial factors. \$\endgroup\$ – Anders Kaseorg Jun 5 '17 at 6:15
  • \$\begingroup\$ It's a nice answer. Thank you. I am hoping someone will dig into the combinatorics to get a significant speed up. \$\endgroup\$ – user9206 Jun 5 '17 at 13:34
  • \$\begingroup\$ @Lembik I’ve fixed a memory wasting bug, done more micro-optimization, and added parallelism. Please retest when you get a chance—I’m hoping to increase my score by 1 or 2. Do you have combinatorial ideas in mind for larger speedups? I haven’t come up with anything. \$\endgroup\$ – Anders Kaseorg Jun 6 '17 at 3:40
  • 1
    \$\begingroup\$ @Lembik There is no formula given at the OEIS entry. (The Mathematica code there also seem to use brute-force.) If you know of one, you might want to tell them about it. \$\endgroup\$ – Christian Sievers Jun 10 '17 at 20:33
6
\$\begingroup\$

C++ using the BuDDy library

A different approach: have a binary formula (as binary decision diagram) that takes the bits of S as input and is true iff that gives some fixed values of Y or N at certain selected positions. If that formula is not constant false, select a free position and recurse, trying both Y and N. If there is no free position, we have found a possible output value. If the formula is constant false, backtrack.

This works relatively reasonable because there are so few possible values so that we can often backtrack early. I tried a similar idea with a SAT solver, but that was less successful.

#include<vector>
#include<iostream>
#include<bdd.h>

// does vars[0..i-1] differ from vars[n-i..n-1] in at least two positions?
bdd cond(int i, int n, const std::vector<bdd>& vars){
  bdd x1 { bddfalse };
  bdd xs { bddfalse };
  for(int k=0; k<i; ++k){
    bdd d { vars[k] ^ vars[n-i+k] };
    xs |= d & x1;
    x1 |= d;
  }
  return xs;
}

void expand(int i, int n, int &c, const std::vector<bdd>& conds, bdd x){
  if (x==bddfalse)
    return;
  if (i==n-2){
    ++c;
    return;
  }

  expand(i+1,n,c,conds, x & conds[2*i]);
  x &= conds[2*i+1];
  expand(i+1,n,c,conds, x);
}

int count(int n){
  if (n==1)   // handle trivial case
    return 1;
  bdd_setvarnum(n-1);
  std::vector<bdd> vars {};
  vars.push_back(bddtrue); // assume first bit is 1
  for(int i=0; i<n-1; ++i)
    if (i%2==0)            // vars in mixed order
      vars.push_back(bdd_ithvar(i/2));
    else
      vars.push_back(bdd_ithvar(n-2-i/2));
  std::vector<bdd> conds {};
  for(int i=n-1; i>1; --i){ // handle long blocks first
    bdd cnd { cond(i,n,vars) };
    conds.push_back( cnd );
    conds.push_back( !cnd );
  }
  int c=0;
  expand(0,n,c,conds,bddtrue);
  return c;
}

int main(void){
  bdd_init(20000000,1000000);
  bdd_gbc_hook(nullptr); // comment out to see GC messages
  for(int n=1; ; ++n){
    std::cout << n << " " << count(n) << "\n" ;
  }
}

To compile with debian 8 (jessie), install libbdd-dev and do g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. It might be useful to increase the first argument to bdd_init even more.

\$\endgroup\$
  • \$\begingroup\$ This looks interesting. What do you get to with this? \$\endgroup\$ – user9206 Jun 11 '17 at 18:24
  • \$\begingroup\$ @Lembik I get 31 in 100s on very old hardware that won't let me answer faster \$\endgroup\$ – Christian Sievers Jun 11 '17 at 18:50
  • \$\begingroup\$ Any help you can give on how to compile this on Windows (e.g. using cygwin) gratefully received. \$\endgroup\$ – user9206 Jun 11 '17 at 18:53
  • \$\begingroup\$ @Lembik I don't know about Windws but github.com/fd00/yacp/tree/master/buddy seems helpful w.r.t. cygwin \$\endgroup\$ – Christian Sievers Jun 11 '17 at 19:25
  • 1
    \$\begingroup\$ Wow, okay, you’ve got me convinced that I need to add this library to my toolkit. Well done! \$\endgroup\$ – Anders Kaseorg Jun 12 '17 at 8:52
4
\$\begingroup\$

Clingo, n ≈ 30 or 31 34

I was a little surprised to see five lines of Clingo code overtake my brute-force Rust solution and come really close to Christian’s BuDDy solution—it looks like it would beat that too with a higher time limit.

corr.lp

{s(2..n)}.
d(K,I) :- K=1..n-2, I=1..n-K, s(I), not s(K+I).
d(K,I) :- K=1..n-2, I=1..n-K, not s(I), s(K+I).
a(K) :- K=1..n-2, {d(K,1..n-K)} 1.
#show a/1.

corr.sh

#!/bin/bash
for ((n=1;;n++)); do
    echo "$n $(clingo corr.lp --project -q -n 0 -c n=$n | sed -n 's/Models *: //p')"
done

plot

\$\endgroup\$
  • \$\begingroup\$ This is great! It looks from your graph that the BuDDy solution suddenly gets worse. Any idea why? \$\endgroup\$ – user9206 Jun 12 '17 at 8:53
  • \$\begingroup\$ @Lembik I haven’t investigated BuDDy enough to be sure, but maybe it runs out of cache at that point? \$\endgroup\$ – Anders Kaseorg Jun 12 '17 at 8:59
  • \$\begingroup\$ Wow! I think a higher first value to bdd_init might help, or allowing to increase the node table more by calling bdd_setmaxincrease with a value much above the default of 50000. - Are you using the changed version of my program? \$\endgroup\$ – Christian Sievers Jun 12 '17 at 9:51
  • \$\begingroup\$ @ChristianSievers I’m using your current version. \$\endgroup\$ – Anders Kaseorg Jun 12 '17 at 10:00
  • 2
    \$\begingroup\$ I do love your graph. \$\endgroup\$ – user9206 Jun 12 '17 at 17:22
2
\$\begingroup\$

Pari/GP, 23

By default, Pari/GP limits its stack size to 8 MB. The first line of the code, default(parisize, "4g"), sets this limit to 4 GB. If it still gives a stackoverflow, you can set it to 8 GB.

default(parisize, "4g")
f(n) = #vecsort([[2 > hammingweight(bitxor(s >> (n-i) , s % 2^i)) | i <- [2..n-1]] | s <- [0..2^(n-1)]], , 8)
for(n = 1, 100, print(n " -> " f(n)))
\$\endgroup\$
  • \$\begingroup\$ Reaches 22 and then gives a stackoverflow. \$\endgroup\$ – user9206 Jun 4 '17 at 17:29
  • \$\begingroup\$ Gets to 24 now. \$\endgroup\$ – user9206 Jun 5 '17 at 8:39
2
\$\begingroup\$

Clojure, 29 in 75 38 seconds, 30 in 80 and 31 in 165

Runtimes from Intel i7 6700K, memory usage is less than 200 MB.

project.clj (uses com.climate/claypoole for multithreading):

(defproject tests "0.0.1-SNAPSHOT"
  :description "FIXME: write description"
  :dependencies [[org.clojure/clojure "1.8.0"]
                 [com.climate/claypoole "1.1.4"]]
  :javac-options ["-target" "1.6" "-source" "1.6" "-Xlint:-options"]
  :aot [tests.core]
  :main tests.core)

Source code:

(ns tests.core
  (:require [com.climate.claypoole :as cp]
            [clojure.set])
  (:gen-class))

(defn f [N]
  (let [n-threads   (.. Runtime getRuntime availableProcessors)
        mask-offset (- 31 N)
        half-N      (quot N 2)
        mid-idx     (bit-shift-left 1 half-N)
        end-idx     (bit-shift-left 1 (dec N))
        lower-half  (bit-shift-right 0x7FFFFFFF mask-offset)
        step        (bit-shift-left 1 12)
        bitcount
          (fn [n]
            (loop [i 0 result 0]
              (if (= i N)
                result
                (recur
                  (inc i)
                  (-> n
                      (bit-xor (bit-shift-right n i))
                      (bit-and (bit-shift-right 0x7FFFFFFF (+ mask-offset i)))
                      Integer/bitCount
                      (< 2)
                      (if (+ result (bit-shift-left 1 i))
                          result))))))]
    (->>
      (cp/upfor n-threads [start (range 0 end-idx step)]
        (->> (for [i      (range start (min (+ start step) end-idx))
                   :when  (<= (Integer/bitCount (bit-shift-right i mid-idx))
                              (Integer/bitCount (bit-and         i lower-half)))]
               (bitcount i))
             (into #{})))
      (reduce clojure.set/union)
      count)))

(defn -main [n]
  (let [n-iters 5]
    (println "Calculating f(n) from 1 to" n "(inclusive)" n-iters "times")
    (doseq [i (range n-iters)]
      (->> n read-string inc (range 1) (map f) doall println time)))
  (shutdown-agents)
  (System/exit 0))

A brute-force solution, each thread goes over a subset of the range (2^12 items) and builds a set of integer values which indicate detected patterns. These are then "unioned" together and thus the distinct count is calculated. I hope the code isn't too tricky to follow eventhough it uses threading macros a lot. My main runs the test a few times to get JVM warmed up.

Update: Iterating over only half of the integers, gets the same result due to symmetry. Also skipping numbers with higher bit count on lower half of the number as they produce duplicates as well.

Pre-built uberjar (v1) (3.7 MB):

$ wget https://s3-eu-west-1.amazonaws.com/nikonyrh-public/misc/so-124424-v2.jar
$ java -jar so-124424-v2.jar 29
Calculating f(n) from 1 to 29 (inclusive) 5 times
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 41341.863703 msecs"
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 37752.118265 msecs"
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 38568.406528 msecs"
[ctrl+c]

Results on different hardwares, expected runtime is O(n * 2^n)?

i7-6700K  desktop: 1 to 29 in  38 seconds
i7-6820HQ laptop:  1 to 29 in  43 seconds
i5-3570K  desktop: 1 to 29 in 114 seconds

You can easily make this single-threaded and avoid that 3rd party dependency by using the standard for:

(for [start (range 0 end-idx step)]
  ... )

Well the built-in pmap also exists but claypoole has more features and tunability.

\$\endgroup\$
  • \$\begingroup\$ Yeah, it makes it trivial to distribute. Would you have time to re-evaluate my solution, I'm quite sure you'd get it up-to 30 now. I do not have further optimizations in sight. \$\endgroup\$ – NikoNyrh Jun 6 '17 at 15:32
  • \$\begingroup\$ Sadly it's a no for 30. Elapsed time: 217150.87386 msecs \$\endgroup\$ – user9206 Jun 6 '17 at 18:14
  • \$\begingroup\$ Ahaa, thanks for giving it a try :D It might have been better to fit a curve on this and interpolate that at which decimal value 120 seconds is spent but even as it is this is a nice challgenge. \$\endgroup\$ – NikoNyrh Jun 6 '17 at 18:29
1
\$\begingroup\$

Mathematica, n=19

press alt+. to abort and the result will be printed

k = 0;
For[n = 1, n < 1000, n++,
Z = Table[HammingDistance[#[[;; i]], #[[-i ;;]]], {i, Length@#}] & /@
Tuples[{0, 1}, n];
Table[If[Z[[i, j]] < 2, Z[[i, j]] = 0, Z[[i, j]] = 1], {i, 
Length@Z}, {j, n}];
k = Length@Union@Z]
Print["f(", n, ")=", k]
\$\endgroup\$
  • \$\begingroup\$ I can't run this so could you explain how it avoids taking exponential time? 2^241 is a very big number! \$\endgroup\$ – user9206 Jun 4 '17 at 3:57
  • \$\begingroup\$ Can you show the output of the code? \$\endgroup\$ – user9206 Jun 4 '17 at 4:17
  • 1
    \$\begingroup\$ I meant f(n)... fixed \$\endgroup\$ – J42161217 Jun 4 '17 at 7:41
1
\$\begingroup\$

Mathematica, 21

f[n_] := Length@
     DeleteDuplicates@
      Transpose@
       Table[2 > Tr@IntegerDigits[#, 2] & /@ 
         BitXor[BitShiftRight[#, n - i], Mod[#, 2^i]], {i, 1, 
         n - 1}] &@Range[0, 2^(n - 1)];
Do[Print[n -> f@n], {n, Infinity}]

For comparison, Jenny_mathy's answer gives n = 19 on my computer.

The slowest part is Tr@IntegerDigits[#, 2] &. It is a shame that Mathematica doesn't have a built-in for Hamming weight.


If you want to test my code, you can download a free trial of Mathematica.

\$\endgroup\$
1
\$\begingroup\$

A C version, using builtin popcount

Works better with clang -O3, but also works if you only have gcc.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

unsigned long pairs(unsigned int n, unsigned long s) { 
  unsigned long result = 0;

  for(int d=1;d<=n;d++) { 
    unsigned long mx = 1 << d;
    unsigned long mask = mx - 1;

    unsigned long diff = (s >> (n - d)) ^ (s & mask);
    if (__builtin_popcountl(diff) <= 1)
      result |= mx;
  } 
  return result;

}

unsigned long f(unsigned long  n) { 
  unsigned long max = 1 << (n - 1);
#define BLEN (max / 2)
  unsigned char * buf = malloc(BLEN);
  memset(buf, 0, BLEN);
  unsigned long long * bufll = (void *) buf;

  for(unsigned long i=0;i<=max;i++) { 
    unsigned int r = pairs(n, i);
    buf[r / 8] |= 1 << (r % 8);
  } 

  unsigned long result = 0;

  for(unsigned long i=0;i<= max / 2 / sizeof(unsigned long long); i++) { 
    result += __builtin_popcountll(bufll[i]);
  } 

  free(buf);

  return result;
}

int main(int argc, char ** argv) { 
  unsigned int n = 1;

  while(1) { 
    printf("%d %ld\n", n, f(n));
    n++;
  } 
  return 0;
}
\$\endgroup\$
  • \$\begingroup\$ It gets to 24 very quickly and then ends. You need to increase the limit. \$\endgroup\$ – user9206 Jun 12 '17 at 13:16
  • \$\begingroup\$ Oh god, I forgot to remove the benchmark code! I'll remove the two offending lines :/ \$\endgroup\$ – bartavelle Jun 12 '17 at 13:18
  • \$\begingroup\$ @Lembik should be fixed now \$\endgroup\$ – bartavelle Jun 12 '17 at 13:20
1
\$\begingroup\$

Haskell, (unofficial n=20)

This is just the naive approach - so far without any optimizations. I wondered how well it would fare against other languages.

How to use it (assuming you have haskell platform installed):

  • Paste the code in one file approx_corr.hs (or any other name, modify following steps accordingly)
  • Navigate to the file and execute ghc approx_corr.hs
  • Run approx_corr.exe
  • Enter the maximal n
  • The result of each computation is displayed, as well as the cumulative real time (in ms) up to that point.

Code:

import Data.List
import Data.Time
import Data.Time.Clock.POSIX

num2bin :: Int -> Int -> [Int]
num2bin 0 _ = []
num2bin n k| k >= 2^(n-1) = 1 : num2bin (n-1)( k-2^(n-1))
           | otherwise  = 0: num2bin (n-1) k

genBinNum :: Int -> [[Int]]
genBinNum n = map (num2bin n) [0..2^n-1]

pairs :: [a] -> [([a],[a])]
pairs xs = zip (prefixes xs) (suffixes xs)
   where prefixes = tail . init . inits 
         suffixes = map reverse . prefixes . reverse 

hammingDist :: (Num b, Eq a) => ([a],[a]) -> b     
hammingDist (a,b) = sum $ zipWith (\u v -> if u /= v then 1 else 0) a b

f :: Int -> Int
f n = length $ nub $ map (map ((<=1).hammingDist) . pairs) $ genBinNum n
--f n = sum [1..n]

--time in milliseconds
getTime = getCurrentTime >>= pure . (1000*) . utcTimeToPOSIXSeconds >>= pure . round


main :: IO()
main = do 
    maxns <- getLine 
    let maxn = (read maxns)::Int
    t0 <- getTime 
    loop 1 maxn t0
     where loop n maxn t0|n==maxn = return ()
           loop n maxn t0
             = do 
                 putStrLn $ "fun eval: " ++ (show n) ++ ", " ++ (show $ (f n)) 
                 t <- getTime
                 putStrLn $ "time: " ++ show (t-t0); 
                 loop (n+1) maxn t0
\$\endgroup\$
  • \$\begingroup\$ The code appears not to give output as it runs. This makes it a little hard to test. \$\endgroup\$ – user9206 Jun 12 '17 at 8:48
  • \$\begingroup\$ Strange, does it compile without error? What happens if you try to compile the program main = putStrLn "Hello World!" ? \$\endgroup\$ – flawr Jun 12 '17 at 8:58
  • \$\begingroup\$ The Data.Bits module might be useful. For your main loop, you could use something like main = do maxn <- getmax; t0 <- gettime; loop 1 where loop n|n==maxn = return () and loop n = do printresult n (f n); t <- gettime; printtime (t-t0); loop (n+1). getmax could for example use getArgs to use the program arguments. \$\endgroup\$ – Christian Sievers Jun 12 '17 at 10:23
  • \$\begingroup\$ @ChristianSievers Thanks a lot!!! I asked this question at stackoverflow, I think it would be great if you could add that there too! \$\endgroup\$ – flawr Jun 12 '17 at 11:52
  • \$\begingroup\$ I don't see how to answer there. You have a similar loop there already, and I didn't say anything about getting the time: that you already had here. \$\endgroup\$ – Christian Sievers Jun 12 '17 at 14:34
1
\$\begingroup\$

A Haskell solution, using popCount and manually managed parallelism

Compile: ghc -rtsopts -threaded -O2 -fllvm -Wall foo.hs

(drop the -llvm if it doesn't work)

Run : ./foo +RTS -N

module Main (main) where

import Data.Bits
import Data.Word
import Data.List
import qualified Data.IntSet as S 
import System.IO
import Control.Monad
import Control.Concurrent
import Control.Exception.Base (evaluate)

pairs' :: Int -> Word64 -> Int
pairs' n s = fromIntegral $ foldl' (.|.) 0 $ map mk [1..n]
  where mk d = let mask = 1 `shiftL` d - 1 
                   pc = popCount $! xor (s `shiftR` (n - d)) (s .&. mask)
               in  if pc <= 1 
                     then mask + 1 
                     else 0 

mkSet :: Int -> Word64 -> Word64 -> S.IntSet
mkSet n a b = S.fromList $ map (pairs' n) [a .. b]

f :: Int -> IO Int
f n 
   | n < 4 = return $ S.size $ mkSet n 0 mxbound
   | otherwise = do
        mvs <- replicateM 4 newEmptyMVar
        forM_ (zip mvs cpairs) $ \(mv,(mi,ma)) -> forkIO $ do
          evaluate (mkSet n mi ma) >>= putMVar mv
        set <- foldl' S.union S.empty <$> mapM readMVar mvs
        return $! S.size set
   where
     mxbound = 1 `shiftL` (n - 1)
     bounds = [0,1 `shiftL` (n - 3) .. mxbound]
     cpairs = zip bounds (drop 1 bounds)

main :: IO()
main = do
    hSetBuffering stdout LineBuffering
    mapM_ (f >=> print) [1..]
\$\endgroup\$
  • \$\begingroup\$ There is a buffering problem it seems in that I don't get any output at all if I run it from the cygwim command line. \$\endgroup\$ – user9206 Jun 12 '17 at 11:18
  • \$\begingroup\$ I just updated my solution, but I don't know if it will help much. \$\endgroup\$ – bartavelle Jun 12 '17 at 11:22
  • \$\begingroup\$ @Lembik Unsure if that is obvious, but that should be compiled with -O3, and might be faster with -O3 -fllvm ... \$\endgroup\$ – bartavelle Jun 12 '17 at 11:26
  • \$\begingroup\$ (And all build files should be removed before recompiling, if not source code change happened) \$\endgroup\$ – bartavelle Jun 12 '17 at 11:28
  • \$\begingroup\$ @Lembik : I introduced parallelism. It should be a bit faster. \$\endgroup\$ – bartavelle Jun 12 '17 at 12:17
0
\$\begingroup\$

Python 2 + pypy, n = 22

Here is a really simple Python solution as a sort of baseline benchmark.

import itertools
def hamming(A, B):
    n = len(A)
    assert(len(B) == n)
    return n-sum([A[i] == B[i] for i in xrange(n)])

def prefsufflist(P):
    n = len(P)
    return [hamming(P[:i], P[n-i:n]) for i in xrange(1,n+1)]

bound = 1
for n in xrange(1,25):
    booleans = set()
    for P in itertools.product([0,1], repeat = n):
        booleans.add(tuple(int(HD <= bound) for HD in prefsufflist(P)))
    print "n = ", n, len(booleans)
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy