У меня есть этот код, который я написал в Python / NumPy

from __future__ import division
import numpy as np
import itertools

n = 6
iters = 1000
firstzero = 0
bothzero = 0
""" The next line iterates over arrays of length n+1 which contain only -1s and 1s """
for S in itertools.product([-1, 1], repeat=n+1):
    """For i from 0 to iters -1 """
    for i in xrange(iters):
        """ Choose a random array of length n.
            Prob 1/4 of being -1, prob 1/4 of being 1 and prob 1/2 of being 0. """
        F = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n)
        """The next loop just makes sure that F is not all zeros."""
        while np.all(F == 0):
            F = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n)
        """np.convolve(F, S, 'valid') computes two inner products between
        F and the two successive windows of S of length n."""
        FS = np.convolve(F, S, 'valid')
        if FS[0] == 0:
            firstzero += 1
        if np.all(FS == 0):
            bothzero += 1

print("firstzero: %i" % firstzero)
print("bothzero: %i" % bothzero)

Он рассчитывает количество раз, когда свертка двух случайных массивов, один из которых длиннее другого, с определенным распределением вероятности, имеет 0 в первой позиции или 0 в обеих позициях.

Я поспорил с другом, который говорит, что Python - ужасный язык для написания кода, который должен быть быстрым. Это займет 9 секунд на моем компьютере. Он говорит, что это можно сделать в 100 раз быстрее, если написать на «правильном языке».

Задача состоит в том, чтобы узнать, действительно ли этот код может быть сделан в 100 раз быстрее на любом языке по вашему выбору. Я протестирую ваш код, и через неделю победит самый быстрый. Если кто-то опускается ниже 0,09, он автоматически выигрывает, а я проигрываю.

Статус

• Python . Алистер Баксон ускоряет в 30 раз! Хотя это и не самое быстрое решение, оно на самом деле мое любимое.
• Октав . 100 раз ускоряется @Thethos.
• Rust . В 500 раз ускоряется @dbaupp.
• С ++ . Гай Сиртон в 570 раз ускоряется.
• C . 727 раз ускоряется @ace.
• С ++ . Невероятно быстро @Stefan.

Самые быстрые решения сейчас слишком быстры, чтобы разумно рассчитывать время. Поэтому я увеличил n до 10 и установил iters = 100000, чтобы сравнить лучшие. По этой мере самыми быстрыми являются.

• C . 7,5 с @ace.
• С ++ . 1s @Stefan.

Моя машина Время будет запущено на моей машине. Это стандартная установка Ubuntu на восьмиъядерный процессор AMD FX-8350. Это также означает, что мне нужно иметь возможность запускать ваш код.

Последующие сообщения Поскольку это соревнование было слишком легким, чтобы получить ускорение x100, я опубликовал продолжение для тех, кто хочет использовать свои знания в области скоростных гуру. Смотрите, насколько медлен Python на самом деле (часть II)?

С ++ немного магии

0,84 мс с простым RNG, 1,67 мс с c ++ 11 std :: knuth

0,16 мс с небольшой алгоритмической модификацией (см. Редактирование ниже)

Реализация Python выполняется на моей установке за 7,97 секунды. Так что это в 9488–4772 раза быстрее, в зависимости от того, какую ГСЧ вы выберете.

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

Компилировать в 64-битные для дополнительных регистров. При использовании простого генератора случайных чисел циклы в convolve () запускаются без доступа к памяти, все переменные сохраняются в регистрах.

Как это работает: вместо того, чтобы хранить Sи Fкак массивы в памяти, он сохраняется как биты в uint32_t.
Например S, используются nмладшие значащие биты, где установленный бит обозначает +1, а неустановленный бит обозначает -1.
Fтребуется как минимум 2 бита для создания распределения [-1, 0, 0, 1]. Это делается путем генерации случайных битов и изучения 16 наименее значимых (вызываемых r) и 16 наиболее значимых (вызываемых l) битов . Если l & ~rмы предположим, что F равен +1, если ~l & rмы предположим, что Fэто -1. В противном случаеF это 0. Это создает распределение, которое мы ищем.

Теперь у нас есть S, posBitsс установленным битом в каждом месте, где F == 1 иnegBits с набором бит на каждом месте , где F == -1.

Мы можем доказать, что F * S(где * обозначает умножение) при условии 1 (S & posBits) | (~S & negBits). Мы также можем сгенерировать аналогичную логику для всех случаев, где значение F * Sравно -1. И, наконец, мы знаем, чтоsum(F * S) значение равно 0 тогда и только тогда, когда в результате есть равное количество -1 и +1. Это очень легко рассчитать, просто сравнив количество +1 бит и -1 бит.

Эта реализация использует 32-битные целые, а максимальная n допустимое значение равно 16. Можно масштабировать реализацию до 31 бита путем изменения кода случайной генерации и до 63 битов, используя uint64_t вместо uint32_t.

редактировать

Следующая сверточная функция:

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

сокращает время выполнения до 0,160-0,161мс. Ручное развертывание петли (не показано выше) составляет 0,150. Менее тривиальный n = 10, iter = 100000 кейс работает менее 250 мс. Я уверен, что смогу получить его за 50 мс, используя дополнительные ядра, но это слишком просто.

Это делается путем освобождения ветви внутреннего цикла и замены F и S цикла.
Если bothZeroне требуется, я могу сократить время выполнения до 0,02 мс, разреживая циклы по всем возможным S-массивам.

Python2.7 + Numpy 1.8.1: 10.242 с

Фортран 90+: 0,029 с 0,003 с 0,022 с 0,010 с

Черт побери, ты проиграл! Здесь тоже нет ни капли распараллеливания, просто прямой Fortran 90+.

РЕДАКТИРОВАТЬ Я взял алгоритм Гая Сиртона для перестановки массива S(хорошая находка: D). Очевидно, у меня также были -g -tracebackактивны флаги компилятора, которые замедляли этот код примерно до 0,017 с. В настоящее время я собираю это как

ifort -fast -o convolve convolve_random_arrays.f90

Для тех, у кого нет ifort, вы можете использовать

gfortran -O3 -ffast-math -o convolve convolve_random_arrays.f90

РЕДАКТИРОВАТЬ 2 : Уменьшение времени выполнения, потому что я делал что-то не так ранее и получил неправильный ответ. Делать это правильно, по-видимому, медленнее. Я до сих пор не могу поверить, что C ++ работает быстрее, чем мой, поэтому я, вероятно, собираюсь потратить некоторое время на этой неделе, пытаясь исправить это дерьмо, чтобы ускорить его.

РЕДАКТИРОВАТЬ 3 : Просто изменив секцию RNG, используя секцию, основанную на RNG BSD (как предложено Sampo Smolander) и исключив постоянное деление на m1, я сократил время выполнения до того же, что и ответ C ++ Гая Сиртона . Использование статических массивов (как предложено Sharpie) понижает время выполнения до уровня времени выполнения C ++! Yay Фортран! : D

РЕДАКТИРОВАТЬ 4 Очевидно, что это не компилируется (с gfortran) и работает правильно (неправильные значения), потому что целые числа выходят за свои пределы. Я внес исправления, чтобы убедиться, что он работает, но для этого требуется, чтобы у него был либо ifort 11+, либо gfortran 4.7+ (или другой компилятор, который позволяет iso_fortran_envи int64тип F2008 ).

Вот код:

program convolve_random_arrays
   use iso_fortran_env
   implicit none
   integer(int64), parameter :: a1 = 1103515245
   integer(int64), parameter :: c1 = 12345
   integer(int64), parameter :: m1 = 2147483648
   real, parameter ::    mi = 4.656612873e-10 ! 1/m1
   integer, parameter :: n = 6
   integer :: p, pmax, iters, i, nil(0:1), seed
   !integer, allocatable ::  F(:), S(:), FS(:)
   integer :: F(n), S(n+1), FS(2)

   !n = 6
   !allocate(F(n), S(n+1), FS(2))
   iters = 1000
   nil = 0

   !call init_random_seed()

   S = -1
   pmax = 2**(n+1)
   do p=1,pmax
      do i=1,iters
         F = rand_int_array(n)
         if(all(F==0)) then
            do while(all(F==0))
               F = rand_int_array(n)
            enddo
         endif

         FS = convolve(F,S)

         if(FS(1) == 0) then
            nil(0) = nil(0) + 1
            if(FS(2) == 0) nil(1) = nil(1) + 1
         endif

      enddo
      call permute(S)
   enddo

   print *,"first zero:",nil(0)
   print *," both zero:",nil(1)

 contains
   pure function convolve(x, h) result(y)
!x is the signal array
!h is the noise/impulse array
      integer, dimension(:), intent(in) :: x, h
      integer, dimension(abs(size(x)-size(h))+1) :: y
      integer:: i, j, r
      y(1) = dot_product(x,h(1:n-1))
      y(2) = dot_product(x,h(2:n  ))
   end function convolve

   pure subroutine permute(x)
      integer, intent(inout) :: x(:)
      integer :: i

      do i=1,size(x)
         if(x(i)==-1) then
            x(i) = 1
            return
         endif
         x(i) = -1
      enddo
   end subroutine permute

   function rand_int_array(i) result(x)
     integer, intent(in) :: i
     integer :: x(i), j
     real :: y
     do j=1,i
        y = bsd_rng()
        if(y <= 0.25) then
           x(j) = -1
        else if (y >= 0.75) then
           x(j) = +1
        else
           x(j) = 0
        endif
     enddo
   end function rand_int_array

   function bsd_rng() result(x)
      real :: x
      integer(int64) :: b=3141592653
      b = mod(a1*b + c1, m1)
      x = real(b)*mi
   end function bsd_rng
end program convolve_random_arrays

Я полагаю, что теперь вопрос заключается в том, перестанете ли вы использовать Python с медленной скоростью, как меласса, и использовать Fortran с быстрым движением электронов;).

Python 2,7 - 0,882 с 0,283 с

(Оригинал ОП: 6.404с)

Редактировать: оптимизация Стивена Румбальски путем предварительного вычисления значений F. С этой оптимизацией cpython превосходит 0,365 Pypy.

import itertools
import operator
import random

n=6
iters = 1000
firstzero = 0
bothzero = 0

choicesF = filter(any, itertools.product([-1, 0, 0, 1], repeat=n))

for S in itertools.product([-1,1], repeat = n+1):
    for i in xrange(iters):
        F = random.choice(choicesF)
        if not sum(map(operator.mul, F, S[:-1])):
            firstzero += 1
            if not sum(map(operator.mul, F, S[1:])):
                bothzero += 1

print "firstzero", firstzero
print "bothzero", bothzero

Оригинальный код OP использует такие крошечные массивы, что использование Numpy бесполезно, как демонстрирует эта чистая реализация на python. Но посмотрите также на эту клочковатую реализацию, которая в три раза быстрее моего кода.

Я также оптимизирую, пропуская остальную часть свертки, если первый результат не равен нулю.

Ржавчина: 0,011 с

Оригинальный Python: 8,3

Прямой перевод оригинального Python.

extern crate rand;

use rand::Rng;

static N: uint = 6;
static ITERS: uint = 1000;

fn convolve<T: Num>(into: &mut [T], a: &[T], b: &[T]) {
    // we want `a` to be the longest array
    if a.len() < b.len() {
        convolve(into, b, a);
        return
    }

    assert_eq!(into.len(), a.len() - b.len() + 1);

    for (n,place) in into.mut_iter().enumerate() {
        for (x, y) in a.slice_from(n).iter().zip(b.iter()) {
            *place = *place + *x * *y
        }
    }
}

fn main() {
    let mut first_zero = 0;
    let mut both_zero = 0;
    let mut rng = rand::XorShiftRng::new().unwrap();

    for s in PlusMinus::new() {
        for _ in range(0, ITERS) {
            let mut f = [0, .. N];
            while f.iter().all(|x| *x == 0) {
                for p in f.mut_iter() {
                    match rng.gen::<u32>() % 4 {
                        0 => *p = -1,
                        1 | 2 => *p = 0,
                        _ => *p = 1
                    }
                }
            }

            let mut fs = [0, .. 2];
            convolve(fs, s, f);

            if fs[0] == 0 { first_zero += 1 }
            if fs.iter().all(|&x| x == 0) { both_zero += 1 }
        }
    }

    println!("{}\n{}", first_zero, both_zero);
}



/// An iterator over [+-]1 arrays of the appropriate length
struct PlusMinus {
    done: bool,
    current: [i32, .. N + 1]
}
impl PlusMinus {
    fn new() -> PlusMinus {
        PlusMinus { done: false, current: [-1, .. N + 1] }
    }
}

impl Iterator<[i32, .. N + 1]> for PlusMinus {
    fn next(&mut self) -> Option<[i32, .. N+1]> {
        if self.done {
            return None
        }

        let ret = self.current;

        // a binary "adder", that just adds one to a bit vector (where
        // -1 is the zero, and 1 is the one).
        for (i, place) in self.current.mut_iter().enumerate() {
            *place = -*place;
            if *place == 1 {
                break
            } else if i == N {
                // we've wrapped, so we want to stop after this one
                self.done = true
            }
        }

        Some(ret)
    }
}

• Составлено с --opt-level=3
• Мой компилятор ржавчины - недавний ночной : ( rustc 0.11-pre-nightly (eea4909 2014-04-24 23:41:15 -0700)чтобы быть точным)

C ++ (VS 2012) - 0,026 с 0,015 с

Python 2.7.6 / Numpy 1.8.1 - 12 с

Ускорение ~ x800.

Разрыв был бы намного меньше, если бы свернутые массивы были очень большими ...

#include <vector>
#include <iostream>
#include <ctime>

using namespace std;

static unsigned int seed = 35;

int my_random()
{
   seed = seed*1664525UL + 1013904223UL; // numerical recipes, 32 bit

   switch((seed>>30) & 3)
   {
   case 0: return 0;
   case 1: return -1;
   case 2: return 1;
   case 3: return 0;
   }
   return 0;
}

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

void convolve(vector<int>& out, const vector<int>& v1, const vector<int>& v2)
{
   for(size_t i = 0; i<out.size(); ++i)
   {
      int result = 0;
      for(size_t j = 0; j<v2.size(); ++j)
      {
         result += v1[i+j]*v2[j];
      }
      out[i] = result;
   }
}

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

void convolve_random_arrays(void)
{
   const size_t n = 6;
   const int two_to_n_plus_one = 128;
   const int iters = 1000;
   int bothzero = 0;
   int firstzero = 0;

   vector<int> S(n+1);
   vector<int> F(n);
   vector<int> FS(2);

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

   for(auto &x : S)
   {
      x = -1;
   }
   for(int i=0; i<two_to_n_plus_one; ++i)
   {
      for(int j=0; j<iters; ++j)
      {
         do
         {
            for(auto &x : F)
            {
               x = my_random();
            }
         } while(allzero(F));
         convolve(FS, S, F);
         if(FS[0] == 0)
         {
            firstzero++;
            if(FS[1] == 0)
            {
               bothzero++;
            }
         }
      }
      advance(S);
   }
   cout << firstzero << endl; // This output can slow things down
   cout << bothzero << endl; // comment out for timing the algorithm
}

Несколько заметок:

• Случайная функция вызывается в цикле, поэтому я выбрал очень легкий линейный конгруэнтный генератор (но щедро посмотрел на MSB).
• Это действительно только отправная точка для оптимизированного решения.
• Не заняло так много времени, чтобы написать ...
• Я перебираю все значения S, принимая S[0]«наименее значимую» цифру.

Добавьте эту основную функцию для отдельного примера:

int main(int argc, char** argv)
{
  for(int i=0; i<1000; ++i) // run 1000 times for stop-watch
  {
      convolve_random_arrays();
  }
}

С

Занимает 0,015 с на моей машине, а исходный код OP занимает ~ 7,7 с. Попытка оптимизации путем генерации случайного массива и свертки в одном и том же цикле, но, похоже, не имеет большого значения.

Первый массив генерируется путем взятия целого числа, записи его в двоичном виде и изменения всех 1 на -1 и всех 0 на 1. Остальные должны быть очень простыми.

Редактировать: вместо того, чтобы иметь nкак int, теперь у нас nесть макроопределенная константа, поэтому мы можем использовать int arr[n];вместо malloc.

Edit2: вместо встроенной rand()функции теперь реализован PRNG ксоршифта. Кроме того, многие условные операторы удаляются при генерации случайного массива.

Составьте инструкции:

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

Код:

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

#define n (6)
#define iters (1000)
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 main() {
    int firstzero=0, bothzero=0;
    int arr[n+1];
    unsigned int i, j;
    x=(int)time(NULL);

    for(i=0; i< 1<<(n+1) ; i++) {
        unsigned int tmp=i;
        for(j=0; j<n+1; j++) {
            arr[j]=(tmp&1)*(-2)+1;
            tmp>>=1;
        }
        for(j=0; j<iters; j++) {
            int randArr[n];
            unsigned int k, flag=0;
            int first=0, second=0;
            do {
                for(k=0; k<n; k++) {
                    randArr[k]=(1-(myRand()&3))%2;
                    flag+=(randArr[k]&1);
                    first+=arr[k]*randArr[k];
                    second+=arr[k+1]*randArr[k];
                }
            } while(!flag);
            firstzero+=(!first);
            bothzero+=(!first&&!second);
        }
    }
    printf("firstzero %d\nbothzero %d\n", firstzero, bothzero);
    return 0;
}

J

Я не рассчитываю выбивать какие-либо скомпилированные языки, и что-то подсказывает мне, что потребуется чудесная машина, чтобы получить менее 0,09 с, но я все равно хотел бы представить этот J, потому что он довольно гладкий.

NB. constants
num =: 6
iters =: 1000

NB. convolve
NB. take the multiplication table                */
NB. then sum along the NE-SW diagonals           +//.
NB. and keep the longest ones                    #~ [: (= >./) #/.
NB. operate on rows of higher dimensional lists  " 1
conv =: (+//. #~ [: (= >./) #/.) @: (*/) " 1

NB. main program
S  =: > , { (num+1) # < _1 1                NB. all {-1,1}^(num+1)
F  =: (3&= - 0&=) (iters , num) ?@$ 4       NB. iters random arrays of length num
FS =: ,/ S conv/ F                          NB. make a convolution table
FB =: +/ ({. , *./)"1 ] 0 = FS              NB. first and both zero
('first zero ',:'both zero ') ,. ":"0 FB    NB. output results

Это займет около 0,5 с на ноутбуке предыдущего десятилетия, всего в 20 раз быстрее, чем Python в ответе. Большая часть времени уходит, convпотому что мы пишем это лениво (мы вычисляем всю свертку) и в полной общности.

Поскольку мы знаем что-то о Sи F, мы можем ускорить процесс, сделав определенные оптимизации для этой программы. Лучшее, что я смог придумать, - это conv =: ((num, num+1) { +//.)@:(*/)"1выбрать конкретно два числа, которые соответствуют от диагональных сумм самым длинным элементам свертки, что примерно вдвое меньше времени.

Perl - в 9,3 раза быстрее ... улучшение на 830%

На моем древнем нетбуке запуск кода OP занимает 53 секунды; Версия Алистера Бакстона занимает около 6,5 секунд, а следующая версия Perl - около 5,7 секунд.

use v5.10;
use strict;
use warnings;

use Algorithm::Combinatorics qw( variations_with_repetition );
use List::Util qw( any sum );
use List::MoreUtils qw( pairwise );

my $n         = 6;
my $iters     = 1000;
my $firstzero = 0;
my $bothzero  = 0;

my $variations = variations_with_repetition([-1, 1], $n+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;
    }

    # The pairwise function doesn't accept array slices,
    # so need to copy into a temp array @S0
    my @S0 = @$S[0..$n-1];

    unless (sum pairwise { $a * $b } @F, @S0)
    {
      $firstzero++;
      my @S1 = @$S[1..$n];  # copy again :-(
      $bothzero++ unless sum pairwise { $a * $b } @F, @S1;
    }
  }
}

say "firstzero ", $firstzero;
say "bothzero ", $bothzero;

Python 2.7 - numpy 1.8.1 с привязками mkl - 0.086 с

(Оригинал ОП: 6,404 с) (чистый питон Бакстона: 0,270 с)

import numpy as np
import itertools

n=6
iters = 1000

#Pack all of the Ses into a single array
S = np.array( list(itertools.product([-1,1], repeat=n+1)) )

# Create a whole array of test arrays, oversample a bit to ensure we 
# have at least (iters) of them
F = np.random.rand(int(iters*1.1),n)
F = ( F < 0.25 )*-1 + ( F > 0.75 )*1
goodrows = (np.abs(F).sum(1)!=0)
assert goodrows.sum() > iters, "Got very unlucky"
# get 1000 cases that aren't all zero
F = F[goodrows][:iters]

# Do the convolution explicitly for the two 
# slots, but on all of the Ses and Fes at the 
# same time
firstzeros = (F[:,None,:]*S[None,:,:-1]).sum(-1)==0
secondzeros = (F[:,None,:]*S[None,:,1:]).sum(-1)==0

firstzero_count = firstzeros.sum()
bothzero_count = (firstzeros * secondzeros).sum()
print "firstzero", firstzero_count
print "bothzero", bothzero_count

Как указывает Бакстон, в исходном коде OP используются такие крошечные массивы, что использование Numpy бесполезно. Эта реализация использует numpy, выполняя все случаи F и S одновременно с ориентацией на массив. Это в сочетании с привязками mkl для python приводит к очень быстрой реализации.

Также обратите внимание, что простая загрузка библиотек и запуск интерпретатора занимает 0,076 секунды, поэтому фактические вычисления занимают ~ 0,01 секунды, что аналогично решению C ++.

MATLAB 0.024s

Компьютер 1

• Исходный код: ~ 3,3 с
• Код Алистара Бакстона: ~ 0.51 с
• Новый код Алистара Бакстона: ~ 0,25 с
• Код Matlab: ~ 0,024 с (Matlab уже запущен)

Компьютер 2

• Исходный код: ~ 6,66 с
• Код Алистара Бакстона: ~ 0,64 с
• Новый код Алистара Бакстона:
• Matlab: ~ 0,07 с (Matlab уже работает)
• Октава: ~ 0,07 с

Я решил дать о, так медленно, Matlab попробовать. Если вы знаете, как, вы можете избавиться от большинства циклов (в Matlab), что делает его довольно быстрым. Однако требования к памяти выше, чем для зацикленных решений, но это не будет проблемой, если у вас нет очень больших массивов ...

function call_convolve_random_arrays
tic
convolve_random_arrays
toc
end

function convolve_random_arrays

n = 6;
iters = 1000;
firstzero = 0;
bothzero = 0;

rnd = [-1, 0, 0, 1];

S = -1 *ones(1, n + 1);

IDX1 = 1:n;
IDX2 = IDX1 + 1;

for i = 1:2^(n + 1)
    F = rnd(randi(4, [iters, n]));
    sel = ~any(F,2);
    while any(sel)
        F(sel, :) = rnd(randi(4, [sum(sel), n]));
        sel = ~any(F,2);
    end

    sum1 = F * S(IDX1)';
    sel = sum1 == 0;
    firstzero = firstzero + sum(sel);

    sum2 = F(sel, :) * S(IDX2)';
    sel = sum2 == 0;
    bothzero = bothzero + sum(sel);

    S = permute(S); 
end

fprintf('firstzero %i \nbothzero %i \n', firstzero, bothzero);

end

function x = permute(x)

for i=1:length(x)
    if(x(i)==-1)
        x(i) = 1;
            return
    end
        x(i) = -1;
end

end

Вот что я делаю:

• используйте функцию Kyle Kanos для перестановки S
• рассчитать все n * iters случайных чисел одновременно
• сопоставьте с 1 по 4 до [-1 0 0 1]
• использовать умножение матриц (поэлементная сумма (F * S (1: 5)) равно умножению матриц F * S (1: 5) '
• для обоих нулей: рассчитывать только те элементы, которые выполняют первое условие

Я предполагаю, что у вас нет matlab, что очень плохо, так как мне бы очень хотелось посмотреть, как он сравнивается ...

(Функция может быть медленнее при первом запуске.)

Юлия: 0,30 с

Op's Python: 21,36 с (Core2 Duo)

71-кратное ускорение

function countconv()                                                                                                                                                           
    n = 6                                                                                                                                                                      
    iters = 1000                                                                                                                                                               
    firstzero = 0                                                                                                                                                              
    bothzero = 0                                                                                                                                                               
    cprod= Iterators.product(fill([-1,1], n+1)...)                                                                                                                             
    F=Array(Float64,n);                                                                                                                                                        
    P=[-1. 0. 0. 1.]                                                                                                                                                                                                                                                                                                             

    for S in cprod                                                                                                                                                             
        Sm=[S...]                                                                                                                                                              
        for i = 1:iters                                                                                                                                                        
            F=P[rand(1:4,n)]                                                                                                                                                  
            while all(F==0)                                                                                                                                                   
                F=P[rand(1:4,n)]                                                                                                                                              
            end                                                                                                                                                               
            if  dot(reverse!(F),Sm[1:end-1]) == 0                                                                                                                           
                firstzero += 1                                                                                                                                                 
                if dot(F,Sm[2:end]) == 0                                                                                                                              
                    bothzero += 1                                                                                                                                              
                end                                                                                                                                                            
            end                                                                                                                                                                
        end                                                                                                                                                                    
    end
    return firstzero,bothzero
end

Я сделал несколько модификаций ответа Джулии от Армана: во-первых, я обернул его в функцию, так как глобальные переменные затрудняют вывод типа Джулии и JIT: глобальная переменная может изменить свой тип в любое время и должна проверяться при каждой операции , Затем я избавился от анонимных функций и массивов. Они на самом деле не нужны, и все еще довольно медленные. Джулия быстрее с низкоуровневыми абстракциями прямо сейчас.

Есть намного больше способов сделать это быстрее, но это делает достойную работу.

Хорошо, я публикую это только потому, что я чувствую, что Java должна быть представлена ​​здесь. Я ужасно разбираюсь с другими языками и признаюсь, что не совсем понял проблему, поэтому мне понадобится некоторая помощь, чтобы исправить этот код. Я украл большую часть примера кода туза C, а затем позаимствовал некоторые фрагменты у других. Я надеюсь, что это не подделка ...

Одна вещь, на которую я хотел бы обратить внимание, состоит в том, что языки, которые оптимизируются во время выполнения, должны запускаться несколько / много раз, чтобы достичь полной скорости. Я думаю, что вполне оправданно брать полностью оптимизированную скорость (или, по крайней мере, среднюю скорость), потому что большинство вещей, которые вам нужны для быстрой работы, будут выполняться несколько раз.

Код все еще нужно исправить, но я все равно запустил его, чтобы посмотреть, сколько раз я получу.

Вот результаты на процессоре Intel (R) Xeon (R) E3-1270 V2 @ 3,50 ГГц на Ubuntu, работающем 1000 раз:

сервер: / tmp # time java8 -cp. тестер

firstzero 40000

оба ноль 20000

Время первого запуска: 41 мс Время последнего запуска: 4 мс

реальный 0m5.014s пользователь 0m4.664s sys 0m0.268s

Вот мой дерьмовый код:

public class Tester 
{
    public static void main( String[] args )
    {
        long firstRunTime = 0;
        long lastRunTime = 0;
        String testResults = null;
        for( int i=0 ; i<1000 ; i++ )
        {
            long timer = System.currentTimeMillis();
            testResults = new Tester().runtest();
            lastRunTime = System.currentTimeMillis() - timer;
            if( i ==0 )
            {
                firstRunTime = lastRunTime;
            }
        }
        System.err.println( testResults );
        System.err.println( "first run time: " + firstRunTime + " ms" );
        System.err.println( "last run time: " + lastRunTime + " ms" );
    }

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

    public String runtest()
    {
        int n = 6;
        int iters = 1000;
        //#define iters (1000)
        //PRNG seeds

        /* xorshift PRNG
         * Taken from https://en.wikipedia.org/wiki/Xorshift#Example_implementation
         * Used under CC-By-SA */

            int firstzero=0, bothzero=0;
            int[] arr = new int[n+1];
            int i=0, j=0;
            x=(int)(System.currentTimeMillis()/1000l);

            for(i=0; i< 1<<(n+1) ; i++) {
                int tmp=i;
                for(j=0; j<n+1; j++) {
                    arr[j]=(tmp&1)*(-2)+1;
                    tmp>>=1;
                }
                for(j=0; j<iters; j++) {
                    int[] randArr = new int[n];
                    int k=0;
                    long flag = 0;
                    int first=0, second=0;
                    do {
                        for(k=0; k<n; k++) {
                            randArr[k]=(1-(myRand()&3))%2;
                            flag+=(randArr[k]&1);
                            first+=arr[k]*randArr[k];
                            second+=arr[k+1]*randArr[k];
                        }
                    } while(allzero(randArr));
                    if( first == 0 )
                    {
                        firstzero+=1;
                        if( second == 0 )
                        {
                            bothzero++;
                        }
                    }
                }
            }
         return ( "firstzero " + firstzero + "\nbothzero " + bothzero + "\n" );
    }

    private boolean allzero(int[] arr)
    {
       for(int x : arr)
       {
          if(x!=0)
          {
             return false;
          }
       }
       return true;
    }

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

И я попытался запустить код Python после обновления Python и установки Python-Numpy, но я получаю это:

server:/tmp# python tester.py
Traceback (most recent call last):
  File "peepee.py", line 15, in <module>
    F = np.random.choice(np.array([-1,0,0,1], dtype=np.int8), size = n)
AttributeError: 'module' object has no attribute 'choice'

Golang версия 45X python на моей машине ниже кодов Golang:

package main

import (
"fmt"
"time"
)

const (
n     = 6
iters = 1000
)

var (
x, y, z, w = 34353, 34353, 57768, 1564 //PRNG seeds
)

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

func main() {
var firstzero, bothzero int
var arr [n + 1]int
var i, j int
x = int(time.Now().Unix())

for i = 0; i < 1<<(n+1); i = i + 1 {
    tmp := i
    for j = 0; j < n+1; j = j + 1 {
        arr[j] = (tmp&1)*(-2) + 1
        tmp >>= 1
    }
    for j = 0; j < iters; j = j + 1 {
        var randArr [n]int
        var flag uint
        var k, first, second int
        for {
            for k = 0; k < n; k = k + 1 {
                randArr[k] = (1 - (myRand() & 3)) % 2
                flag += uint(randArr[k] & 1)
                first += arr[k] * randArr[k]
                second += arr[k+1] * randArr[k]
            }
            if flag != 0 {
                break
            }
        }
        if first == 0 {
            firstzero += 1
            if second == 0 {
                bothzero += 1
            }
        }
    }
}
println("firstzero", firstzero, "bothzero", bothzero)
}

и приведенные ниже коды Python, скопированные сверху:

import itertools
import operator
import random

n=6
iters = 1000
firstzero = 0
bothzero = 0

choicesF = filter(any, itertools.product([-1, 0, 0, 1], repeat=n))

for S in itertools.product([-1,1], repeat = n+1):
    for i in xrange(iters):
        F = random.choice(choicesF)
        if not sum(map(operator.mul, F, S[:-1])):
            firstzero += 1
            if not sum(map(operator.mul, F, S[1:])):
                bothzero += 1

print "firstzero", firstzero
print "bothzero", bothzero

и время ниже:

$time python test.py
firstzero 27349
bothzero 12125

real    0m0.477s
user    0m0.461s
sys 0m0.014s

$time ./hf
firstzero 27253 bothzero 12142

real    0m0.011s
user    0m0.008s
sys 0m0.002s

C # 0.135s

C # на основе простого питона Алистера Бакстона : 0,278 с
Параллелизированный C #: 0,135 с
Python из вопроса: 5.907
с простого питона Алистера: 0,853 с

Я на самом деле не уверен, что эта реализация верна - ее вывод отличается, если вы посмотрите на результаты внизу.

Там, безусловно, более оптимальные алгоритмы. Я просто решил использовать алгоритм, очень похожий на алгоритм Python.

Однопоточный C

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace ConvolvingArrays
{
    static class Program
    {
        static void Main(string[] args)
        {
            int n=6;
            int iters = 1000;
            int firstzero = 0;
            int bothzero = 0;

            int[] arraySeed = new int[] {-1, 1};
            int[] randomSource = new int[] {-1, 0, 0, 1};
            Random rand = new Random();

            foreach (var S in Enumerable.Repeat(arraySeed, n+1).CartesianProduct())
            {
                for (int i = 0; i < iters; i++)
                {
                    var F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    while (!F.Any(f => f != 0))
                    {
                        F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    }
                    if (Enumerable.Zip(F, S.Take(n), (f, s) => f * s).Sum() == 0)
                    {
                        firstzero++;
                        if (Enumerable.Zip(F, S.Skip(1), (f, s) => f * s).Sum() == 0)
                        {
                            bothzero++;
                        }
                    }
                }
            }

            Console.WriteLine("firstzero {0}", firstzero);
            Console.WriteLine("bothzero {0}", bothzero);
        }

        // itertools.product?
        // http://ericlippert.com/2010/06/28/computing-a-cartesian-product-with-linq/
        static IEnumerable<IEnumerable<T>> CartesianProduct<T>
            (this IEnumerable<IEnumerable<T>> sequences)
        {
            IEnumerable<IEnumerable<T>> emptyProduct =
              new[] { Enumerable.Empty<T>() };
            return sequences.Aggregate(
              emptyProduct,
              (accumulator, sequence) =>
                from accseq in accumulator
                from item in sequence
                select accseq.Concat(new[] { item }));
        }
    }
}

Параллельный C #:

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace ConvolvingArrays
{
    static class Program
    {
        static void Main(string[] args)
        {
            int n=6;
            int iters = 1000;
            int firstzero = 0;
            int bothzero = 0;

            int[] arraySeed = new int[] {-1, 1};
            int[] randomSource = new int[] {-1, 0, 0, 1};

            ConcurrentBag<int[]> results = new ConcurrentBag<int[]>();

            // The next line iterates over arrays of length n+1 which contain only -1s and 1s
            Parallel.ForEach(Enumerable.Repeat(arraySeed, n + 1).CartesianProduct(), (S) =>
            {
                int fz = 0;
                int bz = 0;
                ThreadSafeRandom rand = new ThreadSafeRandom();
                for (int i = 0; i < iters; i++)
                {
                    var F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    while (!F.Any(f => f != 0))
                    {
                        F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    }
                    if (Enumerable.Zip(F, S.Take(n), (f, s) => f * s).Sum() == 0)
                    {
                        fz++;
                        if (Enumerable.Zip(F, S.Skip(1), (f, s) => f * s).Sum() == 0)
                        {
                            bz++;
                        }
                    }
                }

                results.Add(new int[] { fz, bz });
            });

            foreach (int[] res in results)
            {
                firstzero += res[0];
                bothzero += res[1];
            }

            Console.WriteLine("firstzero {0}", firstzero);
            Console.WriteLine("bothzero {0}", bothzero);
        }

        // itertools.product?
        // http://ericlippert.com/2010/06/28/computing-a-cartesian-product-with-linq/
        static IEnumerable<IEnumerable<T>> CartesianProduct<T>
            (this IEnumerable<IEnumerable<T>> sequences)
        {
            IEnumerable<IEnumerable<T>> emptyProduct =
              new[] { Enumerable.Empty<T>() };
            return sequences.Aggregate(
              emptyProduct,
              (accumulator, sequence) =>
                from accseq in accumulator
                from item in sequence
                select accseq.Concat(new[] { item }));
        }
    }

    // http://stackoverflow.com/a/11109361/1030702
    public class ThreadSafeRandom
    {
        private static readonly Random _global = new Random();
        [ThreadStatic]
        private static Random _local;

        public ThreadSafeRandom()
        {
            if (_local == null)
            {
                int seed;
                lock (_global)
                {
                    seed = _global.Next();
                }
                _local = new Random(seed);
            }
        }
        public int Next()
        {
            return _local.Next();
        }
        public int Next(int maxValue)
        {
            return _local.Next(maxValue);
        }
    }
}

Тестовый вывод:

Windows (.NET)

C # намного быстрее в Windows. Вероятно, потому что .NET быстрее, чем моно.

Время пользователя и системы не работает (используется git bashдля синхронизации).

$ time /c/Python27/python.exe numpypython.py
firstzero 27413
bothzero 12073

real    0m5.907s
user    0m0.000s
sys     0m0.000s
$ time /c/Python27/python.exe plainpython.py
firstzero 26983
bothzero 12033

real    0m0.853s
user    0m0.000s
sys     0m0.000s
$ time ConvolvingArrays.exe
firstzero 28526
bothzero 6453

real    0m0.278s
user    0m0.000s
sys     0m0.000s
$ time ConvolvingArraysParallel.exe
firstzero 28857
bothzero 6485

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

Linux (моно)

bob@phoebe:~/convolvingarrays$ time python program.py
firstzero 27059
bothzero 12131

real    0m11.932s
user    0m11.912s
sys     0m0.012s
bob@phoebe:~/convolvingarrays$ mcs -optimize+ -debug- program.cs
bob@phoebe:~/convolvingarrays$ time mono program.exe
firstzero 28982
bothzero 6512

real    0m1.360s
user    0m1.532s
sys     0m0.872s
bob@phoebe:~/convolvingarrays$ mcs -optimize+ -debug- parallelprogram.cs
bob@phoebe:~/convolvingarrays$ time mono parallelprogram.exe
firstzero 28857
bothzero 6496

real    0m0.851s
user    0m2.708s
sys     0m3.028s

Haskell: ~ 2000x ускорение на ядро

Скомпилируйте с 'ghc -O3 -funbox-strict-fields -threaded -fllvm' и запустите с '+ RTS -Nk', где k - количество ядер на вашем компьютере.

import Control.Parallel.Strategies
import Data.Bits
import Data.List
import Data.Word
import System.Random

n = 6 :: Int
iters = 1000 :: Int

data G = G !Word !Word !Word !Word deriving (Eq, Show)

gen :: G -> (Word, G)
gen (G x y z w) = let t  = x `xor` (x `shiftL` 11)
                      w' = w `xor` (w `shiftR` 19) `xor` t `xor` (t `shiftR` 8)
                  in (w', G y z w w')  

mask :: Word -> Word
mask = (.&.) $ (2 ^ n) - 1

gen_nonzero :: G -> (Word, G)
gen_nonzero g = let (x, g') = gen g 
                    a = mask x
                in if a == 0 then gen_nonzero g' else (a, g')


data F = F {zeros  :: !Word, 
            posneg :: !Word} deriving (Eq, Show)

gen_f :: G -> (F, G)       
gen_f g = let (a, g')  = gen_nonzero g
              (b, g'') = gen g'
          in  (F a $ mask b, g'')

inner :: Word -> F -> Int
inner s (F zs pn) = let s' = complement $ s `xor` pn
                        ones = s' .&. zs
                        negs = (complement s') .&. zs
                    in popCount ones - popCount negs

specialised_convolve :: Word -> F -> (Int, Int)
specialised_convolve s f@(F zs pn) = (inner s f', inner s f) 
    where f' = F (zs `shiftL` 1) (pn `shiftL` 1)

ss :: [Word]
ss = [0..2 ^ (n + 1) - 1]

main_loop :: [G] -> (Int, Int)
main_loop gs = foldl1' (\(fz, bz) (fz', bz') -> (fz + fz', bz + bz')) . parMap rdeepseq helper $ zip ss gs
    where helper (s, g) = go 0 (0, 0) g
                where go k u@(fz, bz) g = if k == iters 
                                              then u 
                                              else let (f, g') = gen_f g
                                                       v = case specialised_convolve s f
                                                               of (0, 0) -> (fz + 1, bz + 1)
                                                                  (0, _) -> (fz + 1, bz)
                                                                  _      -> (fz, bz)
                                                   in go (k + 1) v g'

seed :: IO G                                        
seed = do std_g <- newStdGen
          let [x, y, z, w] = map fromIntegral $ take 4 (randoms std_g :: [Int])
          return $ G x y z w

main :: IO ()
main = (sequence $ map (const seed) ss) >>= print . main_loop

Рубин

Ruby (2.1.0) 0.277s
Ruby (2.1.1) 0.281s
Python (Алистер Бакстон) 0.330s
Python (alemi) 0.097s

n = 6
iters = 1000
first_zero = 0
both_zero = 0

choices = [-1, 0, 0, 1].repeated_permutation(n).select{|v| [0] != v.uniq}

def convolve(v1, v2)
  [0, 1].map do |i|
    r = 0
    6.times do |j|
      r += v1[i+j] * v2[j]
    end
    r
  end
end

[-1, 1].repeated_permutation(n+1) do |s|
  iters.times do
    f = choices.sample
    fs = convolve s, f
    if 0 == fs[0]
      first_zero += 1
      if 0 == fs[1]
        both_zero += 1
      end
    end
  end
end

puts 'firstzero %i' % first_zero
puts 'bothzero %i' % both_zero

поток не будет полным без PHP

В 6,6 раза быстрее

PHP v5.5.9 - 1.223 0.646 сек;

против

Python v2.7.6 - 8.072 сек

<?php

$n = 6;
$iters = 1000;
$firstzero = 0;
$bothzero = 0;

$x=time();
$y=34353;
$z=57768;
$w=1564; //PRNG seeds

function myRand() {
    global $x;
    global $y;
    global $z;
    global $w;
    $t = $x ^ ($x << 11);
    $x = $y; $y = $z; $z = $w;
    return $w = $w ^ ($w >> 19) ^ $t ^ ($t >> 8);
}

function array_cartesian() {
    $_ = func_get_args();
    if (count($_) == 0)
        return array();
    $a = array_shift($_);
    if (count($_) == 0)
        $c = array(array());
    else
        $c = call_user_func_array(__FUNCTION__, $_);
    $r = array();
    foreach($a as $v)
        foreach($c as $p)
            $r[] = array_merge(array($v), $p);
    return $r;
}

function rand_array($a, $n)
{
    $r = array();
    for($i = 0; $i < $n; $i++)
        $r[] = $a[myRand()%count($a)];
    return $r;
}

function convolve($a, $b)
{
    // slows down
    /*if(count($a) < count($b))
        return convolve($b,$a);*/
    $result = array();
    $w = count($a) - count($b) + 1;
    for($i = 0; $i < $w; $i++){
        $r = 0;
        for($k = 0; $k < count($b); $k++)
            $r += $b[$k] * $a[$i + $k];
        $result[] = $r;
    }
    return $result;
}

$cross = call_user_func_array('array_cartesian',array_fill(0,$n+1,array(-1,1)));

foreach($cross as $S)
    for($i = 0; $i < $iters; $i++){
        while(true)
        {
            $F = rand_array(array(-1,0,0,1), $n);
            if(in_array(-1, $F) || in_array(1, $F))
                break;
        }
        $FS = convolve($S, $F);
        if(0==$FS[0]) $firstzero += 1;
        if(0==$FS[0] && 0==$FS[1]) $bothzero += 1;
    }

echo "firstzero $firstzero\n";
echo "bothzero $bothzero\n";

• Использовал пользовательский генератор случайных чисел (украден из ответа C), PHP один отстой, а числа не совпадают
• convolve Функция немного упростилась, чтобы быть быстрее
• Проверка только для массива с нулями очень оптимизирована (см. $FИ $FSпроверки).

Выходы:

$ time python num.py 
firstzero 27050
bothzero 11990

real    0m8.072s
user    0m8.037s
sys 0m0.024s
$ time php num.php
firstzero 27407
bothzero 12216

real    0m1.223s
user    0m1.210s
sys 0m0.012s

Редактировать. Вторая версия скрипта работает только для 0.646 sec:

<?php

$n = 6;
$iters = 1000;
$firstzero = 0;
$bothzero = 0;

$x=time();
$y=34353;
$z=57768;
$w=1564; //PRNG seeds

function myRand() {
    global $x;
    global $y;
    global $z;
    global $w;
    $t = $x ^ ($x << 11);
    $x = $y; $y = $z; $z = $w;
    return $w = $w ^ ($w >> 19) ^ $t ^ ($t >> 8);
}

function array_cartesian() {
    $_ = func_get_args();
    if (count($_) == 0)
        return array();
    $a = array_shift($_);
    if (count($_) == 0)
        $c = array(array());
    else
        $c = call_user_func_array(__FUNCTION__, $_);
    $r = array();
    foreach($a as $v)
        foreach($c as $p)
            $r[] = array_merge(array($v), $p);
    return $r;
}

function convolve($a, $b)
{
    // slows down
    /*if(count($a) < count($b))
        return convolve($b,$a);*/
    $result = array();
    $w = count($a) - count($b) + 1;
    for($i = 0; $i < $w; $i++){
        $r = 0;
        for($k = 0; $k < count($b); $k++)
            $r += $b[$k] * $a[$i + $k];
        $result[] = $r;
    }
    return $result;
}

$cross = call_user_func_array('array_cartesian',array_fill(0,$n+1,array(-1,1)));

$choices = call_user_func_array('array_cartesian',array_fill(0,$n,array(-1,0,0,1)));

foreach($cross as $S)
    for($i = 0; $i < $iters; $i++){
        while(true)
        {
            $F = $choices[myRand()%count($choices)];
            if(in_array(-1, $F) || in_array(1, $F))
                break;
        }
        $FS = convolve($S, $F);
        if(0==$FS[0]){
            $firstzero += 1;
            if(0==$FS[1])
                $bothzero += 1;
        }
    }

echo "firstzero $firstzero\n";
echo "bothzero $bothzero\n";

Решение F #

Время выполнения составляет 0,030 с при компиляции в x86 на CLR Core i7 4 (8) при 3,4 ГГц

Я понятия не имею, если код правильный.

• Функциональная оптимизация (встроенное сгибание) -> 0,026 с
• Сборка через консольный проект -> 0,022 с
• Добавлен лучший алгоритм генерации массивов перестановок -> 0,018 с
• Mono для Windows -> 0,089 с
• Выполнение сценария Alistair's Python -> 0,259 с

let inline ffoldi n f state =
    let mutable state = state
    for i = 0 to n - 1 do
        state <- f state i
    state

let product values n =
    let p = Array.length values
    Array.init (pown p n) (fun i ->
        (Array.zeroCreate n, i)
        |> ffoldi n (fun (result, i') j ->
            result.[j] <- values.[i' % p]
            result, i' / p
        )
        |> fst
    )

let convolute signals filter =
    let m = Array.length signals
    let n = Array.length filter
    let len = max m n - min m n + 1

    Array.init len (fun offset ->
        ffoldi n (fun acc i ->
            acc + filter.[i] * signals.[m - 1 - offset - i]
        ) 0
    )

let n = 6
let iters = 1000

let next =
    let arrays =
        product [|-1; 0; 0; 1|] n
        |> Array.filter (Array.forall ((=) 0) >> not)
    let rnd = System.Random()
    fun () -> arrays.[rnd.Next arrays.Length]

let signals = product [|-1; 1|] (n + 1)

let firstzero, bothzero =
    ffoldi signals.Length (fun (firstzero, bothzero) i ->
        let s = signals.[i]
        ffoldi iters (fun (first, both) _ ->
            let f = next()
            match convolute s f with
            | [|0; 0|] -> first + 1, both + 1
            | [|0; _|] -> first + 1, both
            | _ -> first, both
        ) (firstzero, bothzero)
    ) (0, 0)

printfn "firstzero %i" firstzero
printfn "bothzero %i" bothzero

Q, 0,296 сегмента

n:6; iter:1000  /parametrization (constants)
c:n#0           /auxiliar constant (sequence 0 0.. 0 (n))
A:B:();         /A and B accumulates results of inner product (firstresult, secondresult)

/S=sequence with all arrays of length n+1 with values -1 and 1
S:+(2**m)#/:{,/x#/:-1 1}'m:|n(2*)\1 

f:{do[iter; F:c; while[F~c; F:n?-1 0 0 1]; A,:+/F*-1_x; B,:+/F*1_x];} /hard work
f'S               /map(S,f)
N:~A; +/'(N;N&~B) / ~A is not A (or A=0) ->bitmap.  +/ is sum (population over a bitmap)
                  / +/'(N;N&~B) = count firstResult=0, count firstResult=0 and secondResult=0

Q - ориентированный на коллекцию язык (kx.com)

Код переписан для использования идиоматического Q, но без других умных оптимизаций

Языки сценариев оптимизируют время программиста, а не время выполнения

• Q не лучший инструмент для этой проблемы

Первая попытка кодирования = не победитель, а разумное время (примерно в 30 раз быстрее)

• вполне конкурентоспособен среди переводчиков
• остановиться и выбрать другую проблему

ПРИМЕЧАНИЯ.-

• программа использует начальное число по умолчанию (повторяемые execs), чтобы выбрать другое начальное число для случайного использования генератора \S seed
• Результат дается в виде последовательности двух целых, поэтому во втором значении есть конечный i-суффикс 27421 12133i -> читается как (27241, 12133)
• Время не считая запуска переводчика. \t sentence время, затраченное на это предложение

Юлия: 12,149 6,929 с

Несмотря на их претензии на скорость , начальное время компиляции JIT сдерживает нас!

Обратите внимание, что следующий код Julia фактически является прямым переводом исходного кода Python (без оптимизации) как демонстрация того, что вы можете легко перенести свой опыт программирования на более быстрый язык;)

require("Iterators")

n = 6
iters = 1000
firstzero = 0
bothzero = 0

for S in Iterators.product(fill([-1,1], n+1)...)
    for i = 1:iters
        F = [[-1 0 0 1][rand(1:4)] for _ = 1:n]
        while all((x) -> round(x,8) == 0, F)
            F = [[-1 0 0 1][rand(1:4)] for _ = 1:n]
        end
        FS = conv(F, [S...])
        if round(FS[1],8) == 0
            firstzero += 1
        end
        if all((x) -> round(x,8) == 0, FS)
            bothzero += 1
        end
    end
end

println("firstzero ", firstzero)
println("bothzero ", bothzero)

редактировать

Бег с n = 8занимает 32,935 с. Учитывая, что сложность этого алгоритма O(2^n), то 4 * (12.149 - C) = (32.935 - C), где Cявляется константой, представляющей время компиляции JIT. Решив, Cмы находим это C = 5.2203, предполагая, что фактическое время выполнения n = 6составляет 6,929 с.

Rust, 6,6 мс, ускорение 1950x

Практически прямой перевод кода Алистера Бакстона на Rust. Я подумал об использовании нескольких ядер с rayon (бесстрашный параллелизм!), Но это не улучшило производительность, возможно, потому что она уже очень быстрая.

extern crate itertools;
extern crate rand;
extern crate time;

use itertools::Itertools;
use rand::{prelude::*, prng::XorShiftRng};
use std::iter;
use time::precise_time_ns;

fn main() {
    let start = precise_time_ns();

    let n = 6;
    let iters = 1000;
    let mut first_zero = 0;
    let mut both_zero = 0;
    let choices_f: Vec<Vec<i8>> = iter::repeat([-1, 0, 0, 1].iter().cloned())
        .take(n)
        .multi_cartesian_product()
        .filter(|i| i.iter().any(|&x| x != 0))
        .collect();
    // xorshift RNG is faster than default algorithm designed for security
    // rather than performance.
    let mut rng = XorShiftRng::from_entropy(); 
    for s in iter::repeat(&[-1, 1]).take(n + 1).multi_cartesian_product() {
        for _ in 0..iters {
            let f = rng.choose(&choices_f).unwrap();
            if f.iter()
                .zip(&s[..s.len() - 1])
                .map(|(a, &b)| a * b)
                .sum::<i8>() == 0
            {
                first_zero += 1;
                if f.iter().zip(&s[1..]).map(|(a, &b)| a * b).sum::<i8>() == 0 {
                    both_zero += 1;
                }
            }
        }
    }
    println!("first_zero = {}\nboth_zero = {}", first_zero, both_zero);

    println!("runtime {} ns", precise_time_ns() - start);
}

И Cargo.toml, так как я использую внешние зависимости:

[package]
name = "how_slow_is_python"
version = "0.1.0"

[dependencies]
itertools = "0.7.8"
rand = "0.5.3"
time = "0.1.40"

Сравнение скорости:

$ time python2 py.py
firstzero: 27478
bothzero: 12246
12.80user 0.02system 0:12.90elapsed 99%CPU (0avgtext+0avgdata 23328maxresident)k
0inputs+0outputs (0major+3544minor)pagefaults 0swaps
$ time target/release/how_slow_is_python
first_zero = 27359
both_zero = 12162
runtime 6625608 ns
0.00user 0.00system 0:00.00elapsed 100%CPU (0avgtext+0avgdata 2784maxresident)k
0inputs+0outputs (0major+189minor)pagefaults 0swaps

6625608 нс - около 6,6 мс. Это означает ускорение в 1950 раз. Здесь возможно много оптимизаций, но я стремился к удобству чтения, а не производительности. Одной из возможных оптимизаций будет использование массивов вместо векторов для хранения вариантов, поскольку они всегда будут иметь nэлементы. Также возможно использовать RNG, отличный от XorShift, поскольку Xorshift работает быстрее, чем CSPRNG HC-128 по умолчанию, но медленнее, чем наивный алгоритм PRNG.