Betrachten Sie eine binäre Zeichenfolge Svon Länge n. Indizieren 1wir von , können wir die Hamming-Entfernungen zwischen S[1..i+1]und S[n-i..n]für alle iin der Reihenfolge von 0bis berechnen n-1. Der Hamming-Abstand zwischen zwei Saiten gleicher Länge ist die Anzahl der Positionen, an denen sich die entsprechenden Symbole unterscheiden. Beispielsweise,

S = 01010

gibt

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

Dies liegt daran, dass 0Übereinstimmungen 0, 01Hamming-Distanz zwei zu 10, 010Übereinstimmungen 010, 0101Hamming-Distanz vier zu 1010 und schließlich sich 01010selbst entsprechen.

Wir interessieren uns jedoch nur für Ausgaben, bei denen der Hamming-Abstand höchstens 1 beträgt. In dieser Aufgabe werden wir also eine Yangeben, wenn die Hamming-Distanz höchstens eins beträgt, und eine Nandere. Also in unserem obigen Beispiel würden wir bekommen

[Y, N, Y, N, Y]

Definieren Sie f(n)die Anzahl der verschiedenen Arrays von Ys und Ns, die Sie erhalten, wenn Sie alle 2^nverschiedenen möglichen Bitfolgen Sder Länge durchlaufen n.

Aufgabe

Zum Erhöhen nab 1sollte Ihr Code ausgegeben werden f(n).

Beispielantworten

Denn n = 1..24die richtigen Antworten sind:

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

Wertung

Ihr Code sollte iterieren, nachdem Sie n = 1die Antwort für jede einzelne nnacheinander gegeben haben. Ich werde den gesamten Lauf zeitlich festlegen und ihn nach zwei Minuten töten.

Ihre Punktzahl ist die höchste, die nSie in dieser Zeit erreichen.

Bei einem Gleichstand gewinnt die erste Antwort.

Wo wird mein Code getestet?

Ich werde Ihren Code auf meinem (etwas alten) Windows 7-Laptop unter Cygwin ausführen. Bitte geben Sie daher jede mögliche Unterstützung, um dies zu vereinfachen.

Mein Laptop hat 8 GB RAM und eine Intel i7 5600U @ 2,6 GHz (Broadwell) -CPU mit 2 Kernen und 4 Threads. Der Befehlssatz umfasst SSE4.2, AVX, AVX2, FMA3 und TSX.

Führende Einträge pro Sprache

• n = 40 in Rust mit CryptoMiniSat von Anders Kaseorg. (In Lubuntu Gast-VM unter Vbox.)
• n = 35 in C ++ unter Verwendung der BuDDy-Bibliothek von Christian Seviers. (In Lubuntu Gast-VM unter Vbox.)
• n = 34 in Clingo von Anders Kaseorg. (In Lubuntu Gast-VM unter Vbox.)
• n = 31 in Rust von Anders Kaseorg.
• n = 29 in Clojure von NikoNyrh.
• n = 29 in C von Bartavelle.
• n = 27 in Haskell von Bartavelle
• n = 24 in Pari / gp von Alephalpha.
• n = 22 in Python 2 + pypy von mir.
• n = 21 in Mathematica von Alephalpha. (Selbst berichtet)

Zukünftige Kopfgelder

Ich werde jetzt eine Prämie von 200 Punkten für jede Antwort geben, die in zwei Minuten auf meinem Computer auf n = 80 steigt .

Rust + CryptoMiniSat , Nr. 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"

Wie es funktioniert

Dies führt eine rekursive Suche durch den Baum aller Teilzuweisungen zu Präfixen des Y / N-Arrays durch, wobei bei jedem Schritt mit einem SAT-Solver geprüft wird, ob die aktuelle Teilzuweisung konsistent ist, und wenn nicht, wird die Suche beschnitten. CryptoMiniSat ist aufgrund seiner speziellen Optimierungen für XOR-Klauseln der richtige SAT-Löser für diesen Job.

Die drei Familien von Einschränkungen sind

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

mit der Ausnahme, dass S ₀ als eine Optimierung zu falsch gezwungen wird, so dass D _{k 0} einfach gleich S _{k ist} .

Rust, n ≈ 30 oder 31 oder 32

Auf meinem Laptop (zwei Kerne, i5-6200U) reicht dies für n = 1, ..., 31 in 53 Sekunden und etwa 2,5 GB Speicher oder für n = 1, ..., 32 in 105 Sekunden und etwa 5 GB der Erinnerung. Kompilieren Sie mit cargo build --releaseund führen Sie es aus 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"

Probieren Sie es online!

Ich habe auch eine etwas langsamere Variante mit sehr viel weniger Speicher.

C ++ mit der BuDDy-Bibliothek

Ein anderer Ansatz: Verwenden Sie eine Binärformel (als binäres Entscheidungsdiagramm ), die die Bits von Sals Eingabe verwendet und true ist, wenn dies einige feste Werte von Yoder Nan bestimmten ausgewählten Positionen ergibt . Wenn diese Formel nicht konstant falsch ist, wählen Sie eine freie Position aus und versuchen Sie es mit beiden Yund N. Wenn keine freie Position vorhanden ist, haben wir einen möglichen Ausgabewert gefunden. Wenn die Formel konstant falsch ist, wird zurückverfolgt.

Dies funktioniert relativ vernünftig, da es so wenige mögliche Werte gibt, dass wir oft früh zurückgehen können. Ich habe eine ähnliche Idee mit einem SAT-Solver versucht, aber das war weniger erfolgreich.

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

Um mit Debian 8 (Jessie) zu kompilieren, installieren Sie libbdd-devund tun Sie g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. Es kann nützlich sein, das erste Argument auf bdd_initnoch mehr zu erhöhen .

Clingo, Nr. 30 oder 31 34

Ich war ein wenig überrascht, als ich sah, dass fünf Zeilen Clingo-Code meine Brute-Force-Rust-Lösung überholten und Christians BuDDy-Lösung sehr nahe kamen - es sieht so aus, als würde dies auch mit einem höheren Zeitlimit übertroffen.

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

Standardmäßig begrenzt Pari / GP die Stapelgröße auf 8 MB. In der ersten Codezeile wird default(parisize, "4g")dieses Limit auf 4 GB festgelegt. Wenn es immer noch einen Stapelüberlauf gibt, können Sie ihn auf 8 GB einstellen.

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 in 75 38 Sekunden, 30 in 80 und 31 in 165

Die Laufzeit von Intel i7 6700K beträgt weniger als 200 MB.

project.clj (verwendet com.climate / claypoole für 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)

Quellcode:

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

Als Brute-Force-Lösung durchläuft jeder Thread eine Teilmenge des Bereichs (2 ^ 12 Elemente) und erstellt eine Reihe von ganzzahligen Werten, die auf erkannte Muster hinweisen. Diese werden dann zusammen "unionisiert" und somit die eindeutige Zählung berechnet. Ich hoffe, der Code ist nicht zu knifflig, obwohl er Threading-Makros verwendet . Ich mainführe den Test einige Male aus, um JVM zum Aufwärmen zu bringen.

Update: Wenn Sie nur die Hälfte der ganzen Zahlen durchlaufen, wird aufgrund der Symmetrie dasselbe Ergebnis erzielt. Überspringen Sie auch Zahlen mit einer höheren Bitanzahl in der unteren Hälfte der Zahl, da sie ebenfalls Duplikate erzeugen.

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

Ergebnisse auf verschiedenen Hardwares, erwartete Laufzeit ist 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

Sie können diesen Single-Thread einfach erstellen und diese Abhängigkeit von Drittanbietern vermeiden, indem Sie den Standard für Folgendes verwenden:

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

Nun, die eingebaute pmap existiert auch, aber Claypoole hat mehr Funktionen und Abstimmbarkeit.

Mathematica, n = 19

drücke alt +. abbrechen und das Ergebnis wird gedruckt

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

f [n_]: = Länge @
     DeleteDuplicates @
      Transponieren@
       Tabelle [2> Tr @ IntegerDigits [#, 2] & / @ 
         BitX oder [BitShiftRight [#, n - i], Mod [#, 2 ^ i], {i, 1, 
         n - 1}] & @ Range [0, 2 ^ (n - 1)];
Mach [Print [n -> f @ n], {n, Infinity}]

Zum Vergleich Jenny_mathy Antwort gibt n = 19auf meinem Computer.

Der langsamste Teil ist Tr@IntegerDigits[#, 2] &. Es ist eine Schande, dass Mathematica kein eingebautes Hamming-Gewicht hat.

Wenn Sie meinen Code testen möchten, können Sie eine kostenlose Testversion von Mathematica herunterladen .

AC-Version mit eingebautem Popcount

Funktioniert besser mit clang -O3, funktioniert aber auch, wenn Sie nur haben 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, (inoffizielles n = 20)

Dies ist nur der naive Ansatz - bisher ohne Optimierungen. Ich fragte mich, wie gut es gegen andere Sprachen abschneiden würde.

Wie man es benutzt (vorausgesetzt, Sie haben eine Hashell-Plattform) installiert):

• Fügen Sie den Code in eine Datei ein approx_corr.hs (oder einen anderen Namen, ändern Sie die folgenden Schritte entsprechend)
• Navigieren Sie zu der Datei und führen Sie sie aus ghc approx_corr.hs
• Lauf approx_corr.exe
• Geben Sie das Maximum ein n
• Das Ergebnis jeder Berechnung sowie die kumulierte Echtzeit (in ms) bis zu diesem Punkt werden angezeigt.

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

Eine Haskell-Lösung, die popCount und manuell verwaltete Parallelität verwendet

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

(Lass fallen -llvm wenn es nicht funktioniert)

Lauf : ./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

Hier ist eine wirklich einfache Python-Lösung als eine Art Referenz.

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)