To kontynuacja tego, jak bardzo wolno jest Python? (Lub jak szybki jest twój język?) .

Okazuje się, że dla mojego ostatniego pytania trochę za łatwo było uzyskać przyspieszenie x100. Dla tych, którzy lubili wyzwania, ale chcą czegoś trudniejszego, w którym mogliby naprawdę wykorzystać swoje umiejętności niskiego poziomu, oto część II. Wyzwanie polega na uzyskaniu przyspieszenia x100 dla następującego kodu python przetestowanego na moim komputerze.

Aby utrudnić, tym razem używam pypy. Obecny czas dla mnie to 1 minuta i 7 sekund przy użyciu pypy 2.2.1.

Zasady

1. Pierwsza osoba, która prześle kod, który mogę uruchomić, jest poprawna i jest 100 razy szybsza na moim komputerze, otrzyma nagrodę w wysokości 50 punktów.
2. Zwycięstwo przyznam najszybszemu kodowi po tygodniu.

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

Dane wyjściowe powinny być podobne do

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

Musisz użyć losowego ziarna w kodzie, a dowolny generator liczb losowych, który jest wystarczająco dobry, aby dać odpowiedzi zbliżone do powyższego, zostanie zaakceptowany.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod.

Objaśnienie kodu

Ten kod iteruje wszystkie tablice S o długości n + m-1, które składają się na -1s i 1s. Dla każdej tablicy S próbkuje 1000 niezerowych losowych tablic F o długości n składającej się z -1,0 lub 1 z prawdopodobieństwem 1/4, 1/2, / 14 pobrania każdej wartości. Następnie oblicza iloczyn wewnętrzny między F i każdym oknem S o długości n, aż znajdzie niezerowy iloczyn wewnętrzny. Dodaje 1 do leadingzerocountskażdej pozycji, w której znalazł zero wewnętrznego produktu.

Status

• Perl . 2,7 razy spowolnienie przez @tobyink. (W porównaniu do pypy nie cpython.)

• J . 39-krotne przyspieszenie o @Eelvex.

• C . 59 razy przyspieszyć @ace.
• Julia . 197 razy szybszy, nie wliczając czasu rozruchu o jeszcze jedną minutę. 8,5 razy szybsze, w tym czas rozruchu (w tym przypadku jest szybszy przy użyciu 4 procesorów niż 8).
• Fortran . 438 razy przyspieszenie o @ semi-extrinsic.
• Rpython . 258 razy przyspieszyć @primo.
• C ++ . 508 razy przyspieszenie @ilmale.

(Przestałem mierzyć czas nowych ulepszeń, ponieważ są one zbyt szybkie, a iteracja była za mała).

Zwrócono uwagę, że czasy poniżej jednej sekundy są niewiarygodne, a także niektóre języki mają koszty początkowe. Argument jest taki, że jeśli chcesz dołączyć, powinieneś również uwzględnić czas kompilacji C / C ++ itp. Oto czasy najszybszego kodu z liczbą iteracji zwiększoną do 100 000.

• Julia . 42 sekundy do @ jeszcze minuty.
• C ++ . 14 sekund autorstwa @GuySirton.
• Fortran . 14s @ semi-extrinsic.
• C ++ . 12s przez @ilmale.
• Rpython . 18s @primo.
• C ++ . 5s @Stefan.

Zwycięzcą jest ... Stefan!

Opublikowane wyzwanie uzupełniające. Jak wysoko możesz iść? (Wyzwanie kodowania + algorytmy) . Ten jest trudniejszy.

C ++ nieco magiczna

~ 16 ms wielowątkowy, 56 ms jednokierunkowy. ~ 4000 przyspieszenia.

(przyspieszenie oparte jest na kodzie wielowątkowym na moim i7-2820QM i 1 min 9 sekundach wymienionych w pytaniu. Ponieważ system OP ma gorszą wydajność pojedynczego wątku niż mój procesor, ale lepszą wydajność wielowątkową, spodziewam się, że ta liczba będzie dokładna)

Część wielowątkowa jest dość nieefektywna z powodu pojawiania się wątków. Prawdopodobnie mógłbym zrobić to lepiej, korzystając z mojej niestandardowej biblioteki zadań, ale ta zawiera błędy w systemach unixowych. Aby uzyskać wyjaśnienie i prawie identyczny kod bez wątków, odwiedź https://codegolf.stackexchange.com/a/26485/20965 .

edytować

Każdemu wątkowi nadałem swój własny RNG i skróciłem długość bitu do 32, co skróciło czas działania o kilka 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;
}

Przykładowe dane wyjściowe:

   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 x 450 x 530

Zamiast tablicy użyłem bitów (i czarnej magii).
Dzięki @ace za szybszą losową funkcję.

Jak to działa: pierwsze 15 bitów liczby całkowitej sreprezentuje tablicę S[15]; zera oznaczają -1, te zera +1. Tablica Fjest budowana w podobny sposób. Ale z dwoma bitami dla każdego symbolu.

• 00 oznacza -1
• 01 i 10 oznaczają 0
• 11 reprezentuje 1

Przyczyna Si Fmają inną reprezentację muszę przeplatają Ssię ze sobą, aby były porównywalne z F.

• 0 (-1) stało się 00 (-1 w reprezentacji F)
• 1 (+1) stał się 11 (+1 w reprezentacji F)

Teraz możemy po prostu użyć Carnota do obliczenia produktu wewnętrznego. Pamiętaj, że jedna zmienna może przyjmować tylko wartość 00 lub 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)

Dla mnie wygląda to nie na xor. :)

Podsumowując, to tylko gra polegająca na zmianie i maskowaniu, nic naprawdę skomplikowanego.

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

Oto przykładowy wynik:

6332350 2525218 1041716 438741 192917 87159 41023 19908 

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

Program został skompilowany z:

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

na Fedorze 20 z gcc 4.8.2 CPU to i7 8core.

Prawdopodobnie uda mi się uzyskać kilka ms poprawiających parametry kompilatora.

Chociaż jest to czas rozwiązania OP na moim komputerze:

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

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

Edytować:

Po prostu dodając openmp i zmieniając kolejność for, mam zysk x3, co prowadzi do poprawy wydajności x450 względem kodu OP. : D W tym przypadku leadingZerotablica musi być atomowa. Losowe globalne ... są losowe, będą bardziej losowe.

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

trzeba dodać -fopenmpdo flagi kompilatora

Edycja: 2 Jako sugestia użytkownika 71404 Zmieniłem funkcje sumOnes i sumArray i teraz jest niesamowicie szybki.

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

Z Openmp jest wolniejszy, ponieważ atomika dodaje zbyt dużo narzutu.

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

Bez atomów jest jeszcze szybciej, ale mam błędny wynik.

2137992 1147218 619297 321243 155815 70946 32919 15579

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

Aby zrozumieć sumArray, należy wziąć pod uwagę, że 16-bitowa reprezentacja i tablica 8 liczb.
00 nie ma 1 i reprezentuje -1
01, a 10 ma jeden 1 i reprezentuje 0
11 ma dwa 1s i reprezentuje 1
Tak więc wbudowane liczenie liczby bitów ustawione na 1 [ http://en.wikipedia.org/wiki/ Hamming_weight] i do każdej grupy usuwamy 1. Cool.

sumOnes to tylko czarna magia.

Oto najnowsze flagi i kod kompilacji.

gcc -std = c ++ 11 -mfpmath = sse -O3 -flto -march = natywny -funroll-loop -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,7 s, 120x szybciej

Jak pokazał user20768, prosty port kodu do Julii jest około dwa razy szybszy niż PyPy. Ale możemy zrobić o wiele więcej.

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

Możesz uruchomić to za pomocą julia -p 8 -e 'require("golf.jl");main()'(8 to liczba procesów, z którymi możesz się pobawić). W najnowszej wersji wstępnej Julii zajmuje to 0,7s vs. 1m22s dla PyPy.

Jeśli masz wystarczającą liczbę rdzeni na swoim komputerze i być może rozpalisz kilka instancji AWS, powinieneś być w stanie trochę się ogolić :)

C, 1,210s

Z kodem OP działającym na moim komputerze 1m45.729s.

Kompilacja:

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

Specjalne podziękowania: @dyp za flagi kompilacji i pomysły na optymalizacje.

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

Przykładowe dane wyjściowe:

6334411 2527506 1042239 439328 192914 87005 40847 19728

Perl

Nie jest to tak szybkie jak rozwiązanie C, ale myślę, że jest dość szybkie w przypadku języka interpretowanego wysokiego poziomu. Zmniejsza to o 40% czas działania implementacji Pythona.

#!/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;

Algorytm :: kombinatoryka jest dostępny w Ubuntu ( sudo apt-get install libalgorithm-combinatorics-perl). Pozostałe używane moduły to podstawowe moduły Perla, więc powinny być już zainstalowane jako część podstawowej instalacji Ubuntu.

Julia: 4,66x wolniej!

Naprawdę zaczynam wątpić w statystyki na ich stronie internetowej ...

Zauważ, że poniższy kod Julii jest faktycznie bezpośrednią transkrypcją kodu Pythona OP bez żadnych optymalizacji. Używam tej time()funkcji, aby wykluczyć wolny czas uruchamiania Julii ...

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

Kod OP w PyPy: 1 m 11.506 s

Wyjście Julii:

6332170
2525472
1041522
438761
193119
86873
40705
19662

RPython 0,187 s (258x szybciej)

Oryginalne źródło w / PyPy2.2.1: 1m 6.718s

Teraz z wątkami, porzucono wsparcie dla standardowego Pythona. Liczbę wątków roboczych można określić jako parametr wiersza poleceń, domyślnie jest to dwa.

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 jest ograniczonym podzbiorem Pythona, który można przetłumaczyć na C, a następnie skompilować za pomocą RPython Toolchain . Jego wyraźnym celem jest pomoc w tworzeniu tłumaczy językowych, ale można go również używać do kompilacji prostych programów, takich jak ten powyżej. Większość „bardziej wyrafinowanych” funkcji Pythona, takich jak itertoolslub nawet mapnie są dostępne.

Aby skompilować, utwórz lokalny klon bieżącego repozytorium pypy i uruchom następujące polecenie:

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

Wynikowy plik wykonywalny zostanie nazwany convolution-club podobny w bieżącym katalogu roboczym.

Sparametryzowałem zmienne wejściowe, więc program powinien zostać uruchomiony jako:

convolution-c 8 8 1000

aby dopasować przykładowy kod.

Uwagi dotyczące wdrażania

S in itertools.product([-1,1], repeat = n+m-1)staje się S in xrange(1<<n+m-1), interpretując Sjako mapę bitową: [ 0, 1] → [-1 , 1]

Podobnie Fjest mapa bitowa z każdych dwóch bitów reprezentujących pojedyncze wartości:
[ 00, 01, 10, 11] → [ -1, 0,0 , 1]

Tablica prawdy służy do wyszukiwania produktu, a nie do wykonywania wielu prób.

Ponieważ używane są 32-bitowe liczby całkowite ze znakiem, nie nmogą być większe niż 15 i n+mnie większe niż 31. Arbitralną obsługę liczb całkowitych można uzyskać za pomocą parametrurpython.rlib.rbigint razie potrzeby moduł .

Pierwsza iteracja pętli iloczyn iloczynu jest rozwijana i łączona z testem nieważności F .

Używany jest domowy PRNG, źródło wymienione. Autor artykułu pokazuje okres 2 ³² -1 i twierdzi, że pomyślnie przeszedł wszystkie testy Dieharda z wyjątkiem jednego, chociaż osobiście tego nie potwierdziłem.

Losowe ziarno zmienia się co milisekundę, co jest tak samo dobre, jak pozwala na to znacznik czasu. Ponadto każdy wątek roboczy ma xorswój identyfikator procesu o tej wartości, aby mieć pewność, że każdy ma inne ziarno.

Przykładowe czasy

2 wątki robocze:

$ 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 wątki robocze:

$ 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 wątków roboczych:

$ 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

Oryginalne źródło 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

Czas dla 100000 iteracji:

$ 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 szybciej) (modyfikacja kodu Armana)

Kod operacji w PyPy: 3 min 1,4 s

Oba wykonano w REPL, nie licząc czasu na załadowanie pakietów.

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 

Występują pewne problemy z kodem Armana, który powoduje, że jest on bardzo wolny: Używa wielu anonimowych funkcji i funkcji wyższego rzędu niepotrzebnie. Aby sprawdzić, czy cały wektor F jest równy zero, dlaczego po prostu nie napisać wszystkich (F == 0) zamiast wszystkich (x-> x == 0, F)? Jest krótszy i dosłownie tysiąc razy szybszy.

Używa także sumy (mapa (*, x, y)) jako iloczynu kropkowego zamiast po prostu kropki (x, y). Pierwsza wersja 650 razy wolniejsza dla wektora 10k podwaja się. A funkcja iloczynu jest zaimplementowana jako pętla for w czystej Julii.

Również interpretacje tablic są powolne. Lepiej jest napisać [0,1,0, -1] [rand (1: 4, n)] zamiast [[-1 0 0 1] [rand (1: 4)] dla j = 1: n] .

Wreszcie, zmienne globalne są złe juju w Julii. Julia jest szybka tylko wtedy, gdy piszesz kod w sposób umożliwiający działanie JIT i wnioskowania typu. Dużą część tego stanowi stabilność typu: kompilator musi mieć pewność, że typ zmiennej nie zmieni się na przykład w pętli.

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

Przykładowe dane wyjściowe:

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

Nimrod kompiluje się do C, dlatego wybór kompilatora C również ma znaczenie.

Używając clang, skompiluj z:

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

Używając gcc, skompiluj z:

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

Pomiń, --passc:-fltojeśli masz starszy kompilator C, który nie obsługuje LTO. Pomiń tę --cc=...opcję, jeśli nie masz nic przeciwko domyślnemu wyborowi kompilatora C. Kod wymaga Nimrod 0.9.4 lub 0.9.5 .

Na moim quadcore iMac (rdzeń i5 2,66 GHz) kod działa w około .15 sekund z gcc 4.9, .16 sekund z clang, w porównaniu do 88 sekund dla PyPy 2.2.1 (tj. Przyspieszenie ponad 500 razy). Niestety nie mam dostępu do maszyny z więcej niż czterema rdzeniami, która ma również zainstalowany PyPy lub w której mógłbym łatwo zainstalować PyPy, chociaż dostaję około .1 sekundy (z dużą ilością szumów pomiarowych) na 64-rdzeniowym AMD Opteron 6376 1,4 GHz (zgodnie z / proc / cpuinfo) z gcc 4.4.6.

Implementacja stara się być wierna oryginalnemu, a nie optymalizować kod kosztem czytelności, nie rezygnując z oczywistych optymalizacji. Co ciekawe, rekurencja w ogonie initVecRand()jest nieco szybsza niż pętla z instrukcją break z gcc i clang. Ręczne rozwijanie jednej iteracji convolvepętli testowej w pętli głównej również spowodowało przyspieszenie, prawdopodobnie z powodu lepszego przewidywania gałęzi.

Jawa

Przetłumaczyłem powyższe rozwiązanie C ++ na Javę:

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

Na moim komputerze otrzymuję następujące dane wyjściowe dla programu 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

Program OPs działa na moim komputerze około 53 sekund:

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

Program c ++ wykonał tylko około 0,15 sekundy:

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

To około 2,5 razy szybciej niż odpowiednie rozwiązanie Java (nie wykluczyłem uruchamiania VM). Te rozwiązania Java są około 142 razy szybsze niż program wykonywany za pomocą PyPy.

Ponieważ byłem osobiście zainteresowany, ustawiłem iters na 100_000 dla Java i C ++, ale współczynnik 2,5 nie spadł na korzyść Javy, jeśli coś się powiększyło.

EDYCJA: Uruchomiłem programy na 64-bitowym komputerze Arch Linux.

EDIT2: Chcę dodać, że zacząłem od zgrubnego tłumaczenia kodu 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));
    }
}

Ten program działał około 3,6 sekundy:

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

Co jest około 14 razy szybsze niż rozwiązanie PyPy. (Wybranie standardowej funkcji losowej zamiast funkcji fastRandom prowadzi do wykonania 5 sekund)

Python 3.5 + numpy 1.10.1, 3,76 sekundy

Testy zostały przeprowadzone na moim Macbooku Pro. Kod OP zajął ~ 6 minut na tej samej maszynie.

Odpowiadam na to pytanie, ponieważ nie mam 10 reputacji i nie mogę odpowiedzieć na część I :-p

W ciągu ostatnich kilku dni starałem się dowiedzieć, jak skutecznie wykonywać ogromne zwoje przy pomocy numpy (bez polegania na pakiecie innej firmy, nawet scipy). Kiedy natknąłem się na tę serię wyzwań podczas moich badań, postanowiłem spróbować. Być może przyjechałem do tej gry za późno, ale oto moja próba użycia Pythona 3.5 i 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)

Wstępnie obliczyłem tablice S i F i spłaszczyłem tablicę S podczas wykonywania splotu, który (na podstawie moich eksperymentów) mógłby skorzystać z prędkości np. Konwolucji. Innymi słowy, ponieważ nie znalazłem procedury wektoryzacji wektoryzacji, podrobiłem wektoryzację kodu przez spłaszczenie całej tablicy i miałem nadzieję, że np. Konwój zrobi dla mnie wektoryzację pod maską, która wydawała się działać. Uwaga: Użyłem mode = „same” i przyciąłem wiodące i końcowe elementy, które były bezużyteczne.

Na moim Macbooku Pro wyniki testu dają 3,76 sekundy . Kiedy uruchomiłem kod OP (zmodyfikowany do Python 3.5), mam około 6 minut . Przyspieszenie jest około 100 razy.

Jedną wadą jest to, że ponieważ tablice S i F mają być przechowywane, zapotrzebowanie na pamięć może stanowić problem, jeśli rozmiary są zbyt duże.

Użyłem tej samej metody w części I i uzyskałem przyspieszenie ~ 60-100x na moim laptopie.

Ponieważ zrobiłem wszystko na moim Macbooku Pro, jeśli ktoś mógłby przetestować mój kod i dać mi znać, jak działa na twoim komputerze, byłbym bardzo wdzięczny!

J, przyspieszenie ^130x ~ 50x?

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

Czasy na losowym debianie:

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

Myślę, że jest miejsce na poprawę.

C ++: x200 (4-rdzeniowy i7, powinien być skalowany do x400 na 8-rdzeniowym)

Próbowanie prostszego C ++ 11 (testowane z VS 2012, gcc i clang) z równoległością.

Aby to skompilować i uruchomić pod Linuksem z gcc 4.8.1:

g ++ -O3-msse-msse2-msse3 -march = natywny -std = c ++ 11 -pthread -Wl, - nie-potrzebny golf.cpp

W systemie Linux potrzebujemy również std::launch::async wymusić wiele wątków. Brakowało mi tego we wcześniejszej wersji.

W programie Visual Studio (2012+) powinno to po prostu działać, ale należy utworzyć kompilację wydania dla pomiaru czasu ...

Na moim starszym dwurdzeniowym i3 działa to w ~ 0,9 sekundy. Na moim czterordzeniowym rdzeniu i7 jest to 0.319s vs. pypy 66 sekund.

Na 8-rdzeniowym i7 powinno to być w zakresie przyspieszenia x400. Przejście na tablice w stylu C przyspieszy to, ale byłem zainteresowany pozostaniem z kontenerami C ++. Dla mnie ciekawie jest zobaczyć przyspieszenie, które można uzyskać, pozostając stosunkowo blisko problematycznej domeny i na stosunkowo wysokim poziomie, w czym myślę, że C ++ jest naprawdę dobry. Na uwagę zasługuje również względna łatwość sparametryzowania przy użyciu konstrukcji C ++ 11.

@ Rozwiązanie bitowe ilmale jest bardzo fajne i działa dla -1/1/0. Można również rzucić SSE i może uzyskać znaczne przyspieszenie.

Poza paralelizacją jest jeszcze jedna „sztuczka”, która zmniejsza liczbę sumowań. Przykładowe wyniki: 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

Dobra, Fortran: Mam to do 106x 155x 160x 316x gdy używam Xorshift RNG i OpenMP na 4-rdzeniowym procesorze i7. Poza tym nie ma wielkich sztuczek. Aby iterator skonstruował S, po prostu używam binarnej reprezentacji 16-bitowej liczby całkowitej i. Zauważysz, że oprócz wbudowanego RNG i „iteratora” / mapowania od i do S, kod jest tak samo wysoki, jak kod Pythona.

Edycja: usunięto „if” w Xorshift, teraz używając „r = abs (w / ...)” zamiast „r = w / ...”. Zwiększa od 106x do 155x.

Edycja2: Generuje 15 razy więcej liczb losowych niż rozwiązanie C ++. Jeśli ktoś ma zerowe narzutowe rozwiązanie do przekształcania losowej liczby int w tablicę zer i jedynek w Fortranie, to jestem cały w uszach. Wtedy moglibyśmy pokonać C ++ :)

Edycja3: Pierwsza edycja wprowadziła błąd, jak zauważył Lembik. Jest to teraz naprawione, z niewielką poprawą przyspieszenia. Spróbuję skorzystać z sugestii Eelvex, aby uzyskać większe przyspieszenie.

Edycja4: Profilowanie wskazało, że konwersja do wartości rzeczywistej i powrót do liczby całkowitej za pomocą nint () była powolna. Zamieniłem to na jedną liczbę całkowitą, która wykonuje zarówno skalowanie, jak i zaokrąglanie, przechodząc od przyspieszenia 160x do 316x.

Połącz z:

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

Przykładowe dane wyjściowe:

czas $ ./a.out
6329624 2524831 1039787 438809 193044 6860 40486 19517
./a.out 1,45s użytkownik 0,00s system 746% procesor 0,192 łącznie

Kod PO:

$ czas pypy golf.py
pypy golf.py 60,68s użytkownik 0,04s system 99% procesor 1: 00.74 łącznie