Ceci est un suivi de la lenteur de Python. (Ou à quelle vitesse votre langue est-elle?) .

Il s’est avéré qu’il était un peu trop facile d’obtenir une accélération x100 pour ma dernière question. Pour la partie qui a aimé le défi mais qui veut quelque chose de plus difficile où ils peuvent vraiment utiliser leurs compétences de bas niveau, voici la partie II. Le défi consiste à obtenir une accélération x100 pour le code python suivant, testé sur mon ordinateur.

Pour rendre les choses plus difficiles, j'utilise pypy cette fois-ci. Pour moi, le timing actuel est de 1 minute et 7 secondes avec Pypy 2.2.1.

Règles

1. La première personne à soumettre du code que je peux exécuter est correcte et x100 fois plus rapide sur ma machine se verra attribuer une prime de 50 points.
2. Je vais attribuer la victoire au code le plus rapide après une semaine.

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 sortie doit être similaire à

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

Vous devez utiliser une graine aléatoire dans votre code et tout générateur de nombre aléatoire suffisant pour donner des réponses proches de ce qui précède sera accepté.

Ma machine Les timings seront exécutés sur ma machine. Il s’agit d’une installation ubuntu standard sur un processeur à huit cœurs AMD FX-8350. Cela signifie également que je dois pouvoir exécuter votre code.

Explication du code

Ce code itère sur tous les tableaux S de longueur n + m-1 constitués de -1 et de 1. Pour chaque tableau S, il échantillonne 1000 tableaux aléatoires non nuls F de longueur n constituée de -1,0 ou 1 avec une probabilité de 1/4, 1/2, / 14 de la prise de chaque valeur. Il calcule ensuite les produits internes entre F et chaque fenêtre de S de longueur n jusqu'à trouver un produit interne non nul. Il ajoute 1 leadingzerocountsà chaque position où il a trouvé un produit intérieur nul.

Statut

• Perl . Ralentissement de 2,7 fois par @tobyink. (Comparé à pypy pas cpython.)

• J . 39 fois plus rapide par @Eelvex.

• C . 59 fois plus rapide par @ace.
• Julia . 197 fois plus rapide sans compter le temps de démarrage de @ une minute de plus. 8,5 fois plus rapide, temps de démarrage compris (l'utilisation de 4 processeurs est plus rapide que 8).
• Fortran . 438 fois plus rapide en @ semi-extrinsèque.
• Rpython . 258 fois plus rapide par @primo.
• C ++ . 508 fois plus rapide par @ilmale.

(J'ai arrêté de chronométrer les nouvelles améliorations car elles sont trop rapides et trop petites.)

Il a été souligné que les durées inférieures à une seconde ne sont pas fiables et que certaines langues ont un coût de démarrage. L'argument est que si vous voulez inclure cela, vous devez également inclure le temps de compilation de C / C ++, etc. Voici les timings pour le code le plus rapide avec le nombre d'itérations porté à 100 000.

• Julia . 42 secondes par @ une minute de plus.
• C ++ . 14 secondes par @GuySirton.
• Fortran . 14s par @ semi-extrinsiques.
• C ++ . 12s par @ilmale.
• Rpython . 18s par @primo.
• C ++ . 5s de @Stefan.

Le gagnant est .. Stefan!

Défi de suivi affiché. À quelle hauteur pouvez-vous aller? (Un défi de codage + algorithmes) . Celui-ci est plus difficile.

C ++ bit magic

~ 16ms en multithread, 56ms en simple. ~ 4000 accélération.

(L'accélération est basée sur le code multithread sur mon i7-2820QM et les 1 min 9 secondes mentionnées dans la question. Étant donné que le système de l'OP a des performances plus médiocres que mon processeur, mais une meilleure perfiormance multi-thread, je m'attends à ce que ce nombre soit exact)

La partie multithread est assez inefficace en raison de la génération de threads. Je pourrais probablement faire mieux en exploitant ma bibliothèque de tâches personnalisée, mais celle-ci présente des bogues sous les systèmes unix. Pour une explication et un code presque identique sans thread, reportez-vous à la page https://codegolf.stackexchange.com/a/26485/20965 .

modifier

J'ai donné à chaque thread son propre RNG et ai réduit la longueur de bit à 32, ce qui a réduit le temps d'exécution de quelques 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;
}

Exemple de sortie:

   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

Au lieu de tableau, j'ai utilisé des bits (et de la magie noire).
Merci @ace pour la fonction aléatoire plus rapide.

Comment ça marche: les 15 premiers bits de l'entier sreprésentent le tableau S[15]; les zéros représentent -1, ceux-ci représentent +1. Le tableau Fest construit de manière similaire. Mais avec deux bits pour chaque symbole.

• 00 représente -1
• 01 et 10 représentent 0
• 11 représentent 1

Causer Set Favoir une représentation différente je dois entrelacer Savec lui-même pour être comparable F.

• 0 (-1) est devenu 00 (-1 dans la représentation de F)
• 1 (+1) est devenu 11 (+1 dans la représentation de F)

Maintenant, nous pouvons simplement utiliser Carnot pour calculer le produit intérieur. Rappelez-vous qu'une variable ne peut prendre que la valeur 00 ou 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)

On dirait un pas xor pour moi. :)

Somme les uns est juste un jeu de décalage et de masque, rien de vraiment complexe.

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

Voici un exemple de sortie:

6332350 2525218 1041716 438741 192917 87159 41023 19908 

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

Le programme a été compilé avec:

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

sur Fedora 20 avec gcc 4.8.2 Le Cpu est un i7 8core.

Je peux probablement gagner quelques ms en peaufinant les paramètres du compilateur.

Bien que ce soit le temps de solution OP sur ma machine:

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

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

Modifier:

Il suffit d’ajouter openmp et de modifier l’ordre des résultats car j’ai un gain x3, ce qui conduit à une amélioration des performances x450 par rapport au code OP. : D Dans ce cas, le leadingZerotableau doit être atomique. Les aléatoires globaux ... sont aléatoires, ils seront plus aléatoires.

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

besoin d'ajouter -fopenmpà l'indicateur du compilateur

Edit: 2 Comme suggéré par user71404, j’ai modifié les fonctions sumOnes et sumArray et c’est maintenant très rapide.

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

Avec openmp, c'est plus lent, car les atomics ajoutent trop de temps système.

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

Sans atomics, c'est encore plus rapide, mais je me trompe de résultat.

2137992 1147218 619297 321243 155815 70946 32919 15579

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

Pour comprendre sumArray, considérons que 16 bits représentent et un tableau de 8 nombres.
00 n'a pas 1 et représente -1
01 et 10 ont un 1 et représentent 0
11 ont deux 1 et représentent 1
Pour que le nombre de bits intégré soit compté à 1 [ http://en.wikipedia.org/wiki/ Hamming_weight] et à chaque groupe, nous supprimons 1. Cool.

sumOnes est juste de la magie noire.

Voici les derniers drapeaux de compilation et le code.

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 plus rapide

Comme le montre user20768, un portage simple du code vers Julia est environ deux fois plus rapide que PyPy. Mais nous pouvons faire beaucoup mieux que cela.

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

Vous pouvez exécuter ceci en utilisant julia -p 8 -e 'require("golf.jl");main()'(le 8 étant le nombre de processus, vous pouvez jouer avec). Sur la dernière version préliminaire de Julia, cela prend 0,7s contre 1m22s pour PyPy.

Si vous avez suffisamment de cœurs sur votre ordinateur et que vous pouvez éventuellement créer quelques instances AWS, vous devriez pouvoir en réduire encore plus :)

C, 1.210s

Avec le code de l'OP qui exécute 1m45.729s sur ma machine.

Compilation:

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

Un merci tout spécial: @dyp pour les drapeaux de compilation et les idées d'optimisation.

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

Exemple de sortie:

6334411 2527506 1042239 439328 192914 87005 40847 19728

Perl

C'est loin d'être aussi rapide que la solution C, mais c'est assez rapide pour un langage interprété de haut niveau, je pense. Cela réduit d'environ 40% le temps d'exécution de l'implémentation 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;

Algorithm :: Combinatorics est disponible dans Ubuntu ( sudo apt-get install libalgorithm-combinatorics-perl). Les autres modules utilisés sont des modules Perl de base, ils devraient donc déjà être installés dans le cadre de l'installation de base d'Ubuntu.

Julia: 4,66x plus lente!

Je commence vraiment à douter des statistiques sur leur site web ...

Notez que le code Julia suivant est en réalité une transcription directe du code Python de l'OP, sans aucune optimisation. J'utilise cette time()fonction pour exclure le temps de démarrage lent 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

Code de l'OP dans PyPy: 1 m 11,506 s

Julia sortie:

6332170
2525472
1041522
438761
193119
86873
40705
19662

RPython 0.187s (258x plus rapide)

Source originale avec PyPy2.2.1: 1m 6.718s

Maintenant, avec le threading, le support de base pour Python standard a été supprimé. Le nombre de threads de travail peut être spécifié en tant que paramètre de ligne de commande. La valeur par défaut est deux.

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 est un sous-ensemble restreint de Python, qui peut être traduit en C puis compilé à l'aide de RPython Toolchain . Son objectif est d'aider à la création d'interprètes de langue, mais il peut également être utilisé pour compiler des programmes simples comme celui ci-dessus. La plupart des fonctionnalités "sophistiquées" de Python, telles que itertoolsou même mapne sont pas disponibles.

Pour compiler, créez un clone local du référentiel pypy actuel et exécutez les opérations suivantes:

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

L'exécutable résultant sera nommé convolution-cou similaire dans le répertoire de travail actuel.

J'ai paramétré les variables d'entrée, le programme doit donc être exécuté comme suit:

convolution-c 8 8 1000

pour correspondre à l'exemple de code.

Notes d'implémentation

S in itertools.product([-1,1], repeat = n+m-1)devient S in xrange(1<<n+m-1), interprétant Scomme un bitmap: [ 0, 1] → [ -1, 1]

De la même manière, Fest également une carte de bits, avec chaque fois deux bits représentant une valeur unique:
[ 00, 01, 10, 11] → [ -1, 0, 0, 1]

Une table de vérité est utilisée pour rechercher le produit plutôt que pour effectuer une réplication multiple.

Etant donné que les entiers signés 32 bits sont utilisés, ils nne doivent pas dépasser 15 et n+m31. Un support arbitraire entier peut être obtenu avec le rpython.rlib.rbigintmodule, si nécessaire.

La première itération de la boucle de produit scalaire est déroulée et associée au test de nullité de F.

Un PRNG homebrew est utilisé, la liste des sources. L'auteur de l'article présente une période de 2 ³² -1 et affirme qu'il réussit tous les tests de Diehard sauf un, bien que je ne l'aie pas personnellement confirmé.

La graine aléatoire change toutes les millisecondes, ce qui est aussi bon d'utiliser un horodatage. De plus, chaque thread de xortravail identifie son identifiant de processus avec cette valeur, afin de s'assurer qu'ils ont chacun une graine différente.

Exemples de timings

2 fils de travail:

$ 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 fils de travail:

$ 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 fils de travail:

$ 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

Source d'origine de l'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

Calendrier pour 100 000 itérations:

$ 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 plus rapide) (modification du code d'Arman)

Code Op dans PyPy: 3 min 1,4 s

Les deux sont effectués dans le REPL, sans compter le temps nécessaire pour charger les packages.

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 

Le code d'Arman pose quelques problèmes en le rendant très lent: il utilise beaucoup de fonctions anonymes et de fonctions supérieures inutilement. Pour tester si tout le vecteur F est nul, pourquoi ne pas simplement écrire tout (F == 0) au lieu de tout (x-> x == 0, F)? Il est plus court et littéralement mille fois plus rapide.

Il utilise également sum (map (*, x, y)) comme produit scalaire au lieu de simplement dot (x, y). La première version 650 fois plus lente pour un vecteur de 10k double. Et la fonction de produit scalaire est implémentée comme une boucle for dans Julia pure.

En outre, la compréhension des tableaux est lente. Il vaut mieux écrire [0,1,0, -1] [rand (1: 4, n)] au lieu de [[-1 -1 0 0 1] [rand (1: 4)] pour j = 1: n] .

Enfin, les variables globales sont mauvaises pour Julia. Julia n'est rapide que si vous écrivez du code de manière à ce que le JIT et l'inférence de type fonctionnent. Une grande partie de ceci est la stabilité du type: Le compilateur doit être capable de s'assurer que le type d'une variable ne changera pas pendant qu'il est dans une boucle, par exemple.

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

Exemple de sortie:

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

Nimrod compile en C, donc le choix du compilateur C pour les questions d’arrière-plan aussi.

En utilisant clang, compilez avec:

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

En utilisant gcc, compilez avec:

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

Omettez --passc:-fltosi vous avez un ancien compilateur C qui ne prend pas en charge LTO. Omettez l' --cc=...option si le choix par défaut du compilateur C vous convient. Le code nécessite Nimrod 0.9.4 ou 0.9.5 .

Sur mon quadcore iMac (2,65 GHz core i5), le code tourne en environ 0,15 seconde avec gcc 4,9, 0,16 seconde avec clang, comparé à 88 secondes pour PyPy 2.2.1 (soit une accélération de plus de 500 fois). Malheureusement, je n'ai pas accès à une machine avec plus de quatre cœurs sur laquelle PyPy est également installé ou sur laquelle je pourrais facilement installer PyPy, bien que je dispose d'environ 0,1 seconde (avec beaucoup de bruit de mesure) sur un processeur AMD 64 cœurs Opteron 6376 1,4 GHz (selon / proc / cpuinfo) avec gcc 4.4.6.

L'implémentation essaie d'être fidèle à l'original plutôt que d'optimiser le code au détriment de la lisibilité, sans renoncer à des optimisations évidentes. Chose intéressante, la récursion de la queue initVecRand()est un peu plus rapide qu'une boucle avec une instruction break avec à la fois gcc et clang. Le déroulement manuel d'une itération de la convolveboucle de test à l'intérieur de la boucle principale a également entraîné une accélération, probablement due à une meilleure prédiction de branche.

Java

J'ai traduit la solution C ++ ci-dessus en 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));
  }
}

Sur ma machine, je reçois la sortie suivante pour le programme 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

Le programme OPs dure environ 53 secondes sur ma machine:

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

Le programme c ++ n'a exécuté que 0,15 seconde environ:

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

C'est environ 2,5 fois plus rapide que la solution Java correspondante (je n'ai pas exclu le démarrage de la machine virtuelle). Cette solution Java est environ 142 fois plus rapide que le programme exécuté avec PyPy.

Étant donné que j'étais personnellement intéressé, je suis iterspassé à 100_000 pour Java et C ++, mais le facteur de 2,5 n'a pas diminué en faveur de Java si quelque chose devenait plus grand.

EDIT: J'ai exécuté les programmes sur un PC Linux Linux 64 bits.

EDIT2: Je veux ajouter que j'ai commencé par une traduction brute du code 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));
    }
}

Ce programme a duré environ 3,6 secondes:

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

Ce qui est environ 14 fois plus rapide que la solution PyPy. (Le choix de la fonction aléatoire standard sur la fonction fastRandom entraîne un temps d'exécution de 5 secondes)

Python 3.5 + numpy 1.10.1, 3,76 secondes

Les tests ont été exécutés sur mon Macbook Pro. Le code d'OP a pris environ 6 minutes sur la même machine.

La raison pour laquelle je réponds à cette question est en fait parce que je n'ai pas 10 réputations et que je ne peux pas répondre à la partie I :-p

Ces derniers jours, j'ai essayé de comprendre comment effectuer des convolutions massives de manière efficace avec numpy (sans faire appel à un package tiers, ni même scipy). Lorsque je suis tombé sur cette série de défis au cours de mes recherches, j'ai décidé de l'essayer. Je suis peut-être arrivé trop tard à ce jeu, mais voici ma tentative d'utilisation de Python 3.5 et de 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)

J'ai pré-calculé les tableaux S et F, et aplati le tableau S tout en effectuant la convolution, ce qui (selon mes expériences) pourrait tirer parti de la vitesse de np.convolve. En d'autres termes, n'ayant pas trouvé de routine de convolution vectorisée, j'ai simulé le code en aplatissant le tableau entier et espéré que np.convolved ferait la vectorisation sous le capot pour moi, ce qui semblait fonctionner. Notez que j'ai utilisé mode = 'same' et coupé les éléments de début et de fin qui étaient inutiles.

Sur mon Macbook Pro, les résultats du test donnent 3,76 secondes . Quand j'ai exécuté le code d'OP (modifié en Python 3.5), j'ai eu environ 6 minutes . L'accélération est environ 100 fois.

Un inconvénient est que, comme les tableaux S et F doivent être stockés, la mémoire requise peut poser problème si les tailles sont trop grandes.

J'ai utilisé la même méthode pour la partie I et j'ai eu une accélération d'environ 60-100x sur mon ordinateur portable.

Comme je faisais tout sur mon Macbook Pro, si quelqu'un pouvait tester mon code et me faire savoir comment ça se passe sur votre machine, je l'apprécierais beaucoup!

J, ^130x ~ 50x accélération?

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

Fois sur un debian aléatoire:

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

Je pense qu'il y a place à l'amélioration.

C ++: x200 (i7 à 4 cœurs, devrait passer à x400 sur 8 cœurs)

Essayer une solution plus simple C ++ 11 (testée avec VS 2012, gcc et clang) avec parallélisation.

Pour que cela soit compilé et exécuté sous Linux avec gcc 4.8.1:

g ++ -O3 -msse -msse2 -msse3 -march = natif -std = c ++ 11 -pthread -Wl, - pas de golf requis.cpp

Sous Linux, nous devons également std::launch::asyncforcer plusieurs threads. Cela me manquait dans une version antérieure.

Dans Visual Studio (2012+), cela devrait fonctionner mais créer une version de compilation pour le timing ...

Sur mon i3 dual core à l'ancienne, cela s'exécute en ~ 0,9 seconde. Sur mon quad core i7, il est 0.319s vs Pypy 66 secondes.

Sur un i7 à 8 cœurs, cela devrait être dans la plage d'accélération x400. Le passage à des tableaux de style C pourrait accélérer les choses, mais je souhaitais rester avec des conteneurs C ++. Pour moi, il est intéressant de voir l'accélération que vous pouvez obtenir tout en restant relativement proche du domaine du problème et à un niveau relativement élevé, quelque chose que je pense que le C ++ est vraiment bon. Il convient également de noter la relative facilité de mise en parallèle utilisant des constructions C ++ 11.

La solution de bit de @ ilmale est très cool et fonctionne pour -1/1/0. On pourrait aussi jeter le SSE sur cela et peut-être obtenir une accélération significative.

Au-delà de la parallélisation, il existe un autre "truc" qui consiste à réduire le nombre de sommations. Exemple de résultats: 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

Fortran: Je l’ai jusqu’à 106x 155x 160x 316x d’ accélération lorsqu’on utilise un RNG Xorshift et OpenMP sur un processeur i7 à 4 cœurs. Autre que cela, il n'y a pas de gros trucs. Pour que l'itérateur construise S, j'utilise simplement la représentation binaire de l'entier 16 bits i. Vous remarquerez que, mis à part le générateur de ressources en ligne et le "itérateur" / mappage de i à S, le code est aussi puissant que le code Python.

Edit: suppression du "if" dans le Xorshift, en utilisant maintenant "r = abs (w / ...)" au lieu de "r = w / ...". Va de 106x à 155x.

Edit2: Ceci génère 15x autant de nombres aléatoires que la solution C ++. Si quelqu'un a une solution zéro-overhead pour convertir un entier aléatoire en un tableau de 0 et de 1 en Fortran, je suis tout ouïe. Ensuite, nous pourrions battre C ++ :)

Edit3: La première édition a introduit un bogue, comme l'a souligné Lembik. Ceci est corrigé maintenant, avec une amélioration minime de l'accélération. J'essaierai d'utiliser la suggestion d'Eelvex pour obtenir plus de rapidité.

Edit4: Profiling a indiqué que la conversion en réel et en entier avec nint () était lente. J'ai remplacé ceci par une division entière effectuant à la fois la mise à l'échelle et l'arrondi, passant de 160x à une accélération de 316x.

Compiler avec:

gfortran -O3 -march = natif -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

Exemple de sortie:

$ time ./a.out
6329624 2524831 1039787 438809 193044 6860 40486 19517
./a.out 1.45s utilisateur 0.00s système 746% cpu 0,192 au total

Code de l'OP:

$ time pypy golf.py
pypy golf.py 60.68s utilisateur 0.04s système 99% cpu 1: 00.74 total