Menemukan perkiraan korelasi

14

Pertimbangkan Spanjang string biner n. Dari pengindeksan 1, kita dapat menghitung jarak Hamming antara S[1..i+1]dan S[n-i..n]untuk semua idalam urutan dari 0ke n-1. Jarak Hamming antara dua string dengan panjang yang sama adalah jumlah posisi di mana simbol yang sesuai berbeda. Sebagai contoh,

S = 01010

memberi

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

Ini karena 0pertandingan 0, 01memiliki jarak Hamming dua 10, 010pertandingan 010, 0101memiliki jarak Hamming empat 1010 dan akhirnya 01010cocok dengan dirinya sendiri.

Kami hanya tertarik pada output di mana jarak Hamming paling banyak 1, namun. Jadi dalam tugas ini kami akan melaporkan Yjika jarak Hamming paling banyak satu dan Nsebaliknya. Jadi dalam contoh kita di atas kita akan dapatkan

[Y, N, Y, N, Y]

Tentukan f(n)jumlah array yang berbeda dari Ys dan Ns yang didapat ketika iterasi semua 2^nstring bit Spanjang yang mungkin berbeda n.

Tugas

Untuk meningkatkan nmulai dari 1, kode Anda harus ditampilkan f(n).

Contoh jawaban

Sebab n = 1..24, jawaban yang benar adalah:

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

Mencetak gol

Kode Anda harus beralih dari n = 1memberikan jawaban untuk masing-masing npada gilirannya. Saya akan mengatur waktu seluruh pelarian, membunuhnya setelah dua menit.

Skor Anda adalah yang tertinggi yang nAnda dapatkan pada waktu itu.

Dalam kasus seri, jawaban pertama menang.

Di mana kode saya akan diuji?

Saya akan menjalankan kode Anda di laptop Windows 7 (agak lama) saya di bawah cygwin. Sebagai hasilnya, tolong berikan bantuan yang Anda bisa untuk mempermudah ini.

Laptop saya memiliki RAM 8GB dan CPU Intel i7 [email protected] GHz (Broadwell) dengan 2 core dan 4 thread. Set instruksi termasuk SSE4.2, AVX, AVX2, FMA3 dan TSX.

Entri terkemuka per bahasa

  • n = 40 di Rust menggunakan CryptoMiniSat, oleh Anders Kaseorg. (Dalam VM tamu Lubuntu di bawah Vbox.)
  • n = 35 dalam C ++ menggunakan perpustakaan BuDDy, oleh Christian Seviers. (Dalam VM tamu Lubuntu di bawah Vbox.)
  • n = 34 di Clingo oleh Anders Kaseorg. (Dalam VM tamu Lubuntu di bawah Vbox.)
  • n = 31 in Rust oleh Anders Kaseorg.
  • n = 29 in Clojure oleh NikoNyrh.
  • n = 29 in C oleh bartavelle.
  • n = 27 in Haskell oleh bartavelle
  • n = 24 in Pari / gp oleh alephalpha.
  • n = 22 in Python 2 + pypy oleh saya.
  • n = 21 dalam Mathematica by alephalpha. (Dilaporkan sendiri)

Karunia masa depan

Sekarang saya akan memberikan hadiah 200 poin untuk setiap jawaban yang mencapai n = 80 pada mesin saya dalam dua menit.


sumber
Apakah Anda tahu beberapa trik yang akan memungkinkan seseorang untuk menemukan algoritma yang lebih cepat daripada brute force yang naif? Jika tidak, tantangan ini adalah "tolong terapkan ini di x86" (atau mungkin jika kami tahu GPU Anda ...).
Jonathan Allan
@ JonathanAllan Sangat mungkin untuk mempercepat pendekatan yang sangat naif. Seberapa cepat Anda bisa mendapatkan saya tidak yakin. Menariknya, jika kita mengubah pertanyaan sehingga Anda mendapatkan Y jika jarak Hamming paling banyak 0 dan N sebaliknya, maka ada rumus form tertutup yang dikenal.
@Lembik Apakah kita mengukur waktu CPU atau waktu nyata?
flawr
@ flawr Saya mengukur waktu nyata tetapi menjalankannya beberapa kali dan mengambil minimum untuk menghilangkan keanehan.

Jawaban:

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"

Bagaimana itu bekerja

Ini melakukan pencarian rekursif melalui pohon dari semua tugas parsial untuk awalan array Y / N, menggunakan SAT solver untuk memeriksa pada setiap langkah apakah tugas parsial saat ini konsisten dan memangkas pencarian jika tidak. CryptoMiniSat adalah pemecah SAT yang tepat untuk pekerjaan ini karena optimasi khusus untuk klausa XOR.

Tiga keluarga kendala adalah

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

kecuali itu, sebagai optimisasi, S 0 dipaksa ke false, sehingga D k 0 sama dengan S k .

Anders Kaseorg
sumber
2
Woohoooooo! :)
Saya masih mencoba untuk mengkompilasi ini di Windows (menggunakan cygwin + gcc). Saya mengkloning cryptominisat dan mengkompilasinya. Tapi saya masih tidak tahu bagaimana cara mengkompilasi kode karat. Ketika saya melakukannya cargo buildsaya mendapatkan--- stderr CMake Error: Could not create named generator Visual Studio 14 2015 Win64
2
@ rahnema1 Terima kasih, tapi sepertinya masalahnya ada pada sistem build CMake dari pustaka C ++ yang disematkan di dalam peti cryptominisat, bukan dengan Rust sendiri.
Anders Kaseorg
1
@Lembik saya mendapatkan 404 dari pasta itu.
Mego
1
@ChristianSievers Pertanyaan bagus. Itu bekerja tetapi tampaknya sedikit lebih lambat (2 × atau lebih). Saya tidak yakin mengapa itu tidak sebaik, jadi mungkin CryptoMiniSat belum dioptimalkan dengan baik untuk beban kerja tambahan semacam itu.
Anders Kaseorg
9

Karat, n ≈ 30 atau 31 atau 32

Di laptop saya (dua core, i5-6200U), ini bisa melalui n = 1,…, 31 dalam 53 detik, menggunakan sekitar 2,5 GiB memori, atau melalui n = 1,…, 32 dalam 105 detik, menggunakan sekitar 5 GiB memori. Kompilasi dengan cargo build --releasedan jalankan 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"

Cobalah online!

Saya juga memiliki varian yang sedikit lebih lambat menggunakan memori yang sangat sedikit.

Anders Kaseorg
sumber
Optimalisasi apa yang Anda gunakan?
1
@Lembik Optimasi terbesar, selain melakukan semuanya dengan aritmatika bitwise dalam bahasa yang dikompilasi, adalah menggunakan hanya nondeterminisme sebanyak yang diperlukan untuk memakukan awalan array Y / N. Saya melakukan pencarian rekursif pada kemungkinan awalan dari array Y / N, membawa serta vektor string yang mungkin mencapai awalan itu, tetapi hanya string yang bagian tengahnya yang tidak teruji diisi dengan nol. Yang mengatakan, ini masih pencarian eksponensial, dan optimasi ini hanya mempercepatnya dengan faktor polinomial.
Anders Kaseorg
Itu jawaban yang bagus. Terima kasih. Saya berharap seseorang akan menggali ke dalam kombinatorik untuk mendapatkan kecepatan yang signifikan.
@Lembik Saya telah memperbaiki bug pemborosan memori, melakukan lebih banyak optimasi mikro, dan menambahkan paralelisme. Harap uji ulang ketika Anda mendapat kesempatan — saya berharap dapat meningkatkan skor saya sebesar 1 atau 2. Apakah Anda memiliki ide kombinatorial dalam pikiran untuk speedup yang lebih besar? Saya belum datang dengan apa pun.
Anders Kaseorg
1
@Lembik Tidak ada rumus yang diberikan pada entri OEIS. (Kode Mathematica di sana juga tampaknya menggunakan brute-force.) Jika Anda mengetahuinya, Anda mungkin ingin memberi tahu mereka tentang hal itu.
Christian Sievers
6

C ++ menggunakan perpustakaan BuDDy

Pendekatan yang berbeda: memiliki rumus biner (sebagai diagram keputusan biner ) yang mengambil bit Ssebagai input dan jika benar yang memberikan beberapa nilai tetap Yatau Npada posisi tertentu yang dipilih. Jika rumus itu tidak konstan salah, pilih posisi bebas dan berulang, coba keduanya Ydan N. Jika tidak ada posisi bebas, kami telah menemukan kemungkinan nilai output. Jika rumusnya konstan salah, mundurlah.

Ini bekerja relatif masuk akal karena ada beberapa nilai yang mungkin sehingga kita dapat sering mundur lebih awal. Saya mencoba ide yang sama dengan pemecah SAT, tetapi itu kurang berhasil.

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

Untuk mengkompilasi dengan debian 8 (jessie), instal libbdd-devdan lakukan g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. Mungkin bermanfaat untuk menambah argumen pertama menjadi bdd_initlebih banyak.

Sievers Kristen
sumber
Ini terlihat menarik. Apa yang Anda dapatkan dengan ini?
@Lembik Saya mendapatkan 31 dalam 100-an pada perangkat keras yang sangat lama yang tidak akan membiarkan saya menjawab lebih cepat
Christian Sievers
Bantuan apa pun yang dapat Anda berikan tentang cara mengompilasi ini di Windows (misalnya menggunakan cygwin) dengan penuh terima kasih.
@Lembik Saya tidak tahu tentang Windws tetapi github.com/fd00/yacp/tree/master/buddy tampaknya membantu wrt cygwin
Christian Sievers
1
Wow, oke, Anda meyakinkan saya bahwa saya perlu menambahkan perpustakaan ini ke toolkit saya. Sudah selesai dilakukan dengan baik!
Anders Kaseorg
4

Clingo, n ≈ 30 atau 31 34

Saya sedikit terkejut melihat lima baris kode Clingo mengambil alih solusi Rust brute-force saya dan mendekati solusi BuDDy Christian — sepertinya itu akan mengalahkan itu juga dengan batas waktu yang lebih tinggi.

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

merencanakan

Anders Kaseorg
sumber
Ini bagus! Terlihat dari grafik Anda bahwa solusi BuDDy tiba-tiba menjadi lebih buruk. Ada yang tahu kenapa?
@Lembik Saya belum menyelidiki BuDDy cukup untuk memastikan, tapi mungkin kehabisan cache pada saat itu?
Anders Kaseorg
Wow! Saya pikir nilai pertama yang lebih tinggi untuk bdd_initmembantu, atau memungkinkan untuk meningkatkan tabel simpul lebih banyak dengan menelepon bdd_setmaxincreasedengan nilai jauh di atas standar 50000. - Apakah Anda menggunakan versi program saya yang berubah?
Christian Sievers
2
Saya suka grafik Anda.
1
Anda mendapatkan peningkatan kinerja yang mengejutkan dengan menggunakan opsi --configuration=crafty( jumpydan trendymemberikan hasil yang serupa).
Christian Sievers
2

Pari / GP , 23

Secara default, Pari / GP membatasi ukuran stack hingga 8 MB. Baris pertama kode default(parisize, "4g"),, menetapkan batas ini menjadi 4 GB. Jika masih memberikan aliran stackover, Anda dapat mengaturnya hingga 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)))
alephalpha
sumber
Mencapai 22 dan kemudian memberikan stackoverflow.
Mendapat 24 sekarang.
2

Clojure, 29 dalam 75 38 detik, 30 dalam 80 dan 31 dalam 165

Runtimes dari Intel i7 6700K , penggunaan memori kurang dari 200 MB.

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

Kode sumber:

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

Solusi brute-force, setiap thread melewati subset rentang (2 ^ 12 item) dan membangun satu set nilai integer yang menunjukkan pola yang terdeteksi. Ini kemudian "disatukan" bersama-sama dan dengan demikian penghitungan yang berbeda dihitung. Saya harap kode ini tidak terlalu sulit untuk diikuti walaupun menggunakan banyak threading macro . Sayamain menjalankan tes beberapa kali untuk mendapatkan pemanasan JVM.

Pembaruan: Iterasi lebih dari setengah bilangan bulat, mendapat hasil yang sama karena simetri. Juga melewatkan angka dengan jumlah bit lebih tinggi pada bagian bawah dari angka saat mereka menghasilkan duplikat juga.

Uberjar pra-dibangun ( 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]

Hasil pada perangkat keras yang berbeda, runtime yang diharapkan adalah 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

Anda dapat dengan mudah membuat ini single-threaded dan menghindari ketergantungan pihak ke-3 itu dengan menggunakan standar untuk:

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

Nah, built-in pmap juga ada tetapi claypoole memiliki lebih banyak fitur dan kemampuan merapikan.

NikoNyrh
sumber
Ya, itu membuatnya sepele untuk didistribusikan. Apakah Anda punya waktu untuk mengevaluasi kembali solusi saya, saya yakin Anda akan mendapatkannya hingga 30 sekarang. Saya tidak memiliki pengoptimalan lebih lanjut yang terlihat.
NikoNyrh
Sayangnya itu tidak untuk 30. Waktu berlalu: 217150.87386 msecs
Ahaa, terima kasih telah mencobanya: D Mungkin lebih baik untuk menyesuaikan kurva ini dan menginterpolasi bahwa di mana nilai desimal 120 detik dihabiskan, tetapi meskipun ini adalah tantangan yang bagus.
NikoNyrh
1

Mathematica, n = 19

tekan alt +. untuk membatalkan dan hasilnya akan dicetak

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
sumber
Saya tidak dapat menjalankan ini, jadi bisakah Anda menjelaskan bagaimana cara menghindari mengambil waktu yang eksponensial? 2 ^ 241 adalah angka yang sangat besar!
Bisakah Anda menunjukkan output dari kode?
1
Maksud saya f (n) ... fix
J42161217
1

Mathematica, 21

f [n_]: = Panjangnya @
     DeleteDuplicates @
      Mengubah urutan@
       Tabel [2> Tr @ IntegerDigits [#, 2] & / @ 
         BitXor [BitShiftRight [#, n - i], Mod [#, 2 ^ i]], {i, 1, 
         n - 1}] & @ Rentang [0, 2 ^ (n - 1)];
Lakukan [Cetak [n -> f @ n], {n, Infinity}]

Sebagai perbandingan, jawaban Jenny_mathy diberikan n = 19di komputer saya.

Bagian paling lambat adalah Tr@IntegerDigits[#, 2] &. Sangat disayangkan bahwa Mathematica tidak memiliki built-in untuk berat Hamming.


Jika Anda ingin menguji kode saya, Anda dapat mengunduh uji coba gratis Mathematica .

alephalpha
sumber
1

Versi AC, menggunakan popcount bawaan

Bekerja dengan lebih baik clang -O3, tetapi juga berfungsi jika Anda hanya memiliki 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
sumber
Menjadi 24 sangat cepat dan kemudian berakhir. Anda perlu menambah batas.
Ya Tuhan, saya lupa menghapus kode benchmark! Saya akan menghapus dua baris yang menyinggung: /
bartavelle
@Lembik harus diperbaiki sekarang
bartavelle
1

Haskell, (tidak resmi n = 20)

Ini hanyalah pendekatan naif - sejauh ini tanpa optimasi. Saya bertanya-tanya seberapa bagus tarifnya terhadap bahasa lain.

Cara menggunakannya (dengan asumsi Anda telah menginstal platform haskell ):

  • Tempel kode dalam satu file approx_corr.hs(atau nama lain, ubah langkah-langkah berikut sesuai)
  • Arahkan ke file dan jalankan ghc approx_corr.hs
  • Lari approx_corr.exe
  • Masukkan maksimal n
  • Hasil dari setiap perhitungan ditampilkan, serta kumulatif waktu nyata (dalam ms) hingga saat itu.

Kode:

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
cacat
sumber
Tampaknya kode tidak memberikan output saat berjalan. Ini membuatnya agak sulit untuk diuji.
Aneh, apakah kompilasi tanpa kesalahan? Apa yang terjadi jika Anda mencoba mengompilasi program main = putStrLn "Hello World!"?
flawr
The Data.BitsModul mungkin berguna. Untuk loop utama Anda, Anda bisa menggunakan sesuatu seperti di main = do maxn <- getmax; t0 <- gettime; loop 1mana loop n|n==maxn = return ()dan loop n = do printresult n (f n); t <- gettime; printtime (t-t0); loop (n+1). getmaxmisalnya bisa digunakan getArgsuntuk menggunakan argumen program.
Christian Sievers
@ChristianSievers Terima kasih banyak !!! Saya mengajukan pertanyaan ini di stackoverflow, saya pikir akan lebih bagus jika Anda bisa menambahkannya di sana juga!
flawr
Saya tidak melihat bagaimana menjawab di sana. Anda memiliki loop serupa di sana, dan saya tidak mengatakan apa-apa tentang mendapatkan waktu: yang sudah Anda miliki di sini.
Christian Sievers
1

Solusi Haskell, menggunakan popCount dan paralelisme yang dikelola secara manual

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

(jatuhkan -llvm jika tidak bekerja)

Lari : ./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
sumber
Ada masalah buffering sepertinya saya tidak mendapatkan output sama sekali jika saya menjalankannya dari baris perintah cygwim.
Saya baru saja memperbarui solusi saya, tetapi saya tidak tahu apakah itu akan banyak membantu.
bartavelle
@Lembik Tidak yakin apakah itu jelas, tetapi itu harus dikompilasi -O3, dan mungkin lebih cepat dengan -O3 -fllvm...
bartavelle
(Dan semua file build harus dihapus sebelum dikompilasi ulang, jika tidak terjadi perubahan kode sumber)
bartavelle
@Lembik: Saya memperkenalkan paralelisme. Seharusnya sedikit lebih cepat.
bartavelle
0

Python 2 + pypy, n = 22

Berikut ini adalah solusi Python yang sangat sederhana sebagai semacam tolok ukur dasar.

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)

sumber