Ahora que todos han desarrollado su (a menudo sorprendente) experiencia de codificación de bajo nivel para ¿Cuán lento es realmente Python? (¿O qué tan rápido es su idioma?) Y ¿Cuán lento es realmente Python (Parte II)? Es hora de un desafío que también ampliará su capacidad para mejorar un algoritmo.

El siguiente código calcula una lista de longitud 9. La posición ien la lista cuenta el número de veces que se encontraron al menos iceros consecutivos al calcular los productos internos entre Fy S. Para hacer esto exactamente, itera sobre todas las listas posibles Fde longitud ny listas Sde longitud 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

La salida es

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

Si aumenta n a 10,12,14,16,18,20 usando este código, muy rápidamente se vuelve demasiado lento.

Reglas

• El desafío es dar la salida correcta para el mayor número posible. Solo los valores pares de n son relevantes.
• Si hay un empate, la victoria va al código más rápido en mi máquina para el n más grande.
• Me reservo el derecho de no probar el código que lleva más de 10 minutos.
• Puede cambiar el algoritmo de la forma que desee, siempre que proporcione la salida correcta. De hecho , tendrá que cambiar el algoritmo para hacer un progreso decente hacia ganar.
• El ganador recibirá una semana desde que se estableció la pregunta.
• La recompensa se otorgará cuando venza, un poco después de que se otorgue el ganador.

Mi máquina Los tiempos se ejecutarán en mi máquina. Esta es una instalación estándar de ubuntu en un procesador AMD FX-8350 de ocho núcleos. Esto también significa que necesito poder ejecutar su código. Como consecuencia, solo use software gratuito fácilmente disponible e incluya instrucciones completas sobre cómo compilar y ejecutar su código.

Estado .

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

¡El ganador fue Reimer Behrends con una entrada en Nimrod!

Como verificación, la salida para n = 22 debería ser [12410090985684467712, 2087229562810269696, 351473149499408384, 59178309967151104, 9975110458933248, 1682628717576192, 284866824372224, 48558946385920, 8416739196928, 1518499004416, 301448822784, 71620493312, 22100246528, 8676573184, 3897278464, 1860960256, 911646720, 451520512, 224785920, 112198656, 56062720, 28031360, 14015680].

La competencia ya terminó pero ... ofreceré 200 puntos por cada envío que aumente n en 2 (dentro de los 10 minutos en mi computadora) hasta que me queden sin puntos. Esta oferta está abierta para siempre .

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)

Compilar con

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

(Nimrod se puede descargar aquí ).

Esto se ejecuta en el tiempo asignado para n = 20 (y para n = 18 cuando solo se usa un solo subproceso, lo que demora aproximadamente 2 minutos en el último caso).

El algoritmo utiliza una búsqueda recursiva, podando el árbol de búsqueda cada vez que se encuentra un producto interno distinto de cero. También reducimos el espacio de búsqueda a la mitad al observar que para cualquier par de vectores (F, -F)solo necesitamos considerar uno porque el otro produce exactamente los mismos conjuntos de productos internos (al negar Stambién).

La implementación utiliza las funciones de metaprogramación de Nimrod para desenrollar / en línea los primeros niveles de la búsqueda recursiva. Esto ahorra un poco de tiempo al usar gcc 4.8 y 4.9 como el backend de Nimrod y una buena cantidad para el sonido metálico.

El espacio de búsqueda podría reducirse aún más al observar que solo necesitamos considerar valores de S que difieran en un número par de las primeras N posiciones de nuestra elección de F. Sin embargo, la complejidad o las necesidades de memoria de eso no escalan para valores grandes de N, dado que el cuerpo del bucle se omite por completo en esos casos.

Tabular donde el producto interno es cero parece ser más rápido que usar cualquier funcionalidad de conteo de bits en el ciclo. Aparentemente acceder a la mesa tiene bastante buena localidad.

Parece que el problema debería ser susceptible a la programación dinámica, considerando cómo funciona la búsqueda recursiva, pero no hay una forma aparente de hacerlo con una cantidad razonable de memoria.

Salidas de ejemplo:

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]

Para comparar el algoritmo con otras implementaciones, N = 16 toma aproximadamente 7.9 segundos en mi máquina cuando se usa un solo hilo y 2.3 segundos cuando se usan cuatro núcleos.

N = 22 tarda unos 15 minutos en una máquina de 64 núcleos con gcc 4.4.6 como backend de Nimrod y desborda enteros de 64 bits leadingZeros[0](posiblemente no sin firmar, no lo he mirado).

Actualización: he encontrado espacio para un par de mejoras más. Primero, para un valor dado de F, podemos enumerar las primeras 16 entradas de los Svectores correspondientes con precisión, porque deben diferir en N/2lugares exactos . Así precomputamos una lista de vectores de bits de tamaño Nque tienen N/2determinados bits y utilizar estos para derivar la parte inicial de Spartir F.

En segundo lugar, podemos mejorar la búsqueda recursiva al observar que siempre conocemos el valor de F[N](ya que el MSB es cero en la representación de bits). Esto nos permite predecir con precisión a qué rama recurrimos desde el producto interno. Si bien eso realmente nos permitiría convertir toda la búsqueda en un bucle recursivo, eso realmente arruina bastante la predicción de rama, por lo que mantenemos los niveles superiores en su forma original. Todavía ahorramos algo de tiempo, principalmente al reducir la cantidad de ramificaciones que estamos haciendo.

Para alguna limpieza, el código ahora usa enteros sin signo y los repara a 64 bits (en caso de que alguien quiera ejecutar esto en una arquitectura de 32 bits).

La aceleración general está entre un factor de x3 y x4. N = 22 todavía necesita más de ocho núcleos para ejecutarse en menos de 10 minutos, pero en una máquina de 64 núcleos ahora se reduce a unos cuatro minutos (con un numThreadsaumento correspondiente en consecuencia). Sin embargo, no creo que haya mucho más margen de mejora sin un algoritmo diferente.

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]

Actualizada nuevamente, haciendo uso de otras posibles reducciones en el espacio de búsqueda. Se ejecuta en aproximadamente 9:49 minutos para N = 22 en mi máquina quadcore.

Actualización final (creo). Mejores clases de equivalencia para las opciones de F, reduciendo el tiempo de ejecución de N = 22 a 3:19 minutos y 57 segundos (editar: accidentalmente lo ejecuté con solo un hilo) en mi máquina.

Este cambio hace uso del hecho de que un par de vectores produce los mismos ceros iniciales si uno puede transformarse en el otro al rotarlo. Desafortunadamente, una optimización de bajo nivel bastante crítica requiere que el bit superior de F en la representación de bits sea siempre el mismo, y mientras usa esta equivalencia recorta un poco el espacio de búsqueda y reduce el tiempo de ejecución en aproximadamente un cuarto sobre el uso de un espacio de estado diferente reducción en F, la sobrecarga de eliminar la optimización de bajo nivel más que compensarla. Sin embargo, resulta que este problema puede eliminarse considerando también el hecho de que F que son inversas entre sí también son equivalentes. Si bien esto aumentó un poco la complejidad del cálculo de las clases de equivalencia, también me permitió retener la optimización de bajo nivel mencionada anteriormente, lo que condujo a una aceleración de aproximadamente x3.

Una actualización más para admitir enteros de 128 bits para los datos acumulados. Para compilar con enteros de 128 bits, necesitará longint.nimdesde aquí y compilar con -d:use128bit. N = 24 todavía lleva más de 10 minutos, pero he incluido el resultado a continuación para los interesados.

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

Creo que la mayoría de las respuestas son mejores que n=16usar un enfoque similar a este, aunque difieren en las simetrías que explotan y la forma en que dividen la tarea entre hilos.

Los vectores definidos en la pregunta se pueden reemplazar con cadenas de bits, y el producto interno con XOR haciendo coincidir la ventana superpuesta y verificando que haya exactamente n/2bits establecidos (y, por lo tanto, n/2bits borrados). Hay n! / ((n/2)!)cadenas de nbits (coeficiente binomial central) con n/2bits establecidos (que yo llamo cadenas equilibradas ), por lo que para cualquier dato Fhay tantas ventanas Sque dan un producto interno cero. Además, la acción de deslizarse a lo Slargo de uno y verificar si todavía podemos encontrar un bit entrante que proporcione un producto interno cero corresponde a buscar un borde en un gráfico cuyos nodos son las ventanas y cuyos bordes unen un nodo ua un nodo vcuyos primeros n-1bits son los ultimosn-1bits de u.

Por ejemplo, con n=6y F=001001obtenemos este gráfico:

y para F=001011obtener este gráfico:

Luego, necesitamos contar para cada uno ide 0a ncuántos caminos de longitud ihay, sumando los gráficos para cada uno F. Creo que la mayoría de nosotros estamos usando la búsqueda en profundidad.

Tenga en cuenta que los gráficos son escasos: es fácil demostrar que cada nodo tiene un grado máximo de 1 y un grado máximo de uno. Eso también significa que las únicas estructuras posibles son cadenas simples y bucles simples. Esto simplifica un poco el DFS.

Exploto un par de simetrías: las cadenas equilibradas están cerradas bajo bit inverso (la ~operación en muchos idiomas de la familia ALGOL) y bajo rotación de bits, por lo que podemos agrupar valores de los Fque están relacionados por estas operaciones y solo hacemos el DFS una vez.

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

En mi Core 2 de 2.5GHz obtengo

# 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

Como la computadora de Lembik tiene 8 núcleos y ejecutó mi programa anterior de un solo subproceso dos veces más rápido que el mío, soy optimista de que se ejecutará n=22en menos de 8 minutos.

do

Básicamente es solo una implementación ligeramente optimizada del algoritmo en la pregunta. Se puede administrar n=12dentro del límite de tiempo.

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

Prueba de funcionamiento n=12, incluida la compilación:

$ 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
$

Comentario: acabo de encender mi cerebro y usé algunos combinatorios simples para calcular que el primer valor siempre será n! / ((n / 2)!)^2 * 2^(n + m - 1). Me parece que debe haber una solución completamente algebraica para este problema.

Java, n=16

Para cualquier valor dado de Fhay \binom{n}{n/2}vectores que tienen un producto interno cero con él. Por lo tanto, podemos construir un gráfico cuyos vértices sean aquellos vectores coincidentes y cuyos bordes correspondan al desplazamiento de S, y luego solo necesitamos contar las rutas de longitud hasta nen el gráfico.

No he intentado microoptimizar esto reemplazando los condicionales con operaciones bit a bit, pero cada incremento doble naumenta el tiempo de ejecución aproximadamente 16 veces, por lo que eso no hará una diferencia suficiente a menos que esté bastante cerca del umbral. En mi máquina, no lo soy.

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("]");
    }
}

En mi Core 2 de 2.5GHz obtengo

$ 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

Multihilo, utilizando un descenso recursivo sin pila. El programa acepta dos argumentos de línea de comando: N y el número de subprocesos de trabajo.

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

Compilar

Haga un clon local del repositorio PyPy usando mercurial, git o lo que prefiera. Escriba el siguiente encantamiento (suponiendo que se nombre el script anterior convolution-high.py):

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

donde %PYPY_REPO%representa una variable de entorno que apunta al repositorio que acaba de clonar. La compilación lleva aproximadamente un minuto.

Tiempos de muestra

N = 16, 4 hilos:

$ 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 hilos:

$ 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 hilos:

$ 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 hilos:

$ 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.5min

Descargo de responsabilidad: mi intención NO es publicar esto como mi propia respuesta, ya que el algoritmo que estoy usando es solo un puerto descarado de la solución RPython de primo . Mi propósito aquí es sólo para mostrar lo que puede hacer en Python si se combinan la magia de numpy y Numba módulos.

Numba explicó en resumen:

Numba es un compilador especializado justo a tiempo que compila código anotado de Python y NumPy en LLVM (a través de decoradores). http://numba.pydata.org/

Actualización 1 : noté después de arrojar los números que simplemente podemos omitir algunos de los números por completo. Entonces ahora maxf se convierte en (1 << n) // 2 y maxs se convierte en maxf 2 **. Esto acelerará el proceso bastante. n = 16 ahora solo toma ~ 48s (antes 4,5 min). También tengo otra idea que intentaré y veré si puedo hacerlo un poco más rápido.

Actualización 2: Algoritmo modificado (solución de primo). Si bien mi puerto aún no admite subprocesos múltiples, es bastante trivial agregarlo. Incluso es posible liberar CPython GIL usando Numba y ctypes. Sin embargo, esta solución también funciona muy rápido en un solo núcleo.

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)

Y finalmente:

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

Esto se ejecuta en mi máquina en 212688ms o ~ 3.5min.

C ++ N = 16

Estoy probando en un EEEPC con un átomo ... mi tiempo no tiene mucho sentido. : D
El átomo resuelve n = 14 en 34 segundos. Yn = 16 en 20 minutos. Quiero probar n = 16 en PC OP. Soy optimista.

La idea es que cada vez que encontremos una solución para una F dada, hemos encontrado una solución 2 ^ i porque podemos cambiar la parte inferior de S que conduce al mismo resultado.

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

compilar:

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

JAVASCRIPT n: 12

En mi computadora tardó 231.242 segundos. En la demostración estoy usando trabajadores web para evitar congelar el navegador. Esto seguramente se puede mejorar aún más con trabajadores paralelos. Sé que JS no tiene ninguna oportunidad en este desafío, ¡pero lo hice por diversión!

Haga clic para ejecutar la demostración en línea

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

Acabo de hacer una versión rápida y sucia en Fortran, sin optimizaciones excepto OpenMP. Debería exprimirse por debajo de los 10 minutos durante n = 12 en la máquina de OP, toma 10:39 en mi máquina, que es un poco más lenta.

Los enteros de 64 bits tienen un impacto negativo en el rendimiento; Supongo que tendría que repensar todo el algoritmo para que esto sea mucho más rápido. No sé si me voy a molestar, creo que preferiré pasar un tiempo pensando en un buen desafío que sea más para mis gustos. Si alguien más quiere tomar esto y ejecutarlo, adelante :)

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

Descargo de responsabilidad: mi intención NO es publicar esto como mi propia respuesta, ya que el algoritmo que estoy usando es descaradamente robado de la respuesta inteligente de Anna Jokela . que fue robado descaradamente de la inteligente respuesta de ilmale .

Además, ni siquiera es válido: tiene imprecisiones causadas por números de coma flotante (sería mejor si Lua admitiera enteros de 64 bits). Sin embargo, todavía lo estoy cargando, solo para mostrar qué tan rápida es esta solución. Es un lenguaje de programación dinámico y, sin embargo, puedo calcular n = 16 en un tiempo razonable (1 minuto en una CPU de 800MHz).

Ejecutar con LuaJIT, el intérprete estándar es demasiado lento.

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