C ++ nieco magiczna
0,84 ms z prostym RNG, 1,67 ms z c ++ 11 std :: knuth
0,16 ms z niewielką modyfikacją algorytmu (patrz edycja poniżej)
Implementacja Pythona działa na moim urządzeniu w 7,97 sekundy. Jest to więc 9488 do 4772 razy szybszy, w zależności od wybranego RNG.
#include <iostream>
#include <bitset>
#include <random>
#include <chrono>
#include <stdint.h>
#include <cassert>
#include <tuple>
#if 0
// C++11 random
std::random_device rd;
std::knuth_b gen(rd());
uint32_t genRandom()
{
return gen();
}
#else
// bad, fast, random.
uint32_t genRandom()
{
static uint32_t seed = std::random_device()();
auto oldSeed = seed;
seed = seed*1664525UL + 1013904223UL; // numerical recipes, 32 bit
return oldSeed;
}
#endif
#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
std::pair<unsigned, unsigned> convolve()
{
const uint32_t n = 6;
const uint32_t iters = 1000;
unsigned firstZero = 0;
unsigned bothZero = 0;
uint32_t S = (1 << (n+1));
// generate all possible N+1 bit strings
// 1 = +1
// 0 = -1
while ( S-- )
{
uint32_t s1 = S % ( 1 << n );
uint32_t s2 = (S >> 1) % ( 1 << n );
uint32_t fmask = (1 << n) -1; fmask |= fmask << 16;
static_assert( n < 16, "packing of F fails when n > 16.");
for( unsigned i = 0; i < iters; i++ )
{
// generate random bit mess
uint32_t F;
do {
F = genRandom() & fmask;
} 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 );
// calculate which bits in the expression S * F evaluate to +1
unsigned firstPosBits = ((s1 & posBits) | (~s1 & negBits));
// idem for -1
unsigned firstNegBits = ((~s1 & posBits) | (s1 & negBits));
if ( popcnt( firstPosBits ) == popcnt( firstNegBits ) )
{
firstZero++;
unsigned secondPosBits = ((s2 & posBits) | (~s2 & negBits));
unsigned secondNegBits = ((~s2 & posBits) | (s2 & negBits));
if ( popcnt( secondPosBits ) == popcnt( secondNegBits ) )
{
bothZero++;
}
}
}
}
return std::make_pair(firstZero, bothZero);
}
int main()
{
typedef std::chrono::high_resolution_clock clock;
int rounds = 1000;
std::vector< std::pair<unsigned, unsigned> > out(rounds);
// do 100 rounds to get the cpu up to speed..
for( int i = 0; i < 10000; i++ )
{
convolve();
}
auto start = clock::now();
for( int i = 0; i < rounds; i++ )
{
out[i] = convolve();
}
auto end = clock::now();
double seconds = std::chrono::duration_cast< std::chrono::microseconds >( end - start ).count() / 1000000.0;
#if 0
for( auto pair : out )
std::cout << pair.first << ", " << pair.second << std::endl;
#endif
std::cout << seconds/rounds*1000 << " msec/round" << std::endl;
return 0;
}
Skompiluj w wersji 64-bitowej dla dodatkowych rejestrów. Podczas korzystania z prostego generatora losowego pętle w trybie splotowym () działają bez dostępu do pamięci, wszystkie zmienne są przechowywane w rejestrach.
Jak to działa: zamiast przechowywać Si Fjako tablice w pamięci, jest zapisywany jako bity w uint32_t.
W przypadku S, gdy nstosowane są najmniej znaczące bity, gdy ustawiony bit oznacza +1 i wyłączony nieco oznacza -1.
Fwymaga co najmniej 2 bitów, aby utworzyć rozkład [-1, 0, 0, 1]. Odbywa się to poprzez generowanie losowych bitów i badanie 16 najmniej znaczących (nazywanych r) i 16 najbardziej znaczących bitów (nazywanych l). Jeśli l & ~rzałożymy, że F wynosi +1, jeśli ~l & rzałożymy, że Fwynosi -1. W przeciwnym razie Fwynosi 0. To generuje rozkład, którego szukamy.
Teraz mamy S, posBitsz ustawionym bitem w każdej lokalizacji, w której F == 1 i negBitsz ustawionym bitem w każdym miejscu, w którym F == -1.
Możemy udowodnić, że F * S(gdzie * oznacza mnożenie) pod warunkiem warunkuje do +1 (S & posBits) | (~S & negBits). Możemy również wygenerować podobną logikę dla wszystkich przypadków, w których F * Swartość wynosi -1. I na koniec wiemy, że sum(F * S)ocenia się na 0 wtedy i tylko wtedy, gdy w wyniku jest równa liczba -1 i + 1. Jest to bardzo łatwe do obliczenia, po prostu porównując liczbę bitów +1 i -1 bitów.
Ta implementacja używa 32-bitowych liczb całkowitych, a maksymalna nakceptowana liczba to 16. Możliwe jest skalowanie implementacji do 31 bitów poprzez modyfikację kodu generowania losowego i do 63 bitów przy użyciu uint64_t zamiast uint32_t.
edytować
Następująca funkcja zwojowa:
std::pair<unsigned, unsigned> convolve()
{
const uint32_t n = 6;
const uint32_t iters = 1000;
unsigned firstZero = 0;
unsigned bothZero = 0;
uint32_t fmask = (1 << n) -1; fmask |= fmask << 16;
static_assert( n < 16, "packing of F fails when n > 16.");
for( unsigned i = 0; i < iters; i++ )
{
// generate random bit mess
uint32_t F;
do {
F = genRandom() & fmask;
} 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+1));
// generate all possible N+1 bit strings
// 1 = +1
// 0 = -1
while ( S-- )
{
// calculate which bits in the expression S * F evaluate to +1
auto firstBits = (S & mask) ^ adjF;
auto secondBits = (S & ( mask << 1 ) ) ^ ( adjF << 1 );
bool a = desiredBits == popcnt( firstBits );
bool b = desiredBits == popcnt( secondBits );
firstZero += a;
bothZero += a & b;
}
}
return std::make_pair(firstZero, bothZero);
}
skraca czas działania do 0.160-0.161ms. Ręczne rozwijanie pętli (nie pokazano powyżej) powoduje, że 0.150. Mniej banalne n = 10, iter = 100000 przypadków działa poniżej 250ms. Jestem pewien, że uda mi się go uzyskać poniżej 50 ms, wykorzystując dodatkowe rdzenie, ale to zbyt łatwe.
Odbywa się to poprzez zwolnienie gałęzi wewnętrznej pętli i zamianę pętli F i S.
Jeśli bothZeronie jest to wymagane, mogę skrócić czas pracy do 0,02 ms, rzadko zapętlając wszystkie możliwe tablice S.