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\$
4
  • \$\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\$ Jun 4, 2017 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, 2017 at 8:16
  • \$\begingroup\$ @Lembik Do we measure CPU time or real time? \$\endgroup\$
    – flawr
    Jun 11, 2017 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, 2017 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\$
17
  • 2
    \$\begingroup\$ Woohoooooo ! :) \$\endgroup\$
    – user9206
    Jun 12, 2017 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, 2017 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\$ Jun 12, 2017 at 21:26
  • 1
    \$\begingroup\$ @Lembik I'm getting a 404 from that paste. \$\endgroup\$
    – user45941
    Jun 14, 2017 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\$ Jun 16, 2017 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\$
9
  • \$\begingroup\$ What optimisations have you used? \$\endgroup\$
    – user9206
    Jun 5, 2017 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\$ Jun 5, 2017 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, 2017 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\$ Jun 6, 2017 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\$ Jun 10, 2017 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\$
9
  • \$\begingroup\$ This looks interesting. What do you get to with this? \$\endgroup\$
    – user9206
    Jun 11, 2017 at 18:24
  • \$\begingroup\$ @Lembik I get 31 in 100s on very old hardware that won't let me answer faster \$\endgroup\$ Jun 11, 2017 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, 2017 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\$ Jun 11, 2017 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\$ Jun 12, 2017 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\$
8
  • \$\begingroup\$ This is great! It looks from your graph that the BuDDy solution suddenly gets worse. Any idea why? \$\endgroup\$
    – user9206
    Jun 12, 2017 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\$ Jun 12, 2017 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\$ Jun 12, 2017 at 9:51
  • 2
    \$\begingroup\$ I do love your graph. \$\endgroup\$
    – user9206
    Jun 12, 2017 at 17:22
  • 1
    \$\begingroup\$ You get a shocking performance boost by using the option --configuration=crafty (jumpy and trendy give similar results). \$\endgroup\$ Nov 10, 2019 at 23:10
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\$
2
  • \$\begingroup\$ Reaches 22 and then gives a stackoverflow. \$\endgroup\$
    – user9206
    Jun 4, 2017 at 17:29
  • \$\begingroup\$ Gets to 24 now. \$\endgroup\$
    – user9206
    Jun 5, 2017 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\$
3
  • \$\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, 2017 at 15:32
  • \$\begingroup\$ Sadly it's a no for 30. Elapsed time: 217150.87386 msecs \$\endgroup\$
    – user9206
    Jun 6, 2017 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, 2017 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\$
3
  • \$\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, 2017 at 3:57
  • \$\begingroup\$ Can you show the output of the code? \$\endgroup\$
    – user9206
    Jun 4, 2017 at 4:17
  • 1
    \$\begingroup\$ I meant f(n)... fixed \$\endgroup\$
    – ZaMoC
    Jun 4, 2017 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\$
3
  • \$\begingroup\$ It gets to 24 very quickly and then ends. You need to increase the limit. \$\endgroup\$
    – user9206
    Jun 12, 2017 at 13:16
  • \$\begingroup\$ Oh god, I forgot to remove the benchmark code! I'll remove the two offending lines :/ \$\endgroup\$
    – bartavelle
    Jun 12, 2017 at 13:18
  • \$\begingroup\$ @Lembik should be fixed now \$\endgroup\$
    – bartavelle
    Jun 12, 2017 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\$
7
  • \$\begingroup\$ The code appears not to give output as it runs. This makes it a little hard to test. \$\endgroup\$
    – user9206
    Jun 12, 2017 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, 2017 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\$ Jun 12, 2017 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, 2017 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\$ Jun 12, 2017 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\$
10
  • \$\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, 2017 at 11:18
  • \$\begingroup\$ I just updated my solution, but I don't know if it will help much. \$\endgroup\$
    – bartavelle
    Jun 12, 2017 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, 2017 at 11:26
  • \$\begingroup\$ (And all build files should be removed before recompiling, if not source code change happened) \$\endgroup\$
    – bartavelle
    Jun 12, 2017 at 11:28
  • \$\begingroup\$ @Lembik : I introduced parallelism. It should be a bit faster. \$\endgroup\$
    – bartavelle
    Jun 12, 2017 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