How Slow Is Python Really (Part II)? This is a follow up to How slow is Python really? (Or how fast is your language?).

It turns out it was a bit too easy to get a x100 speedup for my last question. For those who have enjoyed the challenge but want something harder where they can really use their low level skills, here is part II.  The challenge is to get a x100 speedup for the following python code as tested on my computer.  

To make it more difficult I am using pypy this time. The current timing for me is 1 minute and 7 seconds using pypy 2.2.1.  

Rules 


The first person to submit code which I can run, is correct and is x100 times faster on my machine will be awarded a bounty of 50 points.
I will award the win to the fastest code after a week.


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") --&gt; Ax Ay Bx By Cx Cy Dx Dy
#   itertools.product("AB", repeat=5) --&gt; [
#    ('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 &lt; m and sum(map(operator.mul, F, S[j:j+n])) == 0):
            leadingzerocounts[j] +=1
            j += 1
print leadingzerocounts


The output should be similar to

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


You must use a random seed in your code and any random number generator that is good enough to give answers close to the above will be accepted.

My Machine The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor.  This also means I need to be able to run your code.

Explanation of code

This code iterates over all arrays S of length n+m-1 that are made up for -1s and 1s. For each array S it samples 1000 non-zero random arrays F of length n made up of -1,0 or 1 with probability of 1/4, 1/2,/14 of taking each values. It then computes the inner products between F and each window of S of length n until it finds a non-zero inner product. It adds 1 to leadingzerocounts at each position it found a zero inner product. 

Status 


Perl. 2.7 times slowdown by @tobyink. (Compared to pypy not cpython.)
J. 39 times speedup by @Eelvex.
C. 59 times speed up  by @ace.
Julia. 197 times faster not including start up time by @one-more-minute. 8.5 times speed up including start up time (it's faster using 4 processors in this case than 8).
Fortran. 438 times speed up by @semi-extrinsic. 
Rpython. 258 times speed up by @primo.
C++. 508 times speed up by @ilmale. 


(I stopped timing the new improvements because they are too fast and iters was too small.)



It was pointed out that timings below a second are unreliable and also some languages have a start-up cost. The argument is that if you are to include that you should also include the compilation time of C/C++ etc. Here are the timings for the fastest code with the number of iterations increased to 100,000.


Julia. 42 seconds by @one-more-minute. 
C++. 14 seconds by @GuySirton.
Fortran. 14s by @semi-extrinsic. 
C++. 12s by @ilmale.
Rpython. 18s by @primo.
C++. 5s by @Stefan.


The winner is.. Stefan!

Follow-up challenge posted. How high can you go? (A coding+algorithms challenge) . This one is  harder.
C++ x150 x450 x530

Instead of array I used bits (and dark magic).
Thanks @ace for the faster random function.

How does it work:
the first 15th bits of the integer s represent the array S[15]; the zeroes represent -1, the ones represent +1.
The array F is build in a similar way. But with two bit for each symbol.


00 represent -1
01 and 10 represent 0
11 represent 1


Cause S and F have a different representation I have to interleave S with itself to be comparable with F.


0 (-1) became 00 (-1 in the representation of F)
1 (+1) became 11 (+1 in the representation of F)


Now we can simply use Carnot to compute the inner product.
Remember that one variable can only assume value 00 or 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)

Looks like a not xor to me. :)

Sum the ones is just a game of shift and mask, nothing really complex.

#include &lt;array&gt;
#include &lt;ctime&gt;

// 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 &lt;&lt; 8)) &amp; B[3];
   x = (x | ( x &lt;&lt; 4)) &amp; B[2];
   x = (x | ( x &lt;&lt; 2)) &amp; B[1];
   x = (x | ( x &lt;&lt; 1)) &amp; B[0];
   return x | (x &lt;&lt; 1);
}

inline int32_t sumOnes(int32_t v)
{
   static int b[] = { 1, 0, 0, 1};
   int s = 0;
   for( int i = 0; i &lt; 8; ++i)
   {
      const int a = 3&amp;(v&gt;&gt;(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 &lt; 8; ++i)
   {
      const int a = 3&amp;(v&gt;&gt;(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&lt;&lt;1);
   x = y;
   y = z;
   z = w;
   w = w ^ ( w &gt;&gt; 19) ^ t ^ (t &gt;&gt; 8);
   return w;
}

int main()
{
   std::array&lt;int, 8&gt; leadingZero{0};
   x = static_cast&lt;int32_t&gt;(time(nullptr)); // seed PRNG
   const int maxS = 1 &lt;&lt; 15;
   for(int s = 0; s &lt; maxS; ++s)
   {
      const int32_t x = interleaveBit(s);
      for(int i = 0; i &lt; 1000; ++i)
      {
         int32_t random;
         do
         {
            random = 0xFFFF &amp; myRand();
         }while(sumOnes(random) == 0);
         int j = 7;
         while( j &gt;= 0 )
         {
            const int32_t h = (x &gt;&gt; (j*2));
            const int32_t l = 0xFFFF &amp; (~(random ^ h)); // inner product
            if(sumArray(l) == 0)
            {
               leadingZero[j]++;
            } else
            {
               break;
            }
            j--;
         }

      }
   }
   for(int i = 7; i &gt;= 0; --i)
   {
      printf("%d ", leadingZero[i]);
   }
   printf("\n");
   return 0;
}


Here a sample output:

6332350 2525218 1041716 438741 192917 87159 41023 19908 

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


The program has been compiled with:

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


on Fedora 20 with gcc 4.8.2
The Cpu is a i7 8core.

Probably I can gain some ms tweaking compiler parameters. 

While this is the OP solution time on my machine:

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

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




Edit:

Just adding openmp and change the order of the for I have a x3 gain, leading to a x450 performance improvement against OP code. :D
In this case the leadingZero array must be atomic. The random global... are random, they will be more random.

 #pragma omp parallel for
 for(int i = 0; i &lt; 1000; ++i)
 {
    int32_t random;
    do
    {
       random = 0xFFFF &amp; myRand();
    }while(sumOnes(random) == 0);
    for(int s = 0; s &lt; maxS; ++s)
    {
       const int32_t x = interleaveBit(s);
       int j = 7;
       while( j &gt;= 0 )
       {
          const int32_t h = (x &gt;&gt; (j*2));
          const int32_t l = 0xFFFF &amp; (~(random ^ h)); // inner product
          if( sumArray(l) == 0 )
          {
             leadingZero[j]++;
          } else
          {
             break;
          }
          j--;
       }
    }
 }


need to add -fopenmp to the compiler flag



Edit:2
As suggester by user71404 I changed the sumOnes and sumArray functions and now it's uber fast.

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


With openmp is slower, cause the atomics add too much overhead.

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


Without atomics is even faster, but I get wrong result.

2137992 1147218 619297 321243 155815 70946 32919 15579 

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


To understand sumArray consider that 16 bit represent and array of 8 numbers.
00 have no 1 and represent -1
01 and 10 have one 1 and represent 0
11 have two 1s and represent 1
So that built-in count the number of bit set to 1 [http://en.wikipedia.org/wiki/Hamming_weight] and to each group we remove 1. Cool.

sumOnes is just black magic.

Here the latest compile flags and code.

gcc -std=c++11 -mfpmath=sse -O3 -flto -march=native -funroll-loops -Wall -lstdc++

#include &lt;cstdint&gt;
#include &lt;cstdio&gt;
#include &lt;ctime&gt;

inline int32_t interleaveBit(int32_t x)
{
   static const uint32_t B[] = { 0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF };
   x = (x | ( x &lt;&lt; 8)) &amp; B[3];
   x = (x | ( x &lt;&lt; 4)) &amp; B[2];
   x = (x | ( x &lt;&lt; 2)) &amp; B[1];
   x = (x | ( x &lt;&lt; 1)) &amp; B[0];
   return x | (x &lt;&lt; 1);
}

inline int32_t sumOnes(int32_t v)
{
   /* 0xAAAA == 0b1010 1010 1010 1010 */
   return !!(0xAAAA &amp; (v ^ ~(v &lt;&lt; 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 &lt;&lt; 1);
   x = y;
   y = z;
   z = w;
   w = w ^ ( w &gt;&gt; 19) ^ t ^ (t &gt;&gt; 8);
   return w;
}

int main()
{
   int leadingZero[8] = { 0 };
   x = static_cast&lt;int32_t&gt;(time(nullptr)); // seed PRNG
   const int maxS = 1 &lt;&lt; 15;
   for( int i = 0; i &lt; 1000; ++i )
   {
      int32_t random;
      do
      {
         random = 0xFFFF &amp; myRand();
      } while(sumOnes(random) == 0 );
      for( int s = 0; s &lt; maxS; ++s )
      {
         const int32_t x = interleaveBit(s);
         int j = 7;
         while( j &gt;= 0 )
         {
            const int32_t h = (x &gt;&gt; (j * 2));
            const int32_t l = 0xFFFF &amp; (~(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;
}