Localizando correlações aproximadas

14

Considere uma cadeia Sde comprimento binária n. Indexando de 1, podemos calcular as distâncias de Hamming entre S[1..i+1]e S[n-i..n]para todos ina ordem de 0para n-1. A distância de Hamming entre duas cordas de igual comprimento é o número de posições nas quais os símbolos correspondentes são diferentes. Por exemplo,

S = 01010

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

Isso ocorre porque 0combinações 0, 01tem distância de Hamming 2 a 10, 010combina 010, 0101tem distância de Hamming 4 a 1010 e, finalmente, 01010combina-se.

No entanto, estamos interessados ​​apenas em saídas onde a distância de Hamming é no máximo 1. Portanto, nesta tarefa, reportaremos a Yse a distância de Hamming é no máximo uma e Noutra. Então, no nosso exemplo acima, teríamos

[Y, N, Y, N, Y]

Define f(n)-se o número de matrizes distintas de Ys e Ns que se tem quando iteração sobre todas as 2^nsequências de bits diferentes possíveis Sde comprimento n.

Tarefa

Para aumentar a npartir de 1, seu código deve ser exibido f(n).

Respostas de exemplo

Pois n = 1..24, as respostas corretas são:

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

Pontuação

Seu código deve repetir n = 1a resposta de cada um n. Eu cronometrarei a corrida inteira, matando-a depois de dois minutos.

Sua pontuação é a mais alta que nvocê obtém nesse período.

Em caso de empate, a primeira resposta vence.

Onde meu código será testado?

Vou executar seu código no meu laptop Windows 7 (um pouco antigo) sob o cygwin. Como resultado, forneça qualquer assistência possível para facilitar isso.

Meu laptop possui 8 GB de RAM e uma CPU Intel i7 [email protected] GHz (Broadwell) com 2 núcleos e 4 threads. O conjunto de instruções inclui SSE4.2, AVX, AVX2, FMA3 e TSX.

Entradas principais por idioma

  • n = 40 em Rust usando CryptoMiniSat, de Anders Kaseorg. (Na VM convidada do Lubuntu, no Vbox.)
  • n = 35 em C ++ usando a biblioteca BuDDy, de Christian Seviers. (Na VM convidada do Lubuntu, no Vbox.)
  • n = 34 em Clingo por Anders Kaseorg. (Na VM convidada do Lubuntu, no Vbox.)
  • n = 31 em Rust por Anders Kaseorg.
  • n = 29 in Clojure por NikoNyrh.
  • n = 29 in C por bartavelle.
  • n = 27 in Haskell por bartavelle
  • n = 24 in Pari / gp por alefalpha.
  • n = 22 em Python 2 + pypy por mim.
  • n = 21 no Mathematica por alefalpha. (Auto-relatado)

Recompensas futuras

Agora, darei uma recompensa de 200 pontos por qualquer resposta que chegue a n = 80 na minha máquina em dois minutos.


fonte
Você conhece algum truque que permitirá que alguém encontre um algoritmo mais rápido que uma força bruta ingênua? Caso contrário, esse desafio é "implemente isso em x86" (ou talvez se conheçamos sua GPU ...).
Jonathan Allan
@ JonathanAllan Certamente é possível acelerar uma abordagem muito ingênua. Exatamente o quão rápido você pode chegar, não tenho certeza. Curiosamente, se mudarmos a questão para que você obtenha um Y se a distância de Hamming for no máximo 0 e um N de outra forma, então haverá uma fórmula de forma fechada conhecida.
@Lembik Medimos o tempo da CPU ou em tempo real?
flawr
@flawr Estou medindo o tempo real, mas executando-o algumas vezes e tomando o mínimo para eliminar esquisitices.

Respostas:

9

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

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

Como funciona

Isso faz uma pesquisa recursiva na árvore de todas as atribuições parciais aos prefixos da matriz Y / N, usando um solucionador SAT para verificar a cada passo se a atribuição parcial atual é consistente e poda a pesquisa, se não. O CryptoMiniSat é o solucionador SAT certo para este trabalho devido às suas otimizações especiais para cláusulas XOR.

As três famílias de restrições são

S iS k + iD ki , para 1 ≤ kn - 2, 0 ≤ i ≤ n - k ;
D ki 1D ki 2 ¬ A k , para 1 ≤ kn - 2, 0 ≤ i 1 < i 2 n - k ; ¬ D ki 1 ¬ ⋯ ∨ D ki n - kA k , para 1 ≤ k
- 1∨ ≤ n - 2, 0 ≤ i 1 <⋯ < i n - k - 1n - k ;

exceto que, como uma otimização, S 0 é forçado a falso, de modo que D k 0 é simplesmente igual a S k .

Anders Kaseorg
fonte
2
Woohoooooo! :)
Ainda estou tentando compilar isso no Windows (usando o cygwin + gcc). Eu clonei o cryptominisat e o compilei. Mas ainda não sei como compilar o código de ferrugem. Quando eu cargo buildrecebo--- stderr CMake Error: Could not create named generator Visual Studio 14 2015 Win64
2
@ rahnema1 Obrigado, mas parece que o problema está no sistema de compilação CMake da biblioteca C ++ incorporada na caixa cryptominisat, não na própria Rust.
Anders Kaseorg
1
@Embembik Estou recebendo um 404 dessa pasta.
Mego 14/06
1
@ChristianSievers Boa pergunta. Isso funciona, mas parece ser um pouco mais lento (2 × mais ou menos). Não sei por que não deve ser tão bom, então talvez o CryptoMiniSat não tenha sido bem otimizado para esse tipo de carga de trabalho incremental.
Anders Kaseorg
9

Ferrugem, n ≈ 30 ou 31 ou 32

No meu laptop (dois núcleos, i5-6200U), isso passa por n = 1,…, 31 em 53 segundos, usando cerca de 2,5 GiB de memória, ou por n = 1,…, 32 em 105 segundos, usando cerca de 5 GiB de memória. Compile cargo build --releasee execute 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"

Experimente online!

Eu também tenho uma variante um pouco mais lenta, usando muito menos memória.

Anders Kaseorg
fonte
Quais otimizações você usou?
1
@Lembik A maior otimização, além de fazer tudo com aritmética bit a bit em uma linguagem compilada, é usar apenas o não determinismo necessário para definir um prefixo da matriz Y / N. Faço uma pesquisa recursiva em possíveis prefixos da matriz Y / N, levando um vetor de possíveis strings que atingem esse prefixo, mas apenas as strings cujo meio não examinado é preenchido com zeros. Dito isto, essa ainda é uma pesquisa exponencial e essas otimizações apenas a aceleram por fatores polinomiais.
Anders Kaseorg
É uma boa resposta. Obrigado. Espero que alguém se dedique à combinatória para obter uma velocidade significativa.
@ Lembik Corrigi um bug de perda de memória, fiz mais micro-otimização e adicionei paralelismo. Teste novamente quando tiver uma chance - espero aumentar minha pontuação em 1 ou 2. Você tem ideias combinatórias em mente para acelerar maiores? Eu não inventei nada.
Anders Kaseorg
1
@Lembik Não existe uma fórmula fornecida na entrada OEIS. (O código do Mathematica também parece usar força bruta.) Se você souber de um, poderá contar a eles sobre isso.
Christian Sievers
6

C ++ usando a biblioteca BuDDy

Uma abordagem diferente: tenha uma fórmula binária (como diagrama de decisão binário ) que aceite os bits Scomo entrada e seja verdadeira se der alguns valores fixos de You Nem determinadas posições selecionadas. Se essa fórmula não for constante false, selecione uma posição livre e recorra, tentando ambasY e N. Se não houver posição livre, encontramos um possível valor de saída. Se a fórmula for constante false, retorne.

Isso funciona relativamente razoável, porque há tão poucos valores possíveis, para que possamos voltar sempre cedo. Tentei uma idéia semelhante com um solucionador SAT, mas isso foi menos bem-sucedido.

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

Para compilar com o debian 8 (jessie), instale libbdd-deve faça g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. Pode ser útil aumentar bdd_initainda mais o primeiro argumento .

Peneiradores cristãos
fonte
Isso parece interessante. O que você faz com isso?
@Lembik Recebo 31 em 100s em hardware muito antigo que não vai deixar-me responder mais rápido
Christian Sievers
Qualquer ajuda que você possa dar sobre como compilar isso no Windows (por exemplo, usando o cygwin) recebida com gratidão.
@Lembik Eu não sei sobre Windws mas github.com/fd00/yacp/tree/master/buddy parece wrt útil cygwin
Christian Sievers
1
Uau, ok, você me convenceu de que preciso adicionar esta biblioteca ao meu kit de ferramentas. Bem feito!
Anders Kaseorg
4

Clingo, n. 30 ou 31 34

Fiquei um pouco surpreso ao ver cinco linhas de código Clingo ultrapassando minha solução Rust de força bruta e chegando muito perto da solução BuDDy de Christian - parece que isso também seria melhor com um limite de tempo mais alto.

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

enredo

Anders Kaseorg
fonte
Isso é ótimo! Parece no seu gráfico que a solução BuDDy piora subitamente. Alguma idéia do porquê?
@ Lembik Eu não investiguei o BuDDy o suficiente para ter certeza, mas talvez ele fique sem cache nesse momento?
Anders Kaseorg 12/06
Uau! Eu acho que um primeiro valor mais alto bdd_initpode ajudar, ou permitir aumentar mais a tabela de nós chamando bdd_setmaxincreasecom um valor muito acima do padrão de 50000. - Você está usando a versão alterada do meu programa?
Christian Sievers
2
Eu amo o seu gráfico.
1
Você recebe um aumento chocante no desempenho usando a opção --configuration=crafty( jumpye trendyfornece resultados semelhantes).
Christian Sievers
2

Pari / GP , 23

Por padrão, o Pari / GP limita o tamanho da pilha a 8 MB. A primeira linha do código default(parisize, "4g"),, define esse limite para 4 GB. Se ele ainda der um fluxo de pilha, você poderá configurá-lo para 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)))
alefalpha
fonte
Atinge 22 e, em seguida, fornece um fluxo de pilha.
Chega a 24 agora.
2

Clojure, 29 em 75 38 segundos, 30 em 80 e 31 em 165

Duração de tempos de execução Intel i7 6700K , o uso de memória é inferior a 200 MB.

project.clj (usa com.climate / claypoole para 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)

Código fonte:

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

Uma solução de força bruta, cada encadeamento passa por um subconjunto do intervalo (2 ^ 12 itens) e cria um conjunto de valores inteiros que indicam padrões detectados. Estes são então "unidos" juntos e, portanto, a contagem distinta é calculada. Espero que o código não seja muito complicado de seguir, apesar de usar muito as macros de encadeamento . Minhasmain executa o teste algumas vezes para aquecer a JVM.

Atualização: Iterando apenas a metade dos números inteiros, obtém o mesmo resultado devido à simetria. Também pular números com maior contagem de bits na metade inferior do número, pois eles também produzem duplicatas.

Uberjar pré-construído ( 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]

Resultados em diferentes hardwares, o tempo de execução esperado é 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

Você pode facilmente fazer isso com thread único e evitar essa dependência de terceiros usando o padrão para:

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

Bem, o pmap embutido também existe, mas o claypoole tem mais recursos e capacidade de ajuste.

NikoNyrh
fonte
Sim, torna trivial distribuir. Se você tivesse tempo para reavaliar minha solução, tenho certeza de que conseguiria até 30 agora. Não tenho mais otimizações à vista.
NikoNyrh 06/06
Infelizmente, é um não para 30. Tempo decorrido: 217150.87386 msegs
Ah, obrigado por tentar: D Talvez tenha sido melhor ajustar uma curva e interpolar o valor gasto em 120 segundos, mas mesmo assim, esse é um bom desafio.
NikoNyrh
1

Mathematica, n = 19

pressione alt +. abortar e o resultado será impresso

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]
J42161217
fonte
Não posso executar isso, então você poderia explicar como isso evita um tempo exponencial? 2 ^ 241 é um número muito grande!
Você pode mostrar a saída do código?
1
Eu quis dizer f (n) ... fixa
J42161217
1

Mathematica, 21

f [n_]: = Comprimento @
     DeleteDuplicates @
      Transpor@
       Tabela [2> Tr @ IntegerDigits [#, 2] e / @ 
         BitXor [BitShiftRight [#, n - i], Mod [#, 2 ^ i]], {i, 1, 
         n - 1}] & @ Intervalo [0, 2 ^ (n - 1)];
Faça [Imprimir [n -> f @ n], {n, Infinito}]

Para comparação, a resposta de Jenny_mathyn = 19no meu computador.

A parte mais lenta é Tr@IntegerDigits[#, 2] &. É uma pena que o Mathematica não tenha um peso embutido para Hamming.


Se você quiser testar meu código, pode fazer o download de uma avaliação gratuita do Mathematica .

alefalpha
fonte
1

Versão CA, usando popcount embutido

Funciona melhor com clang -O3, mas também funciona se você tiver 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;
}
bartavelle
fonte
Chega a 24 muito rapidamente e depois termina. Você precisa aumentar o limite.
Oh Deus, eu esqueci de remover o código de referência! Vou remover as duas linhas
incorretas
@Lembik deve ser corrigido agora
bartavelle
1

Haskell, (não oficial n = 20)

Essa é apenas a abordagem ingênua - até agora sem otimizações. Eu me perguntava o quão bem isso se sairia contra outras línguas.

Como usá-lo (supondo que você tenha a plataforma haskell instalada):

  • Cole o código em um arquivo approx_corr.hs(ou qualquer outro nome, modifique as etapas a seguir)
  • Navegue para o arquivo e execute ghc approx_corr.hs
  • Corre approx_corr.exe
  • Digite o valor máximo n
  • O resultado de cada cálculo é exibido, bem como o tempo real acumulado (em ms) até esse ponto.

Código:

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
flawr
fonte
O código parece não fornecer saída enquanto é executado. Isso torna um pouco difícil de testar.
Estranho, ele compila sem erro? O que acontece se você tentar compilar o programa main = putStrLn "Hello World!"?
flawr
O Data.Bitsmódulo pode ser útil. Para o loop principal, você pode usar algo como main = do maxn <- getmax; t0 <- gettime; loop 1onde loop n|n==maxn = return ()e loop n = do printresult n (f n); t <- gettime; printtime (t-t0); loop (n+1). getmaxpoderia, por exemplo, getArgsusar os argumentos do programa.
Christian Sievers
@ChristianSievers Muito obrigado !!! Eu fiz essa pergunta no stackoverflow, acho que seria ótimo se você pudesse adicionar isso também!
flawr
Não vejo como responder lá. Você já tem um loop semelhante, e eu não disse nada sobre ganhar tempo: o que você já teve aqui.
Christian Sievers
1

Uma solução Haskell, usando popCount e paralelismo gerenciado manualmente

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

(largue o -llvmse não funcionar)

Corre : ./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..]
bartavelle
fonte
Parece que existe um problema de armazenamento em buffer que eu não recebo nenhuma saída se eu a executar na linha de comando cygwim.
Acabei de atualizar minha solução, mas não sei se isso ajudará muito.
Bartavelle
@Lembik Não tem certeza se isso é óbvio, mas que deve ser compilado com -O3, e pode ser mais rápido com -O3 -fllvm...
bartavelle
(E todos os arquivos de compilação deve ser removido antes de recompilar, se não alterar o código fonte aconteceu)
bartavelle
@Lembik: eu introduzi o paralelismo. Deve ser um pouco mais rápido.
bartavelle
0

Python 2 + pypy, n = 22

Aqui está uma solução Python realmente simples como uma espécie de referência de linha de base.

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)

fonte