Este es un seguimiento de ¿Qué tan lento es realmente Python? (¿O qué tan rápido es tu idioma?) .

Resulta que fue un poco demasiado fácil obtener una aceleración x100 para mi última pregunta. Para aquellos que han disfrutado el desafío pero quieren algo más difícil donde realmente pueden usar sus habilidades de bajo nivel, aquí está la parte II. El desafío es obtener una aceleración x100 para el siguiente código de Python según lo probado en mi computadora.

Para hacerlo más difícil, estoy usando pypy esta vez. El tiempo actual para mí es de 1 minuto y 7 segundos usando pypy 2.2.1.

Reglas

1. La primera persona en enviar el código que puedo ejecutar, es correcta y es x100 veces más rápido en mi máquina, recibirá una recompensa de 50 puntos.
2. Otorgaré la victoria al código más rápido después de una semana.

import itertools 
import operator 
import random

n = 8 
m  = 8 
iters = 1000  

# creates an array of 0s with length m
# [0, 0, 0, 0, 0, 0, 0, 0]
leadingzerocounts = [0]*m

# itertools.product creates an array of all possible combinations of the 
# args passed to it.
#
# Ex:
#   itertools.product("ABCD", "xy") --> Ax Ay Bx By Cx Cy Dx Dy
#   itertools.product("AB", repeat=5) --> [
#    ('A', 'A', 'A', 'A', 'A'),
#    ('A', 'A', 'A', 'A', 'B'),
#    ('A', 'A', 'A', 'B', 'A'),
#    ('A', 'A', 'A', 'B', 'B'),
#    etc.
#   ]
for S in itertools.product([-1,1], repeat = n+m-1):
    for i in xrange(iters):
        F = [random.choice([-1,0,0,1]) for j in xrange(n)]

        # if the array is made up of only zeros keep recreating it until
        # there is at least one nonzero value.
        while not any(F):
            F = [random.choice([-1,0,0,1]) for j in xrange(n)]

        j = 0
        while (j < m and sum(map(operator.mul, F, S[j:j+n])) == 0):
            leadingzerocounts[j] +=1
            j += 1
print leadingzerocounts

La salida debe ser similar a

[6335185, 2526840, 1041967, 439735, 193391, 87083, 40635, 19694]

Debe usar una semilla aleatoria en su código y se aceptará cualquier generador de números aleatorios que sea lo suficientemente bueno como para dar respuestas cercanas a las anteriores.

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.

Explicación de código

Este código itera sobre todas las matrices S de longitud n + m-1 que están formadas por -1s y 1s. Para cada matriz S, muestra 1000 matrices aleatorias distintas de cero F de longitud n formadas por -1,0 o 1 con probabilidad de 1/4, 1/2, / 14 de tomar cada valor. Luego calcula los productos internos entre F y cada ventana de S de longitud n hasta que encuentra un producto interno distinto de cero. Agrega 1 a leadingzerocountscada posición donde encontró un producto interno cero.

Estado

• Perl . 2.7 veces la ralentización por @tobyink. (En comparación con pypy no cpython).

• J . 39 veces más velocidad por @Eelvex.

• C . 59 veces acelerado por @ace.
• Julia . 197 veces más rápido sin incluir el tiempo de inicio en un minuto más. 8.5 veces más rápido, incluido el tiempo de inicio (en este caso, es más rápido usar 4 procesadores que 8).
• Fortran . 438 veces más rápido por @ semi-extrínseco.
• Rpython . 258 veces acelerado por @primo.
• C ++ . 508 veces acelerado por @ilmale.

(Dejé de cronometrar las nuevas mejoras porque son demasiado rápidas e iters era demasiado pequeño).

Se señaló que los tiempos por debajo de un segundo no son confiables y también algunos idiomas tienen un costo inicial. El argumento es que si va a incluir eso, también debe incluir el tiempo de compilación de C / C ++, etc. Aquí están los tiempos para el código más rápido con el número de iteraciones aumentado a 100,000.

• Julia . 42 segundos por @ un minuto más.
• C ++ . 14 segundos por @GuySirton.
• Fortran . 14s por @ semi-extrinsic.
• C ++ . 12s por @ilmale.
• Rpython . 18s por @primo.
• C ++ . 5s por @Stefan.

El ganador es ... ¡Stefan!

Reto de seguimiento publicado. ¿Qué tan alto puedes llegar? (Un desafío de codificación + algoritmos) . Este es más difícil.

C ++ bit magic

~ 16 ms multiproceso, 56 ms de subproceso único. ~ 4000 aceleración.

(la aceleración se basa en un código multiproceso en mi i7-2820QM y los 1 min 9 segundos mencionados en la pregunta. Dado que el sistema OP tiene peor rendimiento de un solo subproceso que mi CPU pero mejor rendimiento de subprocesos múltiples, espero que este número sea preciso)

La parte multiproceso es bastante ineficiente debido al desove de los hilos. Probablemente podría hacerlo mejor aprovechando mi biblioteca de trabajos personalizada, pero esa tiene errores en sistemas Unix. Para obtener una explicación y un código casi idéntico sin subprocesos, consulte https://codegolf.stackexchange.com/a/26485/20965 .

editar

Le di a cada hilo su propio RNG y reduje la longitud de bits a 32, lo que redujo el tiempo de ejecución en unos pocos ms.

#include <iostream>
#include <bitset>
#include <random>
#include <chrono>
#include <stdint.h>
#include <cassert>
#include <array>
#include <tuple>
#include <memory>
#include <thread>
#include <future>
#include <string.h>


#ifdef _MSC_VER
uint32_t popcnt( uint32_t x ){ return _mm_popcnt_u32(x); }
#else
uint32_t popcnt( uint32_t x ){ return __builtin_popcount(x); }
#endif



void convolve()
{
    static const unsigned threadCount = 32;
    static const unsigned n = 8;
    static const unsigned m = 8;
    static const unsigned totalIters = 1000;
    static_assert( n <= 16, "packing of F fails when n > 16.");
    static uint32_t fmask = (1 << n) -1; fmask |= fmask << 16;

    std::array< uint32_t, m * threadCount > out;
    std::vector< std::future<void> > threads;

    for( int threadId = 0; threadId < threadCount; threadId++)
    {
        threads.emplace_back( std::async( [&, threadId]
        {
            std::random_device rd;
            std::knuth_b gen(rd());
            uint32_t nextRandomNumber = gen();

            const unsigned iters = totalIters / threadCount;

            std::array< uint32_t, m > leadingZeros;
            for( auto& x : leadingZeros )
                x = 0;

            for( unsigned i = 0; i < iters; i++ )
            {
                // generate random bit mess
                uint32_t F;
                do {
                    // this funky looking construction shortens the dependancy chain of F
                    F = nextRandomNumber & fmask;
                    nextRandomNumber = gen();
                } while ( 0 == ((F % (1 << n)) ^ (F >> 16 )) );

                // Assume F is an array with interleaved elements such that F[0] || F[16] is one element
                // here MSB(F) & ~LSB(F) returns 1 for all elements that are positive
                // and  ~MSB(F) & LSB(F) returns 1 for all elements that are negative
                // this results in the distribution ( -1, 0, 0, 1 )
                // to ease calculations we generate r = LSB(F) and l = MSB(F)

                uint32_t r = F % ( 1 << n );
                // modulo is required because the behaviour of the leftmost bit is implementation defined
                uint32_t l = ( F >> 16 ) % ( 1 << n );

                uint32_t posBits = l & ~r;
                uint32_t negBits = ~l & r;
                assert( (posBits & negBits) == 0 );

                uint32_t mask = posBits | negBits;
                uint32_t totalBits = popcnt( mask );
                // if the amount of -1 and +1's is uneven, sum(S*F) cannot possibly evaluate to 0
                if ( totalBits & 1 )
                    continue;

                uint32_t adjF = posBits & ~negBits;
                uint32_t desiredBits = totalBits / 2;

                uint32_t S = (1 << (n + m -1));
                // generate all possible N+1 bit strings
                // 1 = +1
                // 0 = -1
                while ( S-- )
                {
                    for( int shift = 0; shift < m; shift++ )
                    {
                        uint32_t s = (S >> shift) % ( 1 << n );
                        auto firstBits = (s & mask) ^ adjF;

                        if ( desiredBits == popcnt( firstBits ) )
                        {
                            leadingZeros[shift] = leadingZeros[shift] + 1;
                        }
                        else
                        {
                            break;
                        }
                    }
                }
            }

            memcpy( out.data() + (threadId * m), leadingZeros.data(), sizeof( leadingZeros[0] ) * m );
        } ));

    };

    std::array< uint32_t, m > leadingZeros;
    for( auto& x : leadingZeros )
        x = 0;

    for( auto& thread : threads )
    {
        thread.wait();
    }

    for( int i = 0; i < (threadCount * m); i++ )
    {
        leadingZeros[i % m] += out[i];
    }


    for( auto x : leadingZeros )
        std::cout << x << ", ";

    std::cout << std::endl;
}

int main()
{
    typedef std::chrono::high_resolution_clock clock;
    int rounds = 100;

    // do some rounds to get the cpu up to speed..
    for( int i = 0; i < rounds / 10; i++ )
    {
        convolve();
    }


    auto start = clock::now();

    for( int i = 0; i < rounds; i++ )
        convolve();

    auto end = clock::now();
    double seconds = std::chrono::duration_cast< std::chrono::microseconds >( end - start ).count() / 1000000.0;

    std::cout << seconds/rounds*1000 << " msec/round" << std::endl;

    return 0;
}

Salida de muestra:

   6317312, 2515072, 1034368, 434048, 190144, 85200, 39804, 19168,
   6226944, 2481408, 1031168, 438080, 192896, 86816, 40484, 19490,
   6321152, 2524672, 1045376, 442880, 195680, 88464, 41656, 20212,
   6330624, 2517504, 1031104, 430208, 187696, 83976, 38976, 18708,
   6304768, 2510336, 1030720, 433056, 190880, 86824, 40940, 19840,
   6272512, 2494720, 1028160, 432352, 189168, 84752, 39540, 19052,
   6233600, 2507520, 1046912, 447008, 198224, 89984, 42092, 20292,

C ++ x150 x450 x530

En lugar de matriz, utilicé bits (y magia oscura).
Gracias @ace por la función aleatoria más rápida.

Cómo funciona: los primeros 15 bits del número entero srepresentan la matriz S[15]; los ceros representan -1, los que representan +1. La matriz Fse construye de manera similar. Pero con dos bits para cada símbolo.

• 00 representa -1
• 01 y 10 representan 0
• 11 representan 1

Causar Sy Ftener una representación diferente que tengo que intercalar Sconsigo mismo para ser comparable con F.

• 0 (-1) se convirtió en 00 (-1 en la representación de F)
• 1 (+1) se convirtió en 11 (+1 en la representación de F)

Ahora podemos simplemente usar Carnot para calcular el producto interno. Recuerde que una variable solo puede asumir el valor 00 u 11

0. 00 = 11 (-1 * -1 = +1)
0. 01 = 10 (-1 * 0 = 0)
0. 10 = 01 (-1 * 0 = 0)
0. 11 = 00 (-1 * +1 = -1)
1. 00 = 00 (+1 * -1 = -1)
1. 10 = 10 (+1 * 0 = 0)
1. 01 = 01 (+1 * 0 = 0)
1. 11 = 11 (+1 * +1 = +1)

Parece un no xor para mí. :)

Sum the ones es solo un juego de cambio y máscara, nada realmente complejo.

#include <array>
#include <ctime>

// From standford bithacks
// http://graphics.stanford.edu/~seander/bithacks.html
inline int32_t interleaveBit(int32_t x)
{
   static const uint32_t B[] = { 0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF };
   x = (x | ( x << 8)) & B[3];
   x = (x | ( x << 4)) & B[2];
   x = (x | ( x << 2)) & B[1];
   x = (x | ( x << 1)) & B[0];
   return x | (x << 1);
}

inline int32_t sumOnes(int32_t v)
{
   static int b[] = { 1, 0, 0, 1};
   int s = 0;
   for( int i = 0; i < 8; ++i)
   {
      const int a = 3&(v>>(i*2));
      s += b[a];
   }
   return s;
}

inline int32_t sumArray(int32_t v)
{
   static int b[] = { -1, 0, 0, 1};
   int s = 0;
   for( int i = 0; i < 8; ++i)
   {
      const int a = 3&(v>>(i*2));
      s += b[a];
   }
   return s;
}

uint32_t x, y = 24252, z=57768, w=1564; //PRNG seeds

int32_t myRand()
{
   uint32_t t;
   t = x ^ (x<<1);
   x = y;
   y = z;
   z = w;
   w = w ^ ( w >> 19) ^ t ^ (t >> 8);
   return w;
}

int main()
{
   std::array<int, 8> leadingZero{0};
   x = static_cast<int32_t>(time(nullptr)); // seed PRNG
   const int maxS = 1 << 15;
   for(int s = 0; s < maxS; ++s)
   {
      const int32_t x = interleaveBit(s);
      for(int i = 0; i < 1000; ++i)
      {
         int32_t random;
         do
         {
            random = 0xFFFF & myRand();
         }while(sumOnes(random) == 0);
         int j = 7;
         while( j >= 0 )
         {
            const int32_t h = (x >> (j*2));
            const int32_t l = 0xFFFF & (~(random ^ h)); // inner product
            if(sumArray(l) == 0)
            {
               leadingZero[j]++;
            } else
            {
               break;
            }
            j--;
         }

      }
   }
   for(int i = 7; i >= 0; --i)
   {
      printf("%d ", leadingZero[i]);
   }
   printf("\n");
   return 0;
}

Aquí una salida de muestra:

6332350 2525218 1041716 438741 192917 87159 41023 19908 

real 0m0.372s
user 0m0.371s
sys  0m0.001s

El programa ha sido compilado con:

gcc -std=c++11 -O3 -msse4.2 -Wall -lstdc++ 26371.cpp -o fastPy

en Fedora 20 con gcc 4.8.2 La CPU es un i7 8core.

Probablemente pueda obtener algunos parámetros de compilación de ajuste de ms.

Si bien este es el tiempo de solución OP en mi máquina:

time pypy 26371.py
[6330609, 2523914, 1040885, 439303, 192708, 86987, 40710, 19498]

real 0m53.061s
user 0m53.016s
sys  0m0.022s

Editar:

Simplemente agregando openmp y cambiar el orden de la para tengo una ganancia x3, lo que lleva a una mejora del rendimiento x450 frente al código OP. : D En este caso, la leadingZeromatriz debe ser atómica. Los globales aleatorios ... son aleatorios, serán más aleatorios.

 #pragma omp parallel for
 for(int i = 0; i < 1000; ++i)
 {
    int32_t random;
    do
    {
       random = 0xFFFF & myRand();
    }while(sumOnes(random) == 0);
    for(int s = 0; s < maxS; ++s)
    {
       const int32_t x = interleaveBit(s);
       int j = 7;
       while( j >= 0 )
       {
          const int32_t h = (x >> (j*2));
          const int32_t l = 0xFFFF & (~(random ^ h)); // inner product
          if( sumArray(l) == 0 )
          {
             leadingZero[j]++;
          } else
          {
             break;
          }
          j--;
       }
    }
 }

necesita agregar -fopenmpa la bandera del compilador

Editar: 2 Como sugerente por el usuario71404 Cambié las funciones sumOnes y sumArray y ahora es súper rápido.

real  0m0.101s
user  0m0.101s
sys   0m0.000s

Con openmp es más lento, porque los atómicos agregan demasiada sobrecarga.

real  0m0.253s
user  0m1.870s
sys   0m0.001s

Sin atómica es aún más rápido, pero obtengo un resultado incorrecto.

2137992 1147218 619297 321243 155815 70946 32919 15579

real   0m0.048s
user   0m0.338s
sys    0m0.001s

Para comprender sumArray, considere que 16 bits representan y una matriz de 8 números.
00 no tiene 1 y representa -1
01 y 10 tiene un 1 y representa 0
11 tiene dos 1s y representa 1
Para que el recuento incorporado el número de bits establecido en 1 [ http://en.wikipedia.org/wiki/ Hamming_weight] y a cada grupo eliminamos 1. Cool.

sumOnes es solo magia negra.

Aquí las últimas compilaciones de banderas y códigos.

gcc -std = c ++ 11 -mfpmath = sse -O3 -flto -march = native -funroll-loops -Wall -lstdc ++

#include <cstdint>
#include <cstdio>
#include <ctime>

inline int32_t interleaveBit(int32_t x)
{
   static const uint32_t B[] = { 0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF };
   x = (x | ( x << 8)) & B[3];
   x = (x | ( x << 4)) & B[2];
   x = (x | ( x << 2)) & B[1];
   x = (x | ( x << 1)) & B[0];
   return x | (x << 1);
}

inline int32_t sumOnes(int32_t v)
{
   /* 0xAAAA == 0b1010 1010 1010 1010 */
   return !!(0xAAAA & (v ^ ~(v << 1)));
}

inline int32_t sumArray(int32_t v)
{
   return __builtin_popcount(v) - 8;
}

uint32_t x, y = 24252, z = 57768, w = 1564; //PRNG seeds

int32_t myRand()
{
   uint32_t t;
   t = x ^ (x << 1);
   x = y;
   y = z;
   z = w;
   w = w ^ ( w >> 19) ^ t ^ (t >> 8);
   return w;
}

int main()
{
   int leadingZero[8] = { 0 };
   x = static_cast<int32_t>(time(nullptr)); // seed PRNG
   const int maxS = 1 << 15;
   for( int i = 0; i < 1000; ++i )
   {
      int32_t random;
      do
      {
         random = 0xFFFF & myRand();
      } while(sumOnes(random) == 0 );
      for( int s = 0; s < maxS; ++s )
      {
         const int32_t x = interleaveBit(s);
         int j = 7;
         while( j >= 0 )
         {
            const int32_t h = (x >> (j * 2));
            const int32_t l = 0xFFFF & (~(random ^ h)); // inner product
            if( sumArray(l) == 0 )
            {
               leadingZero[j]++;
            } else
            {
               break;
            }
            j--;
         }
      }
   }
   printf("[%d, %d, %d, %d, %d, %d, %d, %d]\n",
      leadingZero[7], leadingZero[6],
      leadingZero[5], leadingZero[4],
      leadingZero[3], leadingZero[2],
      leadingZero[1], leadingZero[0]);
   return 0;
}

Julia: 0.7s, 120x más rápido

Como demostró el usuario 20768, un puerto directo del código para Julia es aproximadamente el doble de rápido que PyPy. Pero podemos hacer mucho mejor que eso.

function pleadingzerocounts(; n = 8,
                              m = 8,
                              iters = 1000)
  @parallel (+) for S = 1:2^(8+8-1)
    leading_counts = zeros(Int, m)
    F = Array(Int, n)
    for i = 1:iters
      flag = 0
      while flag == 0
        for i = 1:n
          v = (1-(rand(Int8)&3))%2
          @inbounds F[i] = v
          flag += v & 1
        end
      end
      for j = 1:m
        sum = 0
        for i = 1:n
          @inbounds sum += S & (1 << (j + i - 2)) > 0 ? F[i] : -F[i]
        end
        sum == 0 ?
          (leading_counts[j] += 1) :
          break
      end
    end
    leading_counts
  end
end

function main()
  # Warm up the JIT
  pleadingzerocounts()
  # Then go for real
  println(@time pleadingzerocounts())
end

Puede ejecutar esto usando julia -p 8 -e 'require("golf.jl");main()'(el 8 es el número de procesos, es posible que desee jugar con él). En la última versión preliminar de Julia, esto toma 0.7s vs. 1m22s para PyPy.

Si tiene suficientes núcleos en su computadora y tal vez active algunas instancias de AWS, debería poder afeitarse un poco más :)

C, 1.210s

Con el código de OP corriendo 1m45.729s en mi máquina.

Compilacion:

gcc -O3 -march=native -fwhole-program -fstrict-aliasing -ftree-vectorize -Wall ./test2.c -o ./test2

Un agradecimiento especial: @dyp para compilar banderas e ideas para optimizaciones.

#include <stdio.h>
#include <time.h>

#define n (8)
#define m (8)
#define iters (1000)
int leadingzerocounts[m]; // declared as global so initialised to 0
unsigned int x,y=34353,z=57768,w=1564; //PRNG seeds

/* xorshift PRNG
 * Taken from https://en.wikipedia.org/wiki/Xorshift#Example_implementation
 * Used under CC-By-SA */
int myRand() {
    unsigned int t;
    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = w ^ (w >> 19) ^ t ^ (t >> 8);
}

int dotproduct(int*F, int*S) {
    unsigned int i;
    int sum=0;
    for(i=0; i<n; i++) {
        sum+=F[i]*S[i];
    }
    return sum;
}

int main() {
    unsigned int i, j, tmp;
    x=(int)time(NULL); //seed PRNG

    int S[n+m-1];
    for(i=0; i<(1<<(n+m-1)); i++) {
        tmp=i;
        for(j=0; j<n+m-1; j++) {
            S[j]=(tmp&1)*(-2)+1;
            tmp>>=1;
        }
        for(j=0; j<iters; j++) {
            int F[n];
            unsigned int k, flag=0;
            do {
                for(k=0; k<n; k++) {
                    F[k]=(1-(myRand()&3))%2;
                    flag+=(F[k]&1);
                }
            } while(!flag);
            for(k=0; k<m&&!dotproduct(F, S+k); k++) {
                leadingzerocounts[k]++;
            }
        }
    }
    for(i=0; i<m; i++) printf("%d ", leadingzerocounts[i]);
    return 0;
}

Salida de muestra:

6334411 2527506 1042239 439328 192914 87005 40847 19728

Perl

Esto no es tan rápido como la solución C, pero creo que es bastante rápido para un lenguaje interpretado de alto nivel. Reduce aproximadamente el 40% del tiempo de ejecución de la implementación de Python.

#!/usr/bin/env perl

use v5.10;
use strict;
use warnings;
use Algorithm::Combinatorics qw( variations_with_repetition );
use List::Util qw( any sum );

use constant {
  N        => 8,
  M        => 8,
  ITERS    => 1000,
};

my @leadingzerocounts;

my $variations = variations_with_repetition([-1, 1], N + M - 1);
while (my $S = $variations->next)
{
  for my $i (1 .. ITERS)
  {
    my @F;
    until (@F and any { $_ } @F)
    {
      @F = map +((-1,0,0,1)[rand 4]), 1..N;
    }

    my $j = 0;
    while ($j < M)
    {
      last if sum map $F[$_]*$S->[$j+$_], 0..N-1;
      $leadingzerocounts[$j++]++;
    }
  }
}

say join ", ", @leadingzerocounts;

El algoritmo :: Combinatoria está disponible en Ubuntu ( sudo apt-get install libalgorithm-combinatorics-perl). Los otros módulos utilizados son módulos centrales de Perl, por lo que ya deberían estar instalados como parte de la instalación básica de Ubuntu.

Julia: ¡4.66 veces más lento!

Realmente estoy empezando a dudar de las estadísticas en su sitio web ...

Tenga en cuenta que el siguiente código de Julia es efectivamente una transcripción directa del código de Python del OP sin ninguna optimización. Uso la time()función para excluir el tiempo de inicio lento de Julia ...

srand(27182818284590)
t = time()

require("Iterators")

n = 8
m = 8
iters = 1000
bothzero = 0
leadingzerocounts = zeros(m)

for S in Iterators.product(fill([-1,1], n+m-1)...)
    for i = 1:iters
        F = [[-1 0 0 1][rand(1:4)] for j = 1:n]
        while all((x) -> x == 0, F)
            F = [[-1 0 0 1][rand(1:4)] for j = 1:n]
        end
        j = 1
        while j <= m && sum(map(*, F, S[j:j+n-1])) == 0
            leadingzerocounts[j] += 1
            j += 1
        end
    end
end

println(leadingzerocounts)

t = time() - t
println("$t seconds")

Julia: 5 m 32.912 s

Código de OP en PyPy: 1 m 11.506 s

Salida de Julia:

6332170
2525472
1041522
438761
193119
86873
40705
19662

RPython 0.187s (258x más rápido)

Fuente original con PyPy2.2.1: 1m 6.718s

Ahora con subprocesos, se ha eliminado el respaldo para Python estándar. El número de subprocesos de trabajo se puede especificar como un parámetro de línea de comando, el valor predeterminado es dos.

from time import time, sleep
from math import fmod

from rpython.rlib.rthread import start_new_thread, allocate_lock, get_ident
class Random:
  __slots__ = ['s']

  def __init__(self, s=1):
    self.s = s

  def init_genrand(self, seed):
    self.s = seed

  def genrand32(self):
    # xorshift PRNG with period 2^32-1
    # see http://core.kmi.open.ac.uk/download/pdf/6250138.pdf
    self.s ^= (self.s << 13)
    self.s ^= (self.s >> 17)
    self.s ^= (self.s << 5)
    return self.s

class ThreadEnv:
  __slots__ = ['n', 'm', 'iters', 'counts', 'running', 'lock']

  def __init__(self):
    self.n = 8
    self.m = 8
    self.iters = 1000
    self.counts = [0]*8
    self.running = 0
    self.lock = None

env = ThreadEnv()
truth = [-1,0,0,1]

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

  if argc == 5:
    num_threads = int(argv[4])
  else:
    num_threads = 2

  env.n = int(argv[1])
  env.m = int(argv[2])
  env.iters = int(argv[3]) // num_threads
  env.counts = [0]*env.m
  env.lock = allocate_lock()

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

  env.running += 1

  # use the main process as a worker
  run()

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

  print env.counts
  return 0

def run():
  n, m, iters = env.n, env.m, env.iters
  counts = [0]*m
  sbits = [0]*(n+m-1)

  random = Random()
  seed = int(fmod(time(), 2147483.648)*1000) ^ get_ident()
  random.init_genrand(seed)

  for S in xrange(1<<n+m-1):
    i = 0
    sbit = 0
    while not sbit:
      sbits[i] ^= 3
      sbit = sbits[i]
      i += 1

    for i in xrange(iters):
      f = 0
      while not f:
        F = random.genrand32()

        G, x = F, 0
        for k in xrange(n):
          x += truth[(G&3)^sbits[k]]
          f |= x
          G >>= 2

      if not x:
        counts[0] += 1
        for j in xrange(1, m):
          G, x = F, 0
          for k in xrange(j, n+j):
            x += truth[(G&3)^sbits[k]]
            G >>= 2
          if x: break
          counts[j] += 1

  # passing True stalls until a lock can be obtained
  env.lock.acquire(True)

  for i in xrange(m):
    env.counts[i] += counts[i]
  env.running -= 1

  env.lock.release()

def target(*args):
  return main, None

RPython es un subconjunto restringido de Python, que puede traducirse a C y luego compilarse utilizando la cadena de herramientas RPython . Su propósito expreso es ayudar en la creación de intérpretes de idiomas, pero también se puede utilizar para compilar programas simples como el anterior. La mayoría de las características 'más elegantes' de Python, como itertoolso incluso mapno están disponibles.

Para compilar, haga un clon local del repositorio actual de pypy y ejecute lo siguiente:

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

El ejecutable resultante se llamará convolution-c o similar en el directorio de trabajo actual.

He parametrizado las variables de entrada, por lo que el programa debe ejecutarse como:

convolution-c 8 8 1000

para que coincida con el código de muestra.

Notas de implementacion

S in itertools.product([-1,1], repeat = n+m-1)se convierte S in xrange(1<<n+m-1), interpretando Scomo un mapa de bits: [ 0, 1] → [ -1,1 ]

Del mismo modo, Fes también un mapa de bits, con cada dos bits que representan un único valor:
[ 00, 01, 10, 11] → [ -1, 0, 0,1 ]

Se utiliza una tabla de verdad para buscar el producto, en lugar de realizar una aplicación múltiple.

Debido a que se utilizan enteros con signo de 32 bits, no npueden ser mayores que 15 n+mni mayores que 31. Se puede lograr un soporte arbitrario de enteros conrpython.rlib.rbigint Si es necesario, módulo.

La primera iteración del bucle dot-producto se desenrolla, y se combina con la prueba de nulidad de F.

Se utiliza un PRNG casero, la fuente está en la lista. El autor del artículo demuestra un período de 2 ³² -1, y afirma que pasa todas las pruebas de Diehard, excepto una, aunque personalmente no lo he confirmado.

La semilla aleatoria cambia cada milisegundo, lo que es tan bueno como lo permita el uso de una marca de tiempo. Además, cada subproceso de trabajo envía xorsu identificación de proceso con este valor, para garantizar que cada uno tenga una semilla diferente.

Tiempos de muestra

2 hilos de trabajo:

$ timeit convolution-c 8 8 1000 2
[6331845, 2526161, 1042330, 440018, 193724, 87147, 40943, 19603]

Elapsed Time:     0:00:00.375
Process Time:     0:00:00.687
System Calls:     6927

4 hilos de trabajo:

$ timeit convolution-c 8 8 1000 4
[6334565, 2527684, 1043502, 440216, 193225, 87398, 40799, 19338]

Elapsed Time:     0:00:00.218
Process Time:     0:00:00.796
System Calls:     3417

8 hilos de trabajo:

$ timeit convolution-c 8 8 1000 8
[6327639, 2522483, 1039869, 437884, 192460, 86771, 40420, 19403]

Elapsed Time:     0:00:00.187
Process Time:     0:00:00.734
System Calls:     3165

Fuente original de OP:

$ timeit pypy convolution-orig.py
[6330610, 2525644, 1041481, 438980, 193001, 86622, 40598, 19449]

Elapsed Time:     0:01:06.718
Process Time:     0:01:06.718
System Calls:     11599808

Tiempo para 100000 iteraciones:

$ timeit convolution-c 8 8 100000 8
[633156171, 252540679, 104129386, 43903716, 19307215, 8709157, 4072133, 1959124]

Elapsed Time:     0:00:16.625
Process Time:     0:01:02.406
System Calls:     171341

Julia: 1 min 21.4s (2.2x más rápido) (modificación del código de Arman)

Código de operación en PyPy: 3 min 1.4 s

Ambos hechos en el REPL, sin incluir el tiempo para cargar paquetes.

function foo()                                                                                                                                                             
    n = 8                                                                                                                                                                  
    m = 8                                                                                                                                                                  
    iters = 1000                                                                                                                                                           
    bothzero = 0                                                                                                                                                           
    leadingzerocounts = zeros(Int,m)                                                                                                                                       
    P=[-1,0,0,1]                                                                                                                                                           

    for S in Iterators.product(fill([-1,1], n+m-1)...)                                                                                                                     
        Sm=[S...]                                                                                                                                                          
        for i = 1:iters                                                                                                                                                    
            F = P[rand(1:4,n)]                                                                                                                                             
            while all(F==0)                                                                                                                                                
                F = P[rand(1:4,n)]                                                                                                                                         
            end                                                                                                                                                            
            j = 1                                                                                                                                                          

            while j <= m && dot(F,Sm[j:j+n-1]) == 0                                                                                                                        
                leadingzerocounts[j] += 1                                                                                                                                  
                j += 1                                                                                                                                                     
            end                                                                                                                                                            
        end                                                                                                                                                                
    end                                                                                                                                                                    

    println(leadingzerocounts)                                                                                                                                             
end 

Hay algunos problemas con el código de Arman que lo hace muy lento: utiliza muchas funciones anónimas y funciones de orden superior innecesariamente. Para probar si todo un vector F es cero, ¿por qué no escribir todo (F == 0) en lugar de todo (x-> x == 0, F)? Es más corto y literalmente mil veces más rápido.

También usa sum (map (*, x, y)) como producto de punto en lugar de simplemente punto (x, y). La primera versión es 650 veces más lenta para un vector de 10k dobles. Y la función del producto punto se implementa como un bucle for en Julia pura.

Además, las comprensiones de matriz son lentas. Es mejor escribir [0,1,0, -1] [rand (1: 4, n)] en lugar de [[-1 0 0 1] [rand (1: 4)] para j = 1: n] .

Finalmente, las variables globales son malas juju en Julia. Julia solo es rápida si escribe código de tal manera que permita que el JIT y la inferencia de tipos funcionen. Una gran parte de esto es la estabilidad de tipo: el compilador debe poder asegurarse de que el tipo de una variable no cambiará dentro de un ciclo, por ejemplo.

Nimrod

import times, locks, strutils, unsigned

const
  N = 8
  M = 8
  iters = 1000
  numThreads = 8

type
  SVec = array[0..N+M-1, int]
  FVec = array[0..N-1, int]
  ComputeThread = TThread[int]

var
  rngSeed = int(epochTime()*1000)
  totalLeadingZeros: array[0..M-1, int]
  lock: TLock

type
  RNGState = object
    x, y, z, w: uint32

proc newRNG(seed: int): RNGState =
  result.x = uint32(seed)

proc random(rng: var RNGState): int =
  let t = rng.x xor (rng.x shl 11)
  rng.x = rng.y; rng.y = rng.z; rng.z = rng.w
  rng.w = rng.w xor (rng.w shr 19) xor t xor (t shr 8)
  result = int(rng.w)

proc initVecRand(v: var FVec, rng: var RNGState) =
  const values = [ -1, 0, 0, 1 ]
  var rnd = rng.random
  var bitAcc = 0
  for i in 0 .. <len(v):
    let val = values[rnd and 3]
    rnd = rnd shr 2
    v[i] = val
    bitAcc = bitAcc or val
  if bitAcc == 0:
    initVecRand(v, rng)

proc convolve(s: SVec, f: FVec, offset: int): int =
  for i in 0 .. <len(f):
    result += s[i+offset]*f[i]

proc iterate(v: var SVec) =
  for i in 0 .. <len(v):
    if v[i] == -1:
      v[i] = 1
      return
    v[i] = -1

proc mainThread(id: int) {.thread.} =
  const numS = 1 shl (N+M-1)
  var
    s: SVec
    f: FVec
    leadingZeros: array[0..M-1, int]
    rng = newRNG(rngSeed + id)
  for k in 0 .. <len(s):
    s[k] = -1
  for i in 1..numS:
    for j in countUp(id, iters, numThreads):
      initVecRand(f, rng)
      if convolve(s, f, 0) == 0:
        leadingZeros[0] += 1
        for k in 1 .. <M:
          if convolve(s, f, k) == 0:
            leadingZeros[k] += 1
          else:
            break
    iterate(s)
  acquire(lock)
  for i in 0 .. <M:
    totalLeadingZeros[i] += leadingZeros[i]
  release(lock)

proc main =
  let startTime = epochTime()
  var threads: array[1..numThreads, ComputeThread]
  initLock(lock)
  for i in 1..numThreads:
    createThread(threads[i], mainThread, i)
  for i in 1..numThreads:
    joinThread(threads[i])
  echo("Leading zeros: ", @totalLeadingZeros)
  let endTime = epochTime()
  echo("Time taken:    ", formatFloat(endTime - startTime, ffDecimal, 3),
       " seconds")

main()

Salida de ejemplo:

Leading zeros: @[6333025, 2525808, 1042466, 439138, 192391, 86751, 40671, 19525]
Time taken:    0.145 seconds

Nimrod compila a C, por lo tanto, la elección del compilador de C para el backend también es importante.

Usando clang, compila con:

nimrod cc --threads:on --cc=clang --passc:-flto -d:release conv.nim

Usando gcc, compila con:

nimrod cc --threads:on --cc=gcc --passc:-flto -d:release conv.nim

Omita --passc:-fltosi tiene un compilador de C más antiguo que no admite LTO. Omita la --cc=...opción si está de acuerdo con la opción predeterminada para el compilador de C. El código requiere Nimrod 0.9.4 o 0.9.5 .

En mi iMac quadcore (2.66 GHz core i5), el código se ejecuta en aproximadamente .15 segundos con gcc 4.9, .16 segundos con clang, en comparación con 88 segundos para PyPy 2.2.1 (es decir, una aceleración de más de 500 veces). Desafortunadamente, no tengo acceso a una máquina con más de cuatro núcleos que también tenga instalado PyPy o donde pueda instalar PyPy fácilmente, aunque obtengo aproximadamente 0,1 segundos (con mucho ruido de medición) en un AMD de 64 núcleos Opteron 6376 1.4 GHz (según / proc / cpuinfo) con gcc 4.4.6.

La implementación intenta ser fiel al código original en lugar de optimizar el código a costa de la legibilidad, sin renunciar a optimizaciones obvias. Curiosamente, la recursión de cola initVecRand()es un poco más rápida que un bucle con una instrucción de interrupción con gcc y clang. Desenrollar manualmente una iteración del convolvebucle de prueba dentro del bucle principal también produjo una aceleración, presumiblemente debido a una mejor predicción de ramificación.

Java

Traduje la solución C ++ anterior a Java:

import java.util.Random;
import java.util.Arrays;

public class Bench2 {
  public static int[] bits = { 0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF };
  public static int[] oneValues = { 1, 0, 0, 1 };
  public static int[] values = { -1, 0, 0, 1 };
  public static int n = 8;
  public static int m = 8;
  public static int iters = 1000;

  private static int x,y=34353,z=57768,w=1564;

  public static void main( String[] args ) {
    x = (int) (System.currentTimeMillis()/1000l);

    int[] leadingzerocounts = new int[ m ];
    Arrays.fill( leadingzerocounts, 0 );

    int maxS = 1 << 15;

    for( int s = 0; s < maxS; s++ ) {
      int x = interleaveBit( s );

      for( int i=0; i<iters; i++ ) {
        int random;

        do {
          random = 0xFFFF & fastRandom( );
        } while( sumOnes( random ) == 0 );

        int j = 7;

        while( j >= 0 ) {
          int h = ( x >> (j*2) );
          int l = 0xFFFF & (~(random ^ h));

          if( sumArray( l ) == 0 ) {
            leadingzerocounts[ j ]++;
          } else {
            break;
          }

          j--;
        }
      }
    }

    for( int i = 7; i >= 0; --i ) {
      System.out.print( leadingzerocounts[ i ] + " " );
    }

    System.out.println( );
  }

  public static int interleaveBit( int x ) {
    x = (x | ( x << 8)) & bits[3];
    x = (x | ( x << 4)) & bits[2];
    x = (x | ( x << 2)) & bits[1];
    x = (x | ( x << 1)) & bits[0];
    return x | (x << 1);
  }

  public static int sumOnes( int v ) {
    return (0xAAAA & (v ^ ~(v << 1)));
    // int s = 0;

    // for( int i = 0; i < 8; ++i ) {
    //   int a = 3 & ( v >> (i*2) );
    //   s += oneValues[ a ];
    // }

    // return s;
  }

  public static int sumArray( int v ) {
    return Integer.bitCount( v ) - 8;
    // int s = 0;

    // for( int i=0; i<8; ++i ) {
    //   int a = 3 & ( v >> (i*2) );
    //   s += values[ a ];
    // }

    // return s;
  }

  public static int fastRandom( ) {
    long t;
    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = (int)( w ^ (w >> 19) ^ t ^ (t >> 8));
  }
}

En mi máquina obtengo el siguiente resultado para el programa java:

time java Bench2
6330616 2524569 1040372 439615 193290 87131 40651 19607 
java Bench2  0.36s user 0.02s system 102% cpu 0.371 total

El programa OP se ejecuta unos 53 segundos en mi máquina:

time pypy start.py
[6330944, 2524897, 1040621, 439317, 192731, 86850, 40830, 19555]
pypy start.py  52.96s user 0.06s system 99% cpu 53.271 total

El programa c ++ se ejecutó solo unos 0,15 segundos:

time ./benchcc
[6112256, 2461184, 1025152, 435584, 193376, 87400, 40924, 19700]
./benchcc  0.15s user 0.00s system 99% cpu 0.151 total

Eso es aproximadamente 2.5 veces más rápido que la solución java correspondiente (no excluí el inicio de VM). Esta solución de Java es aproximadamente 142 veces más rápida que el programa ejecutado con PyPy.

Como estaba personalmente interesado, configuré iters 100_000 para Java y C ++, pero el factor de 2.5 no disminuyó a favor de Java si algo aumentaba.

EDITAR: Ejecuté los programas en una PC Arch Linux de 64 bits.

EDIT2: Quiero agregar que comencé con una traducción aproximada del código de Python:

import java.util.Random;
import java.util.Arrays;

public class Bench {
    public static int[] values = { -1, 0, 0, 1 };
    public static int n = 8;
    public static int m = 8;
    public static int iters = 1000;

    private static int x,y=34353,z=57768,w=1564; 

    public static void main( String[] args ) {
        x = (int) (System.currentTimeMillis()/1000l);

        int[] leadingzerocounts = new int[ m ];
        Arrays.fill( leadingzerocounts, 0 );

        int[] S = new int[ n+m-1 ];
        Arrays.fill( S, -1 );

        do {
            for( int i=0; i<iters; i++ ) {
                int[] F = new int[ n ];

                do {
                    randomArray( F );
                } while( containsOnlyZeros( F ) );

                for( int j=0; j < m && check( F, S, j ); j++ ) {
                    leadingzerocounts[ j ] += 1;
                }
            }
        } while( next( S ) );

        System.out.println( Arrays.toString( leadingzerocounts ) );
    }

    public static void randomArray( int[] F ) {
        for( int i = 0; i<F.length; i++ ) {
            F[ i ] = (1-(fastRandom()&3))%2;
        }
    }

    public static boolean containsOnlyZeros( int[] F ) {
        for( int x : F ) {
            if( x != 0 ) {
                return false;
            }
        }

        return true;
    }

    public static boolean next( int[] S ) {
        for( int i=0; i<S.length; i++ ) {
            if( ( S[ i ] = -S[ i ] ) == 1 ) {
                return true;    
            }
        }

        return false;
    }

    public static boolean check( int[] F, int[] S, int j ) {
      int sum = 0;

      for( int i=0; i<n; i++ ) {
          sum += F[ i ] * S[ j + i ];
      }

      return sum == 0;
    }

    public static int fastRandom( ) {
        long t;
        t = x ^ (x << 11);
        x = y; y = z; z = w;
        return w = (int)( w ^ (w >> 19) ^ t ^ (t >> 8));
    }
}

Este programa ejecutó aproximadamente 3.6 segundos:

time java Bench   
[6330034, 2524369, 1040723, 439261, 193673, 87338, 40840, 19567]
java Bench  3.64s user 0.01s system 101% cpu 3.600 total

Que es aproximadamente 14 veces más rápido que la solución PyPy. (Elegir la función aleatoria estándar sobre la función fastRandom lleva a un tiempo de ejecución de 5 segundos)

Python 3.5 + numpy 1.10.1, 3.76 segundos

Las pruebas se ejecutaron en mi Macbook Pro. El código de OP tomó ~ 6 minutos en la misma máquina.

La razón por la que estoy respondiendo esta pregunta es porque no tengo 10 reputaciones y no puedo responder la Parte I :-p

Durante los últimos días, he estado tratando de descubrir cómo realizar convoluciones masivas de manera eficiente con numpy (sin depender de un paquete de terceros, incluso scipy). Cuando me encontré con esta serie de desafíos durante mi investigación, decidí intentarlo. Puede que haya llegado a este juego demasiado tarde, pero aquí está mi intento de usar Python 3.5 y numpy 1.10.1.

def test_convolv():
    n = 8 
    m  = 8 
    iters = 1000
    ilow = np.ceil(0+n/2).astype(int)
    ihigh = np.ceil(m+n/2).astype(int)

    leadingzerocounts = np.zeros(m)

    # Pre-compute S & F
    S = np.array(list(itertools.product([-1,1], repeat = n+m-1)))
    choicesF = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n*iters).reshape(iters,n)
    imask = ~np.any(choicesF, axis=1)
    while np.any(imask):
        imasksize = np.count_nonzero(imask)
        choicesF[imask,:] = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n*imasksize).reshape(imasksize, n)
        imask = ~np.any(choicesF, axis=1)

    for i in np.arange(iters):
        F = choicesF[i, :]
        # This is where the magic is: by flattening the S array, 
        # I try to take advantage of speed of the np.convolve 
        # (really numpy.multiarray.correlate). 
        FS = (np.convolve(S.reshape(-1), F, 'same').reshape(S.shape))[:, ilow:ihigh]
        jmask_not = (FS[:, 0] != 0)
        leadingzerocounts[0] = leadingzerocounts[0]+np.count_nonzero(~jmask_not)
        for j in np.arange(n-1)+1:
            jmask = (FS[jmask_not, j] != 0)
            leadingzerocounts[j] = leadingzerocounts[j] + np.count_nonzero(~jmask)
            jmask_not[(jmask_not.nonzero()[0])[jmask]] = False

    print(leadingzerocounts)

Precalculé las matrices S y F, y aplané la matriz S mientras realizaba la convolución, que (según mis experimentos) podría aprovechar la velocidad de np.convolve. En otras palabras, como no encontré una rutina de convolución vectorizada, falsifiqué el código vectorizando el conjunto completo y esperé que np.convolved hiciera la vectorización bajo el capó para mí, lo que parecía estar funcionando. Tenga en cuenta que utilicé mode = 'same' y recorté los elementos iniciales y finales que eran inútiles.

En mi Macbook Pro, los resultados de la prueba dan 3,76 segundos . Cuando ejecuté el código de OP (modificado a Python 3.5), obtuve unos 6 minutos . La aceleración es de aproximadamente 100 veces.

Un inconveniente es que debido a que las matrices S y F deben almacenarse, el requisito de memoria puede ser un problema si los tamaños son demasiado grandes.

Utilicé el mismo método para la Parte I y obtuve una aceleración de ~ 60-100x en mi computadora portátil.

Como hice todo en mi Macbook Pro, si alguien pudiera probar mi código y hacerme saber cómo funciona en su máquina, ¡lo agradecería mucho!

J, ^130x ~ 50x aceleración?

n =: m =: 8
len =: 1000
S =: (] - 0 = ])S0=: #:i.2^<:+/n,m
k =: (n#0) -.~ (_1 0 0 1) {~ (n#4) #: i.4^n
sn =: (]-0=])#:i.2^n
ku =: ~. k
M =: 0=+/"1 sn *"1/ ku
fs =: (ku&i.)"1 k
snum =: n #.\"1 S0

run =: 3 : 0
 r =: n#0
 for_t. (snum) do.
   rn =: fs{~? len # #k
   r =: r + +/"1*/\rn{"1 t{M
 end.
 r
)
echo run 0
exit''

Tiempos en un debian aleatorio:

u#>time j slowpy.ijs
6334123 2526955 1041600 440039 193567 87321 40754 19714

real    0m2.453s
user    0m2.368s
sys     0m0.084s


u#>time python slow_pyth.py
[6331017, 2524166, 1041731, 438731, 193599, 87578, 40919, 19705]

real    5m25.541s
user    5m25.548s
sys     0m0.012s

Creo que hay margen de mejora.

C ++: x200 (i7 de 4 núcleos, debe escalar a x400 en 8 núcleos)

Intentando una solución C ++ 11 más sencilla (probado con VS 2012, gcc y clang) con paralelización.

Para que esto se compile y se ejecute en Linux con gcc 4.8.1:

g ++ -O3 -msse -msse2 -msse3 -march = native -std = c ++ 11 -pthread -Wl, - golf.cpp no ​​es necesario

Bajo Linux también necesitamos std::launch::asyncforzar múltiples hilos. Me faltaba eso en una versión anterior.

En Visual Studio (2012+) esto debería funcionar, pero hacer una versión de lanzamiento para el tiempo ...

En mi viejo i3 dual core esto funciona en ~ 0.9 segundos. En mi i7 quad core, esto es 0.319s frente a pypy 66 segundos.

En un i7 de 8 núcleos, esto debería estar en el rango de aceleración x400. Cambiar a matrices de estilo C lo aceleraría, pero estaba interesado en permanecer con contenedores C ++. Para mí es interesante ver la aceleración que puede obtener mientras se mantiene relativamente cerca del dominio del problema y en un nivel relativamente alto, algo en lo que creo que C ++ es realmente bueno. También es de destacar la relativa facilidad de paralelización utilizando construcciones C ++ 11.

La solución de bits de @ ilmale es muy buena y funciona para -1/1/0. También se podría lanzar SSE a esto y tal vez obtener una aceleración significativa.

Más allá de la paralelización, hay otro "truco" que reduce el número de sumas. Resultados de muestra: 6332947 2525357 1041957 438353 193024 87331 40902 19649

#include <vector>
#include <iostream>
#include <thread>
#include <future>
#include <time.h>
#include <ctime>
#include <algorithm>

using namespace std;

// Bring some of these constants out to share
const size_t m = 8;
const int nthreads = 16;
const size_t cn = 15;
const int two_to_cn = 32768;

static unsigned int seed = 35;

int my_random() // not thread safe but let's call that more random!
{
   seed = seed*1664525UL + 1013904223UL; // numberical recipes, 32 bit
   return ((seed>>30)&1)-!!((seed>>30)&2); // Credit to Dave!
}

bool allzero(const vector<int>& T)
{
   for(auto x : T)
   {
      if(x!=0)
      {
         return false;
      }
   }
   return true;
}

// Return the position of the first non-zero element
size_t convolve_until_nonzero(size_t max_n, const vector<int>& v1, const vector<int>& v2)
{
   for(size_t i = 0; i<max_n; ++i)
   {
      int result = 0;
      for(size_t j = 0; j<v2.size(); ++j)
      {
         result += v1[i+j]*v2[j];
      }
      if(result!=0)
      {
         return i;
      }
   }
   return max_n;
}

void advance(vector<int>& v)
{
   for(auto &x : v)
   {
      if(x==-1)
      {
         x = 1;
         return;
      }
      x = -1;
   }
}

vector<int> convolve_random_arrays(vector<int> S, int range)
{
   const int iters = 1000;
   int bothzero = 0;
   int firstzero = 0;

   time_t current_time;
   time(&current_time);
   seed = current_time;


   vector<int> F(m);
   vector<int> leadingzerocounts(m+1);

   for(auto &x: leadingzerocounts)
   {
      x = 0;
   }

   for(int i=0; i<range; ++i)
   {
      for(int j=0; j<iters; ++j)
      {
         do
         {
            for(auto &x : F)
            {
               x = my_random();
            }
         } while(allzero(F));
         leadingzerocounts[convolve_until_nonzero(m, S, F)]++;
      }
      advance(S);
   }

   // Finish adding things up...
   for(int i=m-1; i>0; --i)
   {
      leadingzerocounts[i] += leadingzerocounts[i+1];
   }

   vector<int> withoutfirst(leadingzerocounts.begin()+1, leadingzerocounts.end());
   return withoutfirst;
}

int main(int argc, char* argv[])
{

   vector<int> leadingzerocounts(m);

   for(auto &x: leadingzerocounts)
   {
      x = 0;
   }

   clock_t start = clock();

   vector<int> S(cn);
   for(auto &x : S)
   {
      x = -1;
   }

   vector< future< vector< int > > > fs; // The future results of the threads

   // Go make threads to work on parts of the problem
   for(int i=0; i<nthreads; ++i)
   {
      vector<int> S_reversed = S; // S counts using LSBs but we want the thread start to be in MSBs
      reverse(S_reversed.begin(), S_reversed.end());
      fs.push_back(async(std::launch::async, convolve_random_arrays, S_reversed, two_to_cn/nthreads));
      advance(S);
   }
   // And now collect the data
   for(auto &f : fs)
   {
      vector<int> result = f.get();
      for(int i=0; i<result.size(); ++i)
      {
         leadingzerocounts[i] += result[i];
      }
   }

   for(auto count : leadingzerocounts)
   {
      cout << count << endl;
   }

   return 0;
}

Fortran: 316x

De acuerdo, Fortran: lo tengo hasta 106x 155x 160x 316x cuando uso un Xorshift RNG y OpenMP en una CPU i7 de 4 núcleos. Aparte de eso, no hay grandes trucos. Para que el iterador construya S, solo uso la representación binaria del entero de 16 bits i. Notarás que, aparte del RNG en línea y el "iterador" / mapeo de i a S, el código es tan de alto nivel como el código de Python.

Editar: eliminó el "if" en el Xorshift, ahora usando "r = abs (w / ...)" en lugar de "r = w / ...". Va de 106x a 155x.

Edit2: esto genera 15 veces más números aleatorios que la solución C ++. Si alguien tiene una solución de sobrecarga cero para convertir un int aleatorio en una matriz de 0s y 1s en Fortran, soy todo oídos. Entonces podríamos vencer a C ++ :)

Edit3: La primera edición introdujo un error, como señaló Lembik. Esto se soluciona ahora, con una pequeña mejora en la aceleración. Intentaré usar la sugerencia de Eelvex para obtener más velocidad.

Edit4: el perfil indicó que la conversión a real y de nuevo a entero con nint () fue lenta. Reemplacé esto con una división entera haciendo escala y redondeo, pasando de 160x a 316x de aceleración.

Compilar con:

gfortran -O3 -march = nativo -fopenmp golf.f90

program golf
implicit none
integer, parameter :: m=8, n=8
integer :: F(n), S(m+n-1), leadingzerocounts(m)
integer :: j,k,bindec,enc,tmp,x=123456789,y=362436069,z=521288629,w=88675123
integer*2 :: i
real :: r

leadingzerocounts=0

!$OMP parallel do private(i,enc,j,bindec,S,F,k,tmp,x,y,z,w,r) reduction(+:leadingzerocounts) schedule(dynamic)
do i=0,32766
  enc=i
  ! Short loop to convert i into the array S with -1s and 1s
  do j=16,2,-1
    bindec=2**(j-1)
    if (enc-bindec .ge. 0) then
      S(j-1)=1
      enc=enc-bindec
    else
      S(j-1)=-1
    endif
  end do
  do j=1,1000
    F=0
    do while (.not. any(F /= 0))
      do k=1,n
        ! Start Xorshift RNG
        tmp = ieor(x,ishft(x,11))
        x = y
        y = z
        z = w
        w = ieor(ieor(w,ishft(w,-19)),ieor(tmp,ishft(tmp,-8)))
        ! End Xorshift RNG
        ! Just scale it inside the nint:
        !F(k)=nint(w/2147483648.0)
        ! Scaling by integer division is faster, but then we need the random 
        ! number to be in (-2,2) instead of [-1,1]:
        F(k)=w/1073741824

      end do
    end do
    do k=1,m
      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

Salida de ejemplo:

$ time ./a.out
6329624 2524831 1039787 438809 193044 6860 40486 19517
./a.out 1.45s usuario 0.00s sistema 746% cpu 0.192 total

Código de OP:

$ time pypy golf.py
pypy golf.py 60.68s usuario 0.04s sistema 99% cpu 1: 00.74 total