Teraz, gdy wszyscy rozwinęli swoją (często zadziwiającą) wiedzę na temat kodowania niskiego poziomu na temat tego, jak bardzo wolno jest Python? (Lub jak szybki jest twój język?) I jak powolny jest Python (część II)? nadszedł czas na wyzwanie, które również zwiększy Twoją zdolność do ulepszania algorytmu.

Poniższy kod oblicza listę długości 9. Pozycja ina liście zlicza, ile razy co najmniej ikolejne zera zostały znalezione podczas obliczania produktów wewnętrznych pomiędzy Fi S. Aby to zrobić dokładnie, iteruje wszystkie możliwe listy Fdługości ni listy Sdługości n+m-1.

#!/usr/bin/python
import itertools
import operator

n=8
m=n+1
leadingzerocounts = [0]*m
for S in itertools.product([-1,1], repeat = n+m-1):
    for F in itertools.product([-1,1], repeat = n):
        i = 0
        while (i<m and sum(map(operator.mul, F, S[i:i+n])) == 0):
            leadingzerocounts[i] +=1
            i+=1
print leadingzerocounts

Dane wyjściowe to

[4587520, 1254400, 347648, 95488, 27264, 9536, 4512, 2128, 1064]

Jeśli zwiększysz n do 10,12,14,16,18,20 za pomocą tego kodu, bardzo szybko staje się on zbyt wolny.

Zasady

• Wyzwanie polega na podaniu prawidłowej wartości wyjściowej dla możliwie największej wartości n. Istotne są tylko parzyste wartości n.
• Jeśli jest remis, wygrana trafia do najszybszego kodu na mojej maszynie dla największej liczby n.
• Zastrzegam sobie prawo do nie testowania kodu, który zajmuje więcej niż 10 minut.
• Możesz zmienić algorytm w dowolny sposób, o ile daje on poprawny wynik. W rzeczywistości będziesz musiał zmienić algorytm, aby zrobić porządny postęp w kierunku wygranej.
• Zwycięzca zostanie nagrodzony tydzień po ustaleniu pytania.
• Nagroda zostanie przyznana w odpowiednim czasie, co jest nieco później niż po przyznaniu zwycięzcy.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod. W związku z tym korzystaj tylko z łatwo dostępnego bezpłatnego oprogramowania i dołącz pełne instrukcje, jak skompilować i uruchomić kod.

Status .

• C . n = 12 w 49 sekund przez @Fors
• Java . n = 16 w 3:07 przez @PeterTaylor
• C ++ . n = 16 w 2:21 przez @ilmale
• Rpython . n = 22 w 3:11 przez @primo
• Java . n = 22 w 6:56 przez @PeterTaylor
• Nimrod . n = 24 w 9:28 sekund przez @ReimerBehrends

Zwycięzcą został Reimer Behrends z wpisem do Nimrod!

Jako sprawdzenie, wyjście dla n = 22 powinno wynosić [12410090985684467712, 2087229562810269696, 351473149499408384, 59178309967151104, 9975110458933248, 1682628717576192, 284866824372224, 48558946385920, 8416739196928, 1518499004416, 301448822784, 71620493312, 22100246528, 8676573184, 3897278464, 1860960256, 911646720, 451520512, 224785920, 112198656, 56062720, 28031360, 14015680].

Konkurs się skończył, ale ... Oferuję 200 punktów za każde zgłoszenie, które wzrasta n o 2 (w ciągu 10 minut na moim komputerze), aż skończą mi się punkty. Ta oferta jest otwarta na zawsze .

Nimrod (N = 22)

import math, locks

const
  N = 20
  M = N + 1
  FSize = (1 shl N)
  FMax = FSize - 1
  SStep = 1 shl (N-1)
  numThreads = 16

type
  ZeroCounter = array[0..M-1, int]
  ComputeThread = TThread[int]

var
  leadingZeros: ZeroCounter
  lock: TLock
  innerProductTable: array[0..FMax, int8]

proc initInnerProductTable =
  for i in 0..FMax:
    innerProductTable[i] = int8(countBits32(int32(i)) - N div 2)

initInnerProductTable()

proc zeroInnerProduct(i: int): bool =
  innerProductTable[i] == 0

proc search2(lz: var ZeroCounter, s, f, i: int) =
  if zeroInnerProduct(s xor f) and i < M:
    lz[i] += 1 shl (M - i - 1)
    search2(lz, (s shr 1) + 0, f, i+1)
    search2(lz, (s shr 1) + SStep, f, i+1)

when defined(gcc):
  const
    unrollDepth = 1
else:
  const
    unrollDepth = 4

template search(lz: var ZeroCounter, s, f, i: int) =
  when i < unrollDepth:
    if zeroInnerProduct(s xor f) and i < M:
      lz[i] += 1 shl (M - i - 1)
      search(lz, (s shr 1) + 0, f, i+1)
      search(lz, (s shr 1) + SStep, f, i+1)
  else:
    search2(lz, s, f, i)

proc worker(base: int) {.thread.} =
  var lz: ZeroCounter
  for f in countup(base, FMax div 2, numThreads):
    for s in 0..FMax:
      search(lz, s, f, 0)
  acquire(lock)
  for i in 0..M-1:
    leadingZeros[i] += lz[i]*2
  release(lock)

proc main =
  var threads: array[numThreads, ComputeThread]
  for i in 0 .. numThreads-1:
    createThread(threads[i], worker, i)
  for i in 0 .. numThreads-1:
    joinThread(threads[i])

initLock(lock)
main()
echo(@leadingZeros)

Połącz z

nimrod cc --threads:on -d:release count.nim

(Nimrod można pobrać tutaj .)

Działa to w wyznaczonym czasie dla n = 20 (i dla n = 18, gdy używa się tylko jednego wątku, w tym drugim przypadku zajmuje to około 2 minut).

Algorytm wykorzystuje wyszukiwanie rekurencyjne, przycinając drzewo wyszukiwania za każdym razem, gdy napotkany zostanie niezerowy produkt wewnętrzny. Skróciliśmy również przestrzeń poszukiwań o połowę, obserwując to dla dowolnej pary wektorów(F, -F) musimy wziąć pod uwagę tylko jeden, ponieważ drugi wytwarza dokładnie takie same zestawy produktów wewnętrznych ( Srównież przez negację ).

Wdrożenie wykorzystuje funkcje metaprogramowania Nimrod do rozwijania / wstawiania pierwszych kilku poziomów wyszukiwania rekurencyjnego. To oszczędza trochę czasu, gdy używasz gcc 4.8 i 4.9 jako backendu Nimroda i sporej kwoty na clang.

Przestrzeń wyszukiwania można dodatkowo przyciąć, obserwując, że musimy wziąć pod uwagę tylko wartości S, które różnią się parzystą liczbą pierwszych N pozycji od naszego wyboru F. Jednak złożoność lub potrzeby pamięciowe tego nie są skalowane dla dużych wartości z N, biorąc pod uwagę, że w takich przypadkach ciało pętli jest całkowicie pomijane.

Tabulowanie, gdzie iloczyn wewnętrzny jest zerowy, wydaje się być szybsze niż używanie jakiejkolwiek funkcji liczenia bitów w pętli. Najwyraźniej dostęp do stołu ma całkiem niezłą lokalizację.

Wygląda na to, że problem powinien być podatny na programowanie dynamiczne, biorąc pod uwagę, jak działa wyszukiwanie rekurencyjne, ale nie ma widocznego sposobu, aby to zrobić przy rozsądnej ilości pamięci.

Przykładowe wyniki:

N = 16:

@[55276229099520, 10855179878400, 2137070108672, 420578918400, 83074121728, 16540581888, 3394347008, 739659776, 183838720, 57447424, 23398912, 10749184, 5223040, 2584896, 1291424, 645200, 322600]

N = 18:

@[3341140958904320, 619683355033600, 115151552380928, 21392898654208, 3982886961152, 744128512000, 141108051968, 27588886528, 5800263680, 1408761856, 438001664, 174358528, 78848000, 38050816, 18762752, 9346816, 4666496, 2333248, 1166624]

N = 20:

@[203141370301382656, 35792910586740736, 6316057966936064, 1114358247587840, 196906665902080, 34848574013440, 6211866460160, 1125329141760, 213330821120, 44175523840, 11014471680, 3520839680, 1431592960, 655872000, 317675520, 156820480, 78077440, 39005440, 19501440, 9750080, 4875040]

Dla celów porównania algorytmu z innymi implementacjami, N = 16 zajmuje około 7,9 sekundy na mojej maszynie przy użyciu jednego wątku i 2,3 sekundy przy użyciu czterech rdzeni.

N = 22 zajmuje około 15 minut na 64-rdzeniowej maszynie z gcc 4.4.6, ponieważ backend Nimroda przepełnia 64-bitowe liczby całkowite leadingZeros[0](być może nie podpisane, nie oglądałem).

Aktualizacja: Znalazłem miejsce na kilka ulepszeń. Po pierwsze, dla danej wartości Fmożemy dokładnie wyliczyć pierwsze 16 pozycji odpowiednich Swektorów, ponieważ muszą się one różnić dokładnie w różnych N/2miejscach. Więc wstępnego wyliczenia listy bitowych wektorów wielkości N, które mają N/2bitów zestawu i korzystają z nich czerpać początkową część Sz F.

Po drugie, możemy ulepszyć wyszukiwanie rekurencyjne, obserwując, że zawsze znamy wartość F[N](ponieważ MSB ma zero w reprezentacji bitowej). To pozwala nam dokładnie przewidzieć, w której gałęzi wracamy z produktu wewnętrznego. Chociaż faktycznie pozwoliłoby nam to przekształcić całe wyszukiwanie w pętlę rekurencyjną, tak naprawdę zdarza się, że dość spieszy przewidywanie gałęzi, więc utrzymujemy najwyższe poziomy w oryginalnej formie. Wciąż oszczędzamy trochę czasu, przede wszystkim poprzez zmniejszenie liczby rozgałęzień, które wykonujemy.

W przypadku niektórych porządków kod używa teraz liczb całkowitych bez znaku i naprawia je w wersji 64-bitowej (na wypadek, gdyby ktoś chciał uruchomić to w architekturze 32-bitowej).

Ogólne przyspieszenie wynosi od współczynnika x3 do x4. N = 22 nadal potrzebuje więcej niż ośmiu rdzeni, aby uruchomić się w czasie krótszym niż 10 minut, ale na 64-rdzeniowej maszynie jest to teraz około czterech minut (znumThreads odpowiednio zwiększone). Nie sądzę jednak, aby było dużo więcej miejsca na ulepszenia bez innego algorytmu.

N = 22:

@[12410090985684467712, 2087229562810269696, 351473149499408384, 59178309967151104, 9975110458933248, 1682628717576192, 284866824372224, 48558946385920, 8416739196928, 1518499004416, 301448822784, 71620493312, 22100246528, 8676573184, 3897278464, 1860960256, 911646720, 451520512, 224785920, 112198656, 56062720, 28031360, 14015680]

Zaktualizowano ponownie, wykorzystując dalsze możliwe ograniczenia przestrzeni wyszukiwania. Działa za około 9:49 minut dla N = 22 na mojej maszynie quadcore.

Ostatnia aktualizacja (myślę). Lepsze klasy równoważności dla wyborów F, skrócenie czasu działania dla N = 22 do 3:19 minut 57 sekund (edycja: przypadkowo uruchomiłem to z tylko jednym wątkiem) na mojej maszynie.

Ta zmiana wykorzystuje fakt, że para wektorów wytwarza te same zera wiodące, jeśli jeden można przekształcić w drugi, obracając go. Niestety, dość krytyczna optymalizacja niskiego poziomu wymaga, aby górny bit F w reprezentacji bitów był zawsze taki sam, a podczas korzystania z tej równoważności nieco zmniejszył przestrzeń wyszukiwania i skrócił czas działania o około jedną czwartą w porównaniu z użyciem innej przestrzeni stanów redukcja F, narzut związany z eliminacją optymalizacji niskiego poziomu bardziej niż ją zrównoważył. Okazuje się jednak, że problem ten można wyeliminować również biorąc pod uwagę fakt, że F, które są odwrotne względem siebie, są również równoważne. Wprawdzie nieco to skomplikowało obliczenia klas równoważności, ale pozwoliło mi również zachować wspomnianą optymalizację niskiego poziomu, co doprowadziło do przyspieszenia około x3.

Jeszcze jedna aktualizacja do obsługi 128-bitowych liczb całkowitych dla zgromadzonych danych. Aby skompilować z 128 bitowych liczb całkowitych, trzeba longint.nimz tutaj i skompilować -d:use128bit. N = 24 nadal zajmuje więcej niż 10 minut, ale dla zainteresowanych zainteresowałem się poniższym wynikiem.

N = 24:

@[761152247121980686336, 122682715414070296576, 19793870419291799552, 3193295704340561920, 515628872377565184, 83289931274780672, 13484616786640896, 2191103969198080, 359662314586112, 60521536552960, 10893677035520, 2293940617216, 631498735616, 230983794688, 102068682752, 48748969984, 23993655296, 11932487680, 5955725312, 2975736832, 1487591936, 743737600, 371864192, 185931328, 92965664]

import math, locks, unsigned

when defined(use128bit):
  import longint
else:
  type int128 = uint64 # Fallback on unsupported architectures
  template toInt128(x: expr): expr = uint64(x)

const
  N = 22
  M = N + 1
  FSize = (1 shl N)
  FMax = FSize - 1
  SStep = 1 shl (N-1)
  numThreads = 16

type
  ZeroCounter = array[0..M-1, uint64]
  ZeroCounterLong = array[0..M-1, int128]
  ComputeThread = TThread[int]
  Pair = tuple[value, weight: int32]

var
  leadingZeros: ZeroCounterLong
  lock: TLock
  innerProductTable: array[0..FMax, int8]
  zeroInnerProductList = newSeq[int32]()
  equiv: array[0..FMax, int32]
  fTable = newSeq[Pair]()

proc initInnerProductTables =
  for i in 0..FMax:
    innerProductTable[i] = int8(countBits32(int32(i)) - N div 2)
    if innerProductTable[i] == 0:
      if (i and 1) == 0:
        add(zeroInnerProductList, int32(i))

initInnerProductTables()

proc ror1(x: int): int {.inline.} =
  ((x shr 1) or (x shl (N-1))) and FMax

proc initEquivClasses =
  add(fTable, (0'i32, 1'i32))
  for i in 1..FMax:
    var r = i
    var found = false
    block loop:
      for j in 0..N-1:
        for m in [0, FMax]:
          if equiv[r xor m] != 0:
            fTable[equiv[r xor m]-1].weight += 1
            found = true
            break loop
        r = ror1(r)
    if not found:
      equiv[i] = int32(len(fTable)+1)
      add(fTable, (int32(i), 1'i32))

initEquivClasses()

when defined(gcc):
  const unrollDepth = 4
else:
  const unrollDepth = 4

proc search2(lz: var ZeroCounter, s0, f, w: int) =
  var s = s0
  for i in unrollDepth..M-1:
    lz[i] = lz[i] + uint64(w)
    s = s shr 1
    case innerProductTable[s xor f]
    of 0:
      # s = s + 0
    of -1:
      s = s + SStep
    else:
      return

template search(lz: var ZeroCounter, s, f, w, i: int) =
  when i < unrollDepth:
    lz[i] = lz[i] + uint64(w)
    if i < M-1:
      let s2 = s shr 1
      case innerProductTable[s2 xor f]
      of 0:
        search(lz, s2 + 0, f, w, i+1)
      of -1:
        search(lz, s2 + SStep, f, w, i+1)
      else:
        discard
  else:
    search2(lz, s, f, w)

proc worker(base: int) {.thread.} =
  var lz: ZeroCounter
  for fi in countup(base, len(fTable)-1, numThreads):
    let (fp, w) = fTable[fi]
    let f = if (fp and (FSize div 2)) == 0: fp else: fp xor FMax
    for sp in zeroInnerProductList:
      let s = f xor sp
      search(lz, s, f, w, 0)
  acquire(lock)
  for i in 0..M-1:
    let t = lz[i].toInt128 shl (M-i).toInt128
    leadingZeros[i] = leadingZeros[i] + t
  release(lock)

proc main =
  var threads: array[numThreads, ComputeThread]
  for i in 0 .. numThreads-1:
    createThread(threads[i], worker, i)
  for i in 0 .. numThreads-1:
    joinThread(threads[i])

initLock(lock)
main()
echo(@leadingZeros)

Java (n=22 ?)

Myślę, że większość odpowiedzi jest lepsza niż n=16 podobne podejście, chociaż różnią się symetriami, które wykorzystują i sposobem podziału zadania między wątkami.

Wektory zdefiniowane w pytaniu można zastąpić ciągami bitów, a iloczyn wewnętrzny za pomocą XOR nakładającego się na siebie okna i sprawdzania, czy dokładnie n/2ustawione są bity (a tym samym n/2bity usuwane). Istnieją n! / ((n/2)!)(środkowy współczynnik dwumianowy) ciągi nbitów z n/2zestawem bitów (które nazywam ciągami zrównoważonymi ), więc dla każdego z nich Fistnieje tak wiele okien, Sktóre dają zerowy iloczyn wewnętrzny. Ponadto czynność przesuwania się Swzdłuż jednego i sprawdzania, czy nadal możemy znaleźć przychodzący bit, który daje zerowy iloczyn wewnętrzny, odpowiada szukaniu krawędzi na wykresie, którego węzłami są okna, a których krawędzie łączą węzeł uz węzłem, vktórego pierwsze bity z .n-1 bity są ostatnien-1u

Na przykład za pomocą n=6i F=001001otrzymujemy ten wykres:

i dla F=001011otrzymujemy ten wykres:

Następnie musimy liczyć na siebie iod 0do nilu ścieżki o długości iistnieją zsumowanie na wykresach dla każdego F. Myślę, że większość z nas korzysta z wyszukiwania w głąb.

Zauważ, że wykresy są rzadkie: łatwo jest udowodnić, że każdy węzeł ma stopień co najwyżej 1, a stopień co najwyżej jeden. Oznacza to również, że jedynymi możliwymi strukturami są proste łańcuchy i proste pętle. To trochę upraszcza DFS.

Wykorzystuję kilka symetrii: zbalansowane łańcuchy są zamknięte pod odwrotnym bitem ( ~operacja w wielu językach z rodziny ALGOL) i pod rotacją bitów, dzięki czemu możemy grupować razem wartości, Fktóre są powiązane tymi operacjami i wykonują tylko DFS pewnego razu.

public class CodeGolf26459v8D implements Runnable {
    private static final int NUM_THREADS = 8;

    public static void main(String[] args) {
        v8D(22);
    }

    private static void v8D(int n) {
        int[] bk = new int[1 << n];
        int off = 0;
        for (int i = 0; i < bk.length; i++) {
            bk[i] = Integer.bitCount(i) == n/2 ? off++ : -1;
        }

        int[] fwd = new int[off];
        for (int i = 0; i < bk.length; i++) {
            if (bk[i] >= 0) fwd[bk[i]] = i;
        }

        CodeGolf26459v8D[] runners = new CodeGolf26459v8D[NUM_THREADS];
        Thread[] threads = new Thread[runners.length];
        for (int i = 0; i < runners.length; i++) {
            runners[i] = new CodeGolf26459v8D(n, i, runners.length, bk, fwd);
            threads[i] = new Thread(runners[i]);
            threads[i].start();
        }

        try {
            for (int i = 0; i < threads.length; i++) threads[i].join();
        }
        catch (InterruptedException ie) {
            throw new RuntimeException("This shouldn't be reachable", ie);
        }

        long surviving = ((long)fwd.length) << (n - 1);
        for (int i = 0; i <= n; i++) {
            for (CodeGolf26459v8D runner : runners) surviving -= runner.survival[i];
            System.out.print(i == 0 ? "[" : ", ");
            java.math.BigInteger result = new java.math.BigInteger(Long.toString(surviving));
            System.out.print(result.shiftLeft(n + 1 - i));
        }
        System.out.println("]");
    }

    public final int n;
    protected final int id;
    protected final int numRunners;
    private final int[] bk;
    private final int[] fwd;

    public long[] survival;

    public CodeGolf26459v8D(int n, int id, int numRunners, int[] bk, int[] fwd) {
        this.n = n;
        this.id = id;
        this.numRunners = numRunners;

        this.bk = bk;
        this.fwd = fwd;
    }

    private int dfs2(int[] graphShape, int flip, int i) {
        if (graphShape[i] != 0) return graphShape[i];

        int succ = flip ^ (fwd[i] << 1);
        if (succ >= bk.length) succ ^= bk.length + 1;

        int j = bk[succ];
        if (j == -1) return graphShape[i] = 1;

        graphShape[i] = n + 1; // To detect cycles
        return graphShape[i] = dfs2(graphShape, flip, j) + 1;
    }

    @Override
    public void run() {
        int n = this.n;
        int[] bk = this.bk;
        int[] fwd = this.fwd;

        // NB The initial count is approx 2^(2n - 1.33 - 0.5 lg n)
        // For n=18 we overflow 32-bit
        // 64-bit is good up to n=32.
        long[] survival = new long[n + 1];
        boolean[] visited = new boolean[1 << (n - 1)];
        int th = 0;
        for (int f = 0; f < visited.length; f++) {
            if (visited[f]) continue;

            int m = 1, g = f;
            while (true) {
                visited[g] = true;
                int ng = g << 1;
                if ((ng >> (n - 1)) != 0) ng ^= (1 << n) - 1;
                if (ng == f) break;
                m++;
                g = ng;
            }

            if (th++ % numRunners != id) continue;

            int[] graphShape = new int[fwd.length];
            int flip = (f << 1) ^ f;
            for (int i = 0; i < graphShape.length; i++) {
                int life = dfs2(graphShape, flip, i);
                if (life <= n) survival[life] += m;
            }
        }

        this.survival = survival;
    }
}

Na moim rdzeniu 2 2,5 GHz dostaję

# n=18
$ javac CodeGolf26459v8D.java && time java CodeGolf26459v8D
[3341140958904320, 619683355033600, 115151552380928, 21392898654208, 3982886961152, 744128512000, 141108051968, 27588886528, 5800263680, 1408761856, 438001664, 174358528, 78848000, 38050816, 18762752, 9346816, 4666496, 2333248, 1166624]

real    0m3.131s
user    0m10.133s
sys     0m0.380s

# n=20
$ javac CodeGolf26459v8D.java && time java CodeGolf26459v8D
[203141370301382656, 35792910586740736, 6316057966936064, 1114358247587840, 196906665902080, 34848574013440, 6211866460160, 1125329141760, 213330821120, 44175523840, 11014471680, 3520839680, 1431592960, 655872000, 317675520, 156820480, 78077440, 39005440, 19501440, 9750080, 4875040]

real    1m8.706s
user    4m20.980s
sys     0m0.564s

# n=22
$ javac CodeGolf26459v8D.java && time java CodeGolf26459v8D
[12410090985684467712, 2087229562810269696, 351473149499408384, 59178309967151104, 9975110458933248, 1682628717576192, 284866824372224, 48558946385920, 8416739196928, 1518499004416, 301448822784, 71620493312, 22100246528, 8676573184, 3897278464, 1860960256, 911646720, 451520512, 224785920, 112198656, 56062720, 28031360, 14015680]

real    20m10.654s
user    76m53.880s
sys     0m6.852s

Ponieważ komputer Lembika ma 8 rdzeni i wykonał mój poprzedni program jednowątkowy dwa razy szybciej niż mój, jestem optymistą, że uruchomi się on n=22w mniej niż 8 minut.

do

Jest to w zasadzie tylko nieco zoptymalizowana implementacja algorytmu w pytaniu. Może zarządzać n=12w terminie.

#include <stdio.h>
#include <inttypes.h>

#define n 12
#define m (n + 1)

int main() {
    int i;
    uint64_t S, F, o[m] = {0};
    for (S = 0; S < (1LLU << (n + m - 1)); S += 2)
        for (F = 0; F < (1 << (n - 1)); F++)
            for (i = 0; i < m; i++)
                if (__builtin_popcount(((S >> i) & ((1 << n) - 1)) ^ F) == n >> 1)
                    o[i] += 4;
                else
                    break;
    for (i = 0; i < m; i++)
        printf("%" PRIu64 " ", o[i]);
    return 0;
}

Uruchomienie testowe n=12, w tym kompilacja:

$ clang -O3 -march=native -fstrict-aliasing -ftree-vectorize -Wall fast.c
$ time ./a.out 
15502147584 3497066496 792854528 179535872 41181184 9826304 2603008 883712 381952 177920 85504 42560 21280 
real    0m53.266s
user    0m53.042s
sys     0m0.068s
$

Komentarz: Właśnie włączyłem swój mózg i użyłem prostych kombinacji kombinatoryki, aby obliczyć, że zawsze pierwsza wartość będzie n! / ((n / 2)!)^2 * 2^(n + m - 1). Wydaje mi się, że musi istnieć całkowicie algebraiczne rozwiązanie tego problemu.

Jawa, n=16

Dla dowolnej wartości Fistnieją \binom{n}{n/2}wektory, które mają zerowy iloczyn wewnętrzny. Możemy więc zbudować wykres, którego wierzchołki to pasujące wektory i których krawędzie odpowiadają przesunięciu S, a następnie wystarczy policzyć ścieżki długości do nwykresu.

Nie próbowałem mikrooptymalizować tego poprzez zamianę operacji warunkowych na operacje bitowe, ale każda podwójna inkrementacja nzwiększa czas działania około 16-krotnie, więc nie zrobi to wystarczającej różnicy, chyba że jestem dość blisko progu. Na moim komputerze nie jestem.

public class CodeGolf26459 {

    public static void main(String[] args) {
        v3(16);
    }

    // Order of 2^(2n-1) * n ops
    private static void v3(int n) {
        long[] counts = new long[n+1];
        int mask = (1 << n) - 1;
        for (int f = 0; f < (1 << (n-1)); f++) {
            // Find adjacencies
            long[] subcounts = new long[1 << n];
            for (int g = 0; g < (1 << n); g++) {
                subcounts[g] = Integer.bitCount(f ^ g) == n/2 ? 2 : -1;
            }

            for (int round = 0; round <= n; round++) {
                long count = 0;
                // Extend one bit.
                long[] next = new long[1 << n];
                for (int i = 0; i < (1 << n); i++) {
                    long s = subcounts[i];
                    if (s == -1) next[i] = -1;
                    else {
                        count += s;
                        int j = (i << 1) & mask;
                        if (subcounts[j] >= 0) next[j] += s;
                        if (subcounts[j + 1] >= 0) next[j + 1] += s;
                    }
                }
                counts[round] += count << (n - round);
                subcounts = next;
            }
        }

        System.out.print("[");
        for (long count : counts) System.out.print(count+", ");
        System.out.println("]");
    }
}

Na moim rdzeniu 2 2,5 GHz dostaję

$ javac CodeGolf26459.java && time java -server CodeGolf26459 
[55276229099520, 10855179878400, 2137070108672, 420578918400, 83074121728, 16540581888, 3394347008, 739659776, 183838720, 57447424, 23398912, 10749184, 5223040, 2584896, 1291424, 645200, 322600, ]

real    6m2.663s
user    6m4.631s
sys     0m1.580s

RPython, N = 22 ~ 3: 23

Wielowątkowy, wykorzystujący rekurencyjne zejście bez stosu. Program akceptuje dwa argumenty wiersza poleceń: N i liczbę wątków roboczych.

from time import sleep

from rpython.rlib.rthread import start_new_thread, allocate_lock
from rpython.rlib.rarithmetic import r_int64, build_int, widen
from rpython.rlib.rbigint import rbigint

r_int8 = build_int('r_char', True, 8)

class ThreadEnv:
  __slots__ = ['n', 'counts', 'num_threads',
               'v_range', 'v_num', 'running', 'lock']

  def __init__(self):
    self.n = 0
    self.counts = [rbigint.fromint(0)]
    self.num_threads = 0
    self.v_range = [0]
    self.v_num = 0
    self.running = 0
    self.lock = None

env = ThreadEnv()

bt_bits = 12
bt_mask = (1<<bt_bits)-1
# computed compile time
bit_table = [r_int8(0)]
for i in xrange(1,1<<bt_bits):
  bit_table += [((i&1)<<1) + bit_table[i>>1]]

def main(argv):
  argc = len(argv)
  if argc < 2 or argc > 3:
    print 'Usage: %s N [NUM_THREADS=2]'%argv[0]
    return 1

  if argc == 3:
    env.num_threads = int(argv[2])
  else:
    env.num_threads = 2

  env.n = int(argv[1])
  env.counts = [rbigint.fromint(0)]*env.n
  env.lock = allocate_lock()

  v_range = []
  v_max = 1<<(env.n-1)
  v_num = 0
  v = (1<<(env.n>>1))-1
  while v < v_max:
    v_num += 1
    v_range += [v]
    if v&1:
      # special case odd v
      s = (v+1)&-v
      v ^= s|(s>>1)
    else:
      s = v&-v
      r = v+s
      # s is at least 2, skip two iterations
      i = 3
      s >>= 2
      while s:
        i += 1
        s >>= 1
      v = r|((v^r)>>i)
  env.v_range = v_range
  env.v_num = v_num

  for i in xrange(env.num_threads-1):
    start_new_thread(run,())

  # use the main process as a worker
  run()

  # wait for any laggers
  while env.running:
    sleep(0.05)

  result = []
  for i in range(env.n):
    result += [env.counts[i].lshift(env.n-i+3).str()]
  result += [env.counts[env.n-1].lshift(3).str()]
  print result
  return 0

def run():
  with env.lock:
    v_start = env.running
    env.running += 1

  n = env.n
  counts = [r_int64(0)]*n
  mask = (1<<n)-1
  v_range = env.v_range
  v_num = env.v_num
  z_count = 1<<(n-2)

  for i in xrange(v_start, v_num, env.num_threads):
    v = v_range[i]
    counts[0] += z_count
    counts[1] += v_num
    r = v^(v<<1)
    for w in v_range:
      # unroll counts[2] for speed
      # ideally, we could loop over x directly,
      # rather than over all v, only to throw the majority away
      # there's a 2x-3x speed improvement to be had here...
      x = w^r
      if widen(bit_table[x>>bt_bits]) + widen(bit_table[x&bt_mask]) == n:
        counts[2] += 1
        x, y = v, x
        o, k = 2, 3
        while k < n:
          # x = F ^ S
          # y = F ^ (S<<1)
          o = k
          z = (((x^y)<<1)^y)&mask
          # z is now F ^ (S<<2), possibly xor 1
          # what S and F actually are is of no consequence

          # the compiler hint `widen` let's the translator know
          # to store the result as a native int, rather than a signed char
          bt_high = widen(bit_table[z>>bt_bits])
          if bt_high + widen(bit_table[z&bt_mask]) == n:
            counts[k] += 1
            x, y = y, z
            k += 1
          elif bt_high + widen(bit_table[(z^1)&bt_mask]) == n:
            counts[k] += 1
            x, y = y, z^1
            k += 1
          else: k = n

  with env.lock:
    for i in xrange(n):
      env.counts[i] = env.counts[i].add(rbigint.fromrarith_int(counts[i]))
    env.running -= 1

def target(*args):
  return main, None

Kompilować

Stwórz lokalny klon repozytorium PyPy przy użyciu mercurial, git lub cokolwiek innego. Wpisz następującą inkantację (zakładając, że powyższy skrypt ma nazwę convolution-high.py):

$ pypy %PYPY_REPO%/rpython/bin/rpython --thread convolution-high.py

gdzie %PYPY_REPO%reprezentuje zmienną środowiskową wskazującą na właśnie sklonowane repozytorium. Kompilacja zajmuje około minuty.

Przykładowe czasy

N = 16, 4 wątki:

$ timeit convolution-high-c 16 4
[55276229099520, 10855179878400, 2137070108672, 420578918400, 83074121728, 16540581888, 3394347008, 739659776, 183838720, 57447424, 23398912, 10749184, 5223040, 2584896, 1291424, 645200, 322600]
Elapsed Time:     0:00:00.109
Process Time:     0:00:00.390

N = 18, 4 wątki:

$ timeit convolution-high-c 18 4
[3341140958904320, 619683355033600, 115151552380928, 21392898654208, 3982886961152, 744128512000, 141108051968, 27588886528, 5800263680, 1408761856, 438001664, 174358528, 78848000, 38050816, 18762752, 9346816, 4666496, 2333248, 1166624]
Elapsed Time:     0:00:01.250
Process Time:     0:00:04.937

N = 20, 4 wątki:

$ timeit convolution-high-c 20 4
[203141370301382656, 35792910586740736, 6316057966936064, 1114358247587840, 196906665902080, 34848574013440, 6211866460160, 1125329141760, 213330821120, 44175523840, 11014471680, 3520839680, 1431592960, 655872000, 317675520, 156820480, 78077440, 39005440, 19501440, 9750080, 4875040]
Elapsed Time:     0:00:15.531
Process Time:     0:01:01.328

N = 22, 4 wątki:

$ timeit convolution-high-c 22 4
[12410090985684467712, 2087229562810269696, 351473149499408384, 59178309967151104, 9975110458933248, 1682628717576192, 284866824372224, 48558946385920, 8416739196928, 1518499004416, 301448822784, 71620493312, 22100246528, 8676573184, 3897278464, 1860960256, 911646720, 451520512, 224785920, 112198656, 56062720, 28031360, 14015680]
Elapsed Time:     0:03:23.156
Process Time:     0:13:25.437

Python 3.3, N = 20, 3,5 min

Oświadczenie: moim zamiarem NIE jest publikowanie tego jako mojej własnej odpowiedzi, ponieważ algorytm, którego używam, jest tylko bezwstydnym portem z rozwiązania RPython primo . Moim celem tutaj jest pokazanie, co możesz zrobić w Pythonie, jeśli połączysz magię Numpy i Numba modułów .

Numba wyjaśnił w skrócie:

Numba jest specjalizującym się w czasie kompilatorem, który kompiluje adnotacje w języku Python i NumPy do LLVM (przez dekoratory). http://numba.pydata.org/

Aktualizacja 1 : Zauważyłem po rzuceniu liczbami, że możemy po prostu całkowicie pominąć niektóre z nich. Więc teraz maxf staje się (1 << n) // 2, a maxs staje się maxf 2 **. Przyspieszy to nieco proces. n = 16 zajmuje teraz tylko ~ 48 s (w porównaniu z 4,5 min). Mam też inny pomysł, który spróbuję sprawdzić, czy uda mi się przyspieszyć.

Aktualizacja 2: Zmieniony algorytm (rozwiązanie primo). Chociaż mój port nie obsługuje jeszcze wielowątkowości, dodanie go jest dość proste. Możliwe jest nawet wydanie CPython GIL przy użyciu Numba i ctypów. To rozwiązanie działa jednak bardzo szybko również na jednym rdzeniu!

import numpy as np
import numba as nb

bt_bits = 11
bt_mask = (1 << bt_bits) - 1
bit_table = np.zeros(1 << bt_bits, np.int32)

for i in range(0, 1 << bt_bits):
    bit_table[i] = ((i & 1) << 1) + bit_table[i >> 1]

@nb.njit("void(int32, int32, int32, int32, int64[:], int64[:])")
def run(n, m, start, re, counts, result):
    mask = (1 << n) - 1

    v_max = (1 << n) // 2
    rr = v_max // 2

    v = (1 << (n >> 1)) - 1
    while v < v_max:
        s = start

        while s < rr:
            f = v ^ s
            counts[0] += 8
            t = s << 1
            o, j = 0, 1

            while o < j and j < m:
                o = j
                w = (t ^ f) & mask
                bt_high = bit_table[w >> bt_bits]

                if bt_high + bit_table[w & bt_mask] == n:
                    counts[j] += 8
                    t <<= 1
                    j += 1
                elif bt_high + bit_table[(w ^ 1) & bt_mask] == n:
                    counts[j] += 8
                    t = (t | 1) << 1
                    j += 1
                    s += re

            s = v & -v
            r = v + s
            o = v ^ r
            o = (o >> 2) // s
            v = r | o

    for e in range(m):
        result[e] += counts[e] << (n - e)

I w końcu:

if __name__ == "__main__":
    n = 20
    m = n + 1

    result = np.zeros(m, np.int64)
    counts = np.zeros(m, np.int64)

    s1 = time.time() * 1000
    run(n, m, 0, 1, counts, result)
    s2 = time.time() * 1000

    print(result)
    print("{0}ms".format(s2 - s1))

Działa to na moim komputerze w 212688ms lub ~ 3.5min.

C++ N=16

I'm testing on a EEEPC with an atom.. my time don't make a lot of sense. :D
The atom solve n=14 in 34 seconds. And n=16 in 20 minutes. I want to test n = 16 on OP pc. I'm optimistic.

The ideas is that every time we find a solution for a given F we have found 2^i solution because we can change the lower part of S leading to the same result.

#include <stdio.h>
#include <cinttypes>
#include <cstring>

int main()
{
   const int n = 16;
   const int m = n + 1;
   const uint64_t maxS = 1ULL << (2*n);
   const uint64_t maxF = 1ULL << n;
   const uint64_t mask = (1ULL << n)-1;
   uint64_t out[m]={0};
   uint64_t temp[m] = {0};
   for( uint64_t F = 0; F < maxF; ++F )
   {
      for( uint64_t S = 0; S < maxS; ++S )
      {
         int numSolution = 1;
         for( int i = n; i >= 0; --i )
         {
            const uint64_t window = S >> i;
            if( __builtin_popcount( mask & (window ^ F) ) == (n / 2) )
            {
               temp[i] += 1;
            } else {
               numSolution = 1 << i;
               S += numSolution - 1;
               break;
            }
         }
         for( int i = n; i >= 0; --i )
         {
            out[i] += temp[i]*numSolution;
            temp[i] = 0;
         }
      }
   }
   for( int i = n; i >= 0; --i )
   {
      uint64_t x = out[i];
      printf( "%lu ", x );
   }
   return 0;
}

to compile:

gcc 26459.cpp -std=c++11 -O3 -march=native -fstrict-aliasing -ftree-vectorize -Wall -pedantic -o 26459

JAVASCRIPT n:12

In my computer it took 231.242 seconds. In the Demo I'm using webworkers to prevent freezing the browser. This sure can be further improved with parallel workers. I know JS don't stand a chance in this challenge but I did it for fun!

Click to run the online Demo

var n = 8;        
var m = n + 1;
var o = [];
var popCount = function(bits) {
  var SK5  = 0x55555555,
      SK3  = 0x33333333,
      SKF0 = 0x0f0f0f0f,
      SKFF = 0xff00ff;

  bits -= (bits >> 1) & SK5;
  bits  = (bits & SK3) + ((bits >> 2) & SK3);
  bits  = (bits & SKF0) + ((bits >> 4) & SKF0);
  bits += bits >> 8;

  return (bits + (bits >> 15)) & 63;
};
for(var S = 0; S < (1 << n + m - 1); S += 2){
  for(var F = 0; F < (1 << n - 1); F += 1){
    for (var i = 0; i < m; i++){
      var c = popCount(((S >> i) & ((1 << n) - 1)) ^ F);
      if(c == n >> 1){
        if(!o[i]) o[i] = 0;
        o[i] += 4;
      } else break;
    }
  }
}
return o;

Fortran: n=12

I just made a quick'n'dirty version in Fortran, no optimizations except OpenMP. It should squeeze in at just below 10 minutes for n=12 on the OPs machine, it takes 10:39 on my machine which is sligthly slower.

64-bit integers have a negative performance impact indeed; guess I would have to rethink the whole algorithm for this to be much faster. Don't know if I'm going to bother, I think I'll rather spend some time thinking up a good challenge that's more to my own tastes. If anyone else wants to take this and run with it, go ahead :)

program golf
use iso_fortran_env
implicit none
integer, parameter ::  n=12
integer :: F(n), S(2*n)
integer(int64) :: leadingzerocounts(n+1)
integer :: k
integer(int64) :: i,j,bindec,enc

leadingzerocounts=0

!$OMP parallel do private(i,enc,j,bindec,S,F,k) reduction(+:leadingzerocounts) schedule(dynamic)
do i=0,2**(2*n)-1
  enc=i
  ! Short loop to convert i into the array S with -1s and 1s
  do j=2*n,1,-1
    bindec=2**(j-1)
    if (enc-bindec .ge. 0) then
      S(j)=1
      enc=enc-bindec
    else
      S(j)=-1
    endif
  end do
  do j=0,2**(n)-1
    ! Convert j into the array F with -1s and 1s
    enc=j
    do k=n,1,-1
      bindec=2**(k-1)
      if (enc-bindec .ge. 0) then
        F(k)=1
        enc=enc-bindec
      else
        F(k)=-1
      endif
    end do
    ! Compute dot product   
    do k=1,n+1
      if (dot_product(F,S(k:k+n-1)) /= 0) exit
      leadingzerocounts(k)=leadingzerocounts(k)+1
    end do
  end do
end do
!$OMP end parallel do

print *, leadingzerocounts

end

Lua: n = 16

Disclaimer: my intention is NOT to post this as my own answer, since the algorithm I'm using is shamelessly stolen from Anna Jokela's clever answer. which was shamelessly stolen from ilmale's clever answer.

Besides, it's not even valid - it has inaccuracies caused by floating point numbers (it would be better if Lua would support 64-bit integers). However, I'm still uploading it, just to show how fast this solution is. It's a dynamic programming language, and yet I can calculate n = 16 in reasonable time (1 minute on 800MHz CPU).

Run with LuaJIT, standard interpreter is too slow.

local bit = require "bit"
local band = bit.band
local bor = bit.bor
local bxor = bit.bxor
local lshift = bit.lshift
local rshift = bit.rshift

-- http://stackoverflow.com/a/11283689/736054
local function pop_count(w)
    local b1 = 1431655765
    local b2 = 858993459
    local b3 = 252645135
    local b7 = 63

    w = band(rshift(w, 1), b1) + band(w, b1)
    w = band(rshift(w, 2), b2) + band(w, b2)
    w = band(w + rshift(w, 4), b3)
    return band(rshift(w, 24) + rshift(w, 16) + rshift(w, 8) + w, b7)
end

local function gen_array(n, value)
    value = value or 0
    array = {}
    for i = 1, n do
        array[i] = value
    end
    return array
end

local n = 16
local u = math.floor(n / 2)
local m = n + 1
local maxf = math.floor(lshift(1, n) / 2)
local maxs = maxf ^ 2
local mask = lshift(1, n) - 1

local out = gen_array(m, 0)
local temp = gen_array(m, 0)


for f = 0, maxf do
    local s = 0
    while s <= maxs do
        local num_solution = 1

        for i = n, 0, -1 do
            if pop_count(band(mask, bxor(rshift(s, i), f))) == u then
                temp[i + 1] = temp[i + 1] + 8
            else
                num_solution = lshift(1, i)
                s = s + num_solution - 1
                break
            end
        end

        for i = 1, m do
            out[i] = out[i] + temp[i] * num_solution
            temp[i] = 0
        end

        s = s + 1
    end
end

for i = m, 1, -1 do
    print(out[i])
end