Rozważ ciąg binarny So długości n. Indeksując od 1, możemy obliczyć odległości Hamminga pomiędzy S[1..i+1]i S[n-i..n]dla wszystkich iw kolejności od 0do n-1. Odległość Hamminga między dwoma strunami o równej długości jest liczbą pozycji, w których odpowiednie symbole są różne. Na przykład,

S = 01010

daje

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

Wynika to z tego, że 0mecze 0, 01odległość Hamminga dwa do 10, 010mecze 010, 0101odległość Hamminga cztery do, 1010 a ostatecznie 01010dopasowuje się.

Interesują nas jednak tylko wyjścia, w których odległość Hamminga wynosi najwyżej 1. W tym zadaniu Yzgłosimy, czy odległość Hamminga wynosi co najwyżej jeden, a w Nprzeciwnym razie. W powyższym przykładzie otrzymalibyśmy

[Y, N, Y, N, Y]

Zdefiniuj f(n)liczbę odrębnych tablic Ys i Ns, które otrzymujesz podczas iteracji wszystkich 2^nróżnych możliwych ciągów bitów Sdługości n.

Zadanie

Aby zwiększyć, nzaczynając od 1, twój kod powinien wypisywać f(n).

Przykładowe odpowiedzi

Dla n = 1..24, poprawne odpowiedzi są:

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

Punktacja

Twój kod powinien powtarzać się od n = 1udzielenia odpowiedzi po nkolei. Poświęcę cały bieg, zabijając go po dwóch minutach.

Twój wynik jest najwyższy nw tym czasie.

W przypadku remisu pierwsza odpowiedź wygrywa.

Gdzie będzie testowany mój kod?

Uruchomię twój kod na moim (nieco starym) laptopie z systemem Windows 7 pod cygwin. W związku z tym proszę udzielić wszelkiej możliwej pomocy, aby ułatwić to.

Mój laptop ma 8 GB pamięci RAM i procesor Intel i7 [email protected] GHz (Broadwell) z 2 rdzeniami i 4 wątkami. Zestaw instrukcji obejmuje SSE4.2, AVX, AVX2, FMA3 i TSX.

Wiodące wpisy według języka

• n = 40 w Rust za pomocą CryptoMiniSat, autor: Anders Kaseorg. (W maszynie wirtualnej gościa Lubuntu pod Vbox.)
• n = 35 w C ++ przy użyciu biblioteki BuDDy, Christian Seviers. (W maszynie wirtualnej gościa Lubuntu pod Vbox.)
• n = 34 w Clingo autorstwa Andersa Kaseorga. (W maszynie wirtualnej gościa Lubuntu pod Vbox.)
• n = 31 w Rust przez Andersa Kaseorga.
• n = 29 w Clojure autorstwa NikoNyrh.
• n = 29 w C według bartavelle.
• n = 27 w Haskell według bartavelle
• n = 24 w Pari / gp przez alephalpha.
• n = 22 w Python 2 + pypy przeze mnie.
• n = 21 w Mathematica przez alephalpha. (Własne zgłoszenie)

Przyszłe nagrody

Dam teraz nagrodę w wysokości 200 punktów za każdą odpowiedź, która na mojej maszynie osiągnie n = 80 w ciągu dwóch minut.

Rdza + 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 <[email protected]>"]

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

Jak to działa

Dokonuje to rekurencyjnego przeszukiwania drzewa wszystkich częściowych przypisań do przedrostków tablicy Y / N, używając solwera SAT do sprawdzania na każdym kroku, czy bieżące częściowe przypisanie jest spójne i przycina wyszukiwanie, jeśli nie. CryptoMiniSat jest właściwym narzędziem SAT do tego zadania dzięki specjalnej optymalizacji klauzul XOR.

Trzy rodziny ograniczeń to

S _i ⊕ S _{k + i} ⊕ D _ki , dla 1 ≤ k ≤ n - 2, 0 ≤ i ≤ n - k ;
D _{ki 1} ∨ D _{ki 2} ∨ ¬ A _k , dla 1 ≤ k ≤ n - 2, 0 ≤ i ₁ < i ₂ ≤ n - k ;
¬ D _{ki 1} ∨ ⋯ ∨ ¬ D _{ki n - k - 1}∨ A _k , dla 1 ≤ k ≤ n - 2, 0 ≤ i ₁ <⋯ < i _{n - k - 1} ≤ n - k ;

z wyjątkiem tego, że jako optymalizacja S ₀ jest zmuszone do fałszu, tak że D _{k 0} jest po prostu równe S _k .

Rdza, n ≈ 30 lub 31 lub 32

Na moim laptopie (dwa rdzenie, i5-6200U), to przechodzi przez n = 1,…, 31 w 53 sekundach, używając około 2,5 GiB pamięci lub przez n = 1,…, 32 w 105 sekund, używając około 5 GiB pamięciowy. Kompiluj cargo build --releasei uruchamiaj 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 <[email protected]>"]

[dependencies]
rayon = "0.7.0"

Wypróbuj online!

Mam też nieco wolniejszy wariant, który wykorzystuje znacznie mniej pamięci.

C ++ przy użyciu biblioteki BuDDy

Inne podejście: zastosuj wzór binarny (jako binarny schemat decyzyjny ), który przyjmuje bity Sjako dane wejściowe i jest prawdziwy iff, który daje pewne ustalone wartości Ylub Nw niektórych wybranych pozycjach. Jeśli ta formuła nie jest stała fałszywe, wybierz wolną pozycję i rekursja, próbując zarówno Yi N. Jeśli nie ma wolnej pozycji, znaleźliśmy możliwą wartość wyjściową. Jeśli formuła ma stałą wartość false, cofnij się.

Działa to względnie rozsądnie, ponieważ istnieje tak mało możliwych wartości, że często możemy wcześnie wycofać się. Wypróbowałem podobny pomysł z solwerem SAT, ale to się nie powiodło.

#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" ;
  }
}

Aby skompilować z debianem 8 (jessie), zainstaluj libbdd-devi wykonaj g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. Przydatne może być zwiększenie pierwszego argumentu do bdd_initjeszcze większego.

Clingo, n ≈ 30 lub 31 34

Byłem trochę zaskoczony, gdy pięć linii kodu Clingo wyprzedziło moje brutalne rozwiązanie Rust i zbliżyłem się bardzo do rozwiązania BuDDy Christiana - wygląda na to, że to też by go pobiło z wyższym limitem czasowym.

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

Pari / GP , 23

Domyślnie Pari / GP ogranicza swój stos do 8 MB. Pierwszy wiersz kodu default(parisize, "4g")określa ten limit na 4 GB. Jeśli nadal daje przepływ stosu, możesz ustawić go na 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)))

Clojure, 29 w 75 38 sekund, 30 w 80 i 31 w 165

Runtimes z Intel i7 6700K , zużycie pamięci jest mniejsze niż 200 MB.

project.clj (używa com.climate / claypoole do wielowątkowości):

(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)

Kod źródłowy:

(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))

Rozwiązanie brute-force, każdy wątek przekracza podzbiór zakresu (2 ^ 12 elementów) i buduje zestaw wartości całkowitych wskazujących wykryte wzorce. Są one następnie „łączone” razem, a zatem obliczana jest odrębna liczba. Mam nadzieję, że kod nie jest zbyt trudny do śledzenia, mimo że używa często makr wątków . Mój maintest przeprowadza kilka razy, aby rozgrzać JVM.

Aktualizacja: Iteracja tylko połowy liczb całkowitych daje ten sam wynik dzięki symetrii. Pomijanie liczb o większej liczbie bitów również w dolnej połowie liczby, ponieważ produkują również duplikaty.

Wbudowany 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]

Wyniki dla różnych programów sprzętowych, oczekiwany czas działania to 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

Możesz łatwo uczynić to jedno-wątkowym i uniknąć zależności zewnętrznych, używając standardu dla:

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

Cóż, wbudowane pmap również istnieje, ale claypoole ma więcej funkcji i możliwości dostrajania.

Mathematica, n = 19

naciśnij alt +. przerwać, a wynik zostanie wydrukowany

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]

Mathematica, 21 lat

f [n_]: = długość @
     DeleteDuplicates @
      Transponować@
       Tabela [2> Tr @ IntegerDigits [#, 2] & / @ 
         BitXor [BitShiftRight [#, n - i], Mod [#, 2 ^ i]], {i, 1, 
         n - 1}] & @ Range [0, 2 ^ (n - 1)];
Wykonaj [Drukuj [n -> f @ n], {n, Infinity}]

Dla porównania odpowiedź Jenny_mathy podaje się n = 19na moim komputerze.

Najwolniejsza część to Tr@IntegerDigits[#, 2] &. Szkoda, że ​​Mathematica nie ma wbudowanej wagi Hamminga.

Jeśli chcesz przetestować mój kod, możesz pobrać bezpłatną wersję próbną Mathematica .

Wersja AC z wbudowanym popcount

Działa lepiej clang -O3, ale działa również, jeśli masz 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;
}

Haskell, (nieoficjalne n = 20)

To tylko naiwne podejście - jak dotąd bez żadnych optymalizacji. Zastanawiałem się, jak dobrze poradzi sobie z innymi językami.

Jak z niego korzystać (zakładając, że masz zainstalowaną platformę Haskell ):

• Wklej kod do jednego pliku approx_corr.hs(lub innej nazwy, odpowiednio zmodyfikuj następujące kroki)
• Przejdź do pliku i uruchom ghc approx_corr.hs
• Biegać approx_corr.exe
• Wpisz maksimum n
• Wyświetlany jest wynik każdego obliczenia, a także skumulowany czas rzeczywisty (w ms) do tego momentu.

Kod:

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

Rozwiązanie Haskell, wykorzystujące popCount i ręcznie zarządzaną równoległość

Skompilować: ghc -rtsopts -threaded -O2 -fllvm -Wall foo.hs

(upuść, -llvmjeśli to nie działa)

Biegać : ./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..]

Python 2 + pypy, n = 22

Oto naprawdę proste rozwiązanie Python jako rodzaj podstawowego testu porównawczego.

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)