Numpy Project Euler 1 Generalized Optimization

Question

I've started to learn Numpy and I'm looking for some ways to write this. I've written a generalization of Euler 1. It takes a list of divisors and a number, for example [3, 5] and 1000 in Euler 1.

Naively in pure python:

def subn1(divisors, n):
    return sum(i for i in range(1,n) if not all(i % d for d in divisors))

This runs in about 2.5 seconds for range(2,20), 1000000.

My first and best so far numpy attempt looks like this:

def subn2(divisors, n):
    a = np.arange(1,n) 
    b = np.zeros(a.shape[0], dtype=bool) 
    for d in divisors:
        b += a % d == 0
    return np.sum(a[b])

And runs in around 0.45 seconds for range(2,20), 1000000.

My third attempt was to remove the foor loop and use pure numpy, however it lost in the speed department by a small margin and uses more memory.

def subn3(divisors, n):
    nums = np.arange(1,n)     
    divs = np.array([divisors]).T
    return np.sum(nums[np.logical_or.reduce(np.mod(nums, divs) == 0, axis=0)])

This runs in around .5 seconds for range(2,20), 100000.

Is there a faster way to write it in 'pure' numpy or are foor loops not to shy away from?

Note: I'm aware that you can optimize it by reducing the divisor list, so there's no need to comment on that :)

Divakar · Answer 1 · 2015-12-06T11:26:19.107

One NumPythonic vectorized approach with broadcasting -

def subn_broadcasting(divisors,n):
    a = np.arange(1,n)
    return (a[(a % np.array(divisors)[:,None] == 0).any(0)]).sum()

Runtime tests and verify -

In [14]: # Inputs
    ...: n = 1000
    ...: divisors = range(2,20)
    ...: 

In [15]: print subn1(divisors, n)
    ...: print subn2(divisors, n)
    ...: print subn3(divisors, n)
    ...: print subn_broadcasting(divisors, n)
    ...: 
416056
416056
416056
416056

In [16]: %timeit subn1(divisors, n)
    ...: %timeit subn2(divisors, n)
    ...: %timeit subn3(divisors, n)
    ...: %timeit subn_broadcasting(divisors, n)
    ...: 
1000 loops, best of 3: 1.39 ms per loop
1000 loops, best of 3: 480 µs per loop
1000 loops, best of 3: 434 µs per loop
1000 loops, best of 3: 428 µs per loop

Hmm, doesn't look like much of an improvement there over the n2 and n3 versions.

score 0 · Answer 2 · answered Dec 06 '15 at 12:42

You can use np.where, like this:

def subn4(divisors, n):
    a = np.arange(np.min(divisors),n) 
    b = np.zeros(a.shape[0], dtype=bool) 
    for d in divisors:
        b += a % d == 0
    return np.sum(a[np.where(b)])

def subn4_(divisors, n):
    a = np.arange(1,n) 
    b = np.zeros(a.shape[0], dtype=bool) 
    for d in divisors:
        b += a % d == 0
    return np.sum(a[np.where(b)])

The tests, like suggested before:

%timeit subn1(divisors, n)
%timeit subn2(divisors, n)
%timeit subn3(divisors, n)
%timeit subn_broadcasting(divisors, n)
%timeit subn4(divisors, n)
%timeit subn4_(divisors, n)


1 loops, best of 3: 596 ms per loop
10 loops, best of 3: 30.1 ms per loop
10 loops, best of 3: 32 ms per loop
10 loops, best of 3: 31.9 ms per loop
10 loops, best of 3: 28.2 ms per loop
10 loops, best of 3: 27.4 ms per loop

Numpy Project Euler 1 Generalized Optimization

2 Answers2