Quelle est la lenteur de Python (Partie II)?

52

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.

Communauté
la source
3
une explication de ce que le code est supposé réaliser serait bien, afin que nous puissions le réécrire et pas simplement le porter
Einacio
6
" La première personne à soumettre un code que je peux exécuter est correcte et x100 fois plus rapide sur ma machine gagne immédiatement et la compétition se ferme. " Quel est le but de la clôture de la compétition de cette manière? Pourquoi ne pas utiliser une date limite comme la plupart des autres, afin que nous puissions la voir réduite dans d’autres langues?
GrovesNL
5
@Einacio C'est une bonne idée. J'ai changé les règles et j'espère que cela ne les dérangera pas.
1
@Lembik J'ai amélioré ma version Fortran, la rendant deux fois plus rapide encore sur ma machine. Pourriez-vous chronométrer à nouveau? :)
semi-extrinsèque
1
@ Fait semi-extrinsèque.

Réponses:

13

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,
Stefan
la source
La sortie n'est pas correcte je pense, peut-être qu'il y a un bug? Comparez à ce qui est dans la question. En particulier, la dernière colonne devrait être un nombre proche de 19180.
@ Lembik je peux voir ce que tu veux dire. Je pense que la sortie aléatoire n’est pas assez aléatoire, ce qui crée parfois une sortie géniale. Avec le générateur aléatoire C ++ 11, cela fonctionne bien. Je corrigerai le code plus tard aujourd'hui.
Stefan
J'obtiens Stefan.cpp: 104: 101: erreur: 'memcpy' n'a pas été déclaré dans cette étendue memcpy (out.data () + (threadId * m), LeadingZeros.data (), sizeof (LeadingZeros [0]) * m )
Je pense que je dois inclure string.h. Essaye encore.
Stefan
Vous compilez ceci avec g ++ -O3 -std = c ++ 0x -pthread -Wl, - non-selon les besoins Stefan.cpp -o Stefan
16

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;
}
ilmale
la source
Maintenant, j'ai hâte de tester ça! Malheureusement, ce ne sera pas avant quelques heures.
1
Ce qui suit était dans une suggestion de modification, mais je pense que cela pourrait correspondre mieux à un commentaire. Vous pouvez remplacer votre sumOnes, sumArray par ce qui suit (semble me donner une vitesse 2x plus rapide que la version openmp). 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; }cela a été suggéré par @ user71404
ace_HongKongIndependence
@ user71404: profileur a dit que le programme avait passé tout son temps dans ces deux fonctions, mais j'étais vraiment fatigué hier, je pense à quelque chose de mieux que ça. Je vais essayer ce soir (UTC). Merci.
ilmale
Souhaitez-vous changer le deuxième extrait de code pour qu'il soit la copie complète et le code pouvant être collé? Je dois faire quelque chose de mal dans mes tentatives pour obtenir votre code openmp afin que cela aide beaucoup.
Agréable. Je pensais que cela pourrait être fait avec des opérations de bits.
Guy Sirton
10

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

une minute de plus
la source
Je suis sûr que vous mesurez le mauvais timing. Python avec Pypy est également un langage basé sur JIT, mais les synchronisations effectuées par l'OP incluent le temps de compilation JIT. Vous l'excluez. J'ai installé la dernière version de Julia git et testé votre code. Sur ma machine, votre commande avec 8 processus prend 6,6 secondes, mais elle affiche le "temps écoulé 0,588 .. secondes".
semi-extrinsèque le
Le minutage Python inclut le démarrage de PyPy et le réchauffement de JIT, mais cela prend au plus quelques secondes. La différence d’une minute d’exécution est négligeable. Je suis heureux si le PO change le chronométrage de Python (cela ne changera rien), mais inclure le temps de démarrage de Julia ne serait pas raisonnable.
une minute de plus le
J'ai demandé au PO dans les commentaires à la question initiale, et il a déclaré: "Les horaires doivent tout inclure pour les langues de JIT." Il a également déclaré qu'il créerait un nouveau défi dans lequel les solutions prendraient beaucoup plus de temps qu'une seconde, laissant Julia dans la compétition.
semi-extrinsèque le
Dans ce cas, la solution optimale consiste à utiliser un algorithme série, qui prend environ 2 secondes. Je publierais le code, mais ce concours est maintenant quelque peu redondant, car tout le monde sait déjà que le C ++ démarre plus rapidement que tout le reste.
une minute de plus le
Je viens de publier ma solution Fortran, je ne vois donc pas pourquoi vous ne devriez pas publier la plus rapide des solutions Julia (si vous avez déjà le code).
semi-extrinsèque le
5

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
ace_HongKongIndependence
la source
1
Intéressant en effet, je peux faire des observations similaires en supprimant tous ces indicateurs d'optimisation. Essayez -march=native -fwhole-program -fstrict-aliasing -ftree-vectorizeBtw. Je suis descendu à moins de 4 secondes en utilisant du C ++ 11, dont un MT19937 plus un uniform_int_distribution.
jeudi
1
1.119s faisant une accélération d'environ 59!
1
@ace Oui, je voulais juste le signaler. C'était plus simple pour moi d'essayer quelques-uns des PRNG de bibliothèques standard en C ++. Notez que vous pouvez utiliser un résultat entier 32 bits d'un PRNG pour générer 8 entrées F.
jeudi
1
Puisque nest égal à 8, vous pouvez probablement utiliser AVX (ou 2 * SSE) pour calculer le produit dot avec un Sstockage approprié .
Michael M.
2
Version SSE, petite accélération: gist.github.com/anonymous/11394210 (n'oubliez pas d'inclure smmintrin.h)
Michael M.
5

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.

tobyink
la source
1
Cela n'affectera pas la vitesse, mais sa 0..N-1portée est la dernière map, non? Vous avez oublié de use warnings? :-) Bien que la logique dans OP soit déroutante, la fenêtre glissante n’atteint jamais le dernier élément de S.
user2846289
Ah Je pensais juste que les tableaux avaient des tailles incompatibles, alors j'ai désactivé, warningsce qui permettait de traiter les éléments manquants comme nuls. N-1améliore cela. Et cela améliore réellement très légèrement la vitesse - il est maintenant environ 40% plus rapide que l’implémentation Python.
tobyink le
Je pense que votre code nécessite une version très moderne de List :: Util. Sur Ubuntu 14.04 je reçois "tout" n’est pas exporté par le module List :: Util
Oh oui, c'est vrai - vous aurez probablement besoin d'installer List :: Util sur CPAN. anyOn peut également trouver List :: MoreUtils, qui, bien qu’il ne s’agisse pas d’un module de base, est l’un des modules CPAN les plus couramment utilisés.
tobyink
4

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
agar agile
la source
7
+1 pour votre <s> impudeur </ s>.
ace_HongKongIndependence
Les variables globales, les importations et la compréhension des tableaux sont lents. Ce n’est pas ainsi que l’on écrit en général un programme Julia pour la performance.
Alex A.
4

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 32 -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
primo
la source
Je n'ai jamais vu de programme rpython auparavant. C'est bien. Existe-t-il un programme équivalent en python pur que pypy peut exécuter en 1.03?
@ Lembik, j'aimerais bien en voir un. Je pensais que 4.7 était assez bon, considérant que ma première tentative de python pur était d'environ 15 ans.
primo
Oui, désolé pour le retard. Je n'ai pas encore le code en cours d'exécution, mais je le ferai dès que possible.
Vous devriez essayer d'ajouter un JIT. Maintenant que ce serait rapide!
Kirbyfan64sos
@ Lembik, merci pour la mention;) Par curiosité, a-t-il fonctionné le plus rapidement avec 4 threads de production, ou 8?
Primo
3

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.

utilisateur20768
la source
Merci! Je vois qu'il me reste encore beaucoup à apprendre sur le langage Julia avant de pouvoir faire des déclarations sur sa vitesse :) Vraiment heureux de constater que quelques corrections triviales apportées à mon code augmentent son temps d'exécution de plusieurs fois.
agar agar
2

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.

Reimer Behrends
la source
Comment obtenez-vous Nimrod pour Ubuntu?
@Lembik Une recherche rapide sur Google vous donnerait nimrod-lang.org/download.html
ace_HongKongIndependence
@ace J'avais aussi inclus le lien dans mon message (bien qu'il soit difficile de voir en bleu sur noir maintenant que je le regarde).
Reimer Behrends
Vous pourriez accélérer encore davantage en augmentant la taille de la graine à 128 bits: xorshift.di.unimi.it
user60561
2

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)

Dinfuehr
la source
2

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!

Essayer. Trop dur
la source
1

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.

Eelvex
la source
Le script Python est censé être exécuté en utilisant pypy, non python, ce qui explique pourquoi votre script semble accélérer 130 fois.
ace_HongKongIndependence
@ace oui j'ai remarqué mais je ne peux pas installer pypy: - / Je pense que l'ordre de grandeur restera.
Eelvex le
Pas nécessairement ... i.imgur.com/n566hzw.png
ace_HongKongIndependence
En effet, pas nécessairement.
Eelvex le
Quel problème avez-vous installer pypy?
1

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;
}
Guy Sirton
la source
1

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

semi-extrinsèque
la source
Ce que j’utilisais dans J était une liste de pré-construction de 4 ^ n nombres en base-4, puis convertie en triadique et en excluant 0. Le RNG choisit simplement dans cette liste.
Eelvex
Je ne suis pas sûr que votre code est correct. Pour 100 000 itérations, je reçois 633140285 271390368 118307997 52751245 23725837 10744292 4944464 2388125 mais je pense que cela devrait être plus proche de 633170604 252560981 104156146 43911426 19316309 8713324 4073378 19594. C'est un compromis
1
Ah, merci, @Lembik, mon édition d'accélération (suppression de l'instruction if) était en effet un bogue. J'ai mis à jour mon code, il devrait donc être correct maintenant. Je vais essayer de poster une version en utilisant la suggestion de Eelvex plus tard.
semi-extrinsèque
Cela a également accéléré il semble!
Oui, légère accélération je suppose. J'ai réalisé que j'ajoutais 1,0 et que je soustrayais 1,0 dans une boucle étroite.
semi-extrinsèque