https://twitter.com/loicaroyer/status/1146173537658404864
Given a list of integers u=[a_0, ..., a_m] find the set of indices {i_0, ..., i_n} for that list such that the product u[i_0]* ... *u[i_n] is the closest to a given integer N. The shortest and most #elegant solution wins. (standard libs allowed)
import numpy as np
import itertools as it# Setup inputs
# Use same inputs as randomly generated by @james_t_webber
# https://twitter.com/james_t_webber/status/1146203142079410178
m = 20
N = 797
u = np.array([5,9,19,22,25,27,28,31,33,34,36,43,44,68,72,81,86,92,96,98])
print('N: ',N)
print('u: ',u)N: 797
u: [ 5 9 19 22 25 27 28 31 33 34 36 43 44 68 72 81 86 92 96 98]
Calculate all the products and select the closest one. This takes on the order of seconds.
# Brute Force code by @james_t_webber adapted into a function
# This version computes the products of all combinations of u
# without repetition. Changed m = u.size
# https://twitter.com/james_t_webber/status/1146210405368266752
def closest_product_brute_force(N,u):
i = min(
map(np.array, it.product([False, True],repeat=u.size)),
key=lambda i:np.abs(N-u[i].prod())
)
i = np.nonzero(i)[0]
return i%timeit closest_product_brute_force(N,u)
i = closest_product_brute_force(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())1 loop, best of 3: 7.91 s per loop
i: [ 3 10]
u[i]: [22 36]
u[i].prod(): 792
Rather than computing all the product of all the combinations, we just need to compute the products of the integers closest to N. Once we find such a product we can stop.
To do this, we iterate over the integers closest to N starting from N working outwards.
- Call each integer M.
- Figure out which elements of u are divisors of M.
- Compute the possible combinations of those divisors.
- Compute the product of those divisor combinations.
- Stop (short-circuit) if the product equals M
- Choose M one step farther from N
To implement this we will use a special type of iterable called a generator. Generators use lazy evaluation which defers evaluation until the next element of the iterable is needed. This will help us define the iteration without having to actually compute anything until we need it.
This takes 100s of microseconds, a 10^4 speed-up over the brute force method for the given input.
def closest_product_no_repetition(N,u):
# Setup generator for product targets "M"
Mg = it.chain.from_iterable(
(
(N-y,N+y) if y != 0 else (N,)
for y in range(min(abs(N-u))+1)
)
)
# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
lambda f: f[0].size > 0,
(
(u[M % u == 0],M) for M in Mg
)
)
# For each tuple of (divisors,M),
# generate combinations of divisors
divisor_combinations = it.chain.from_iterable(zip(
it.chain.from_iterable(
(it.combinations(d,r+1) for r in range(d.size))
),
it.repeat(M)
) for d,M in divisors)
# We want to use it.takewhile, but we need info on the divisor combination that met the condition
# Use tee to duplicate the iterator and advance it by one
g = it.tee(divisor_combinations)
u_i = next(g[1])
g = zip(*g)
# Up to this point no products have been calculated, we have just set up iterators
# Now, start iterating through M, divisors, and divisior_combinations
# until the product equals M
out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
if len(out) > 0:
u_i = out[-1][1][0]
else:
u_i = list(u_i)[0]
u_i = np.array(u_i)
u_i = u_i[u_i != 1]
i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
return i%timeit i = closest_product_no_repetition(N,u)
i = closest_product_no_repetition(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())1000 loops, best of 3: 190 µs per loop
i: [3, 10]
u[i]: [22 36]
u[i].prod(): 792
It is not clear from the original problem description if elements of u can be repeated. This is easily addressed by just taking all the combinations with repetition.
def closest_product_with_repetition(N,u):
# Setup generator for product targets "M"
Mg = it.chain.from_iterable(
(
(N-y,N+y) if y != 0 else (N,)
for y in range(min(abs(N-u))+1)
)
)
# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
lambda f: f[0].size > 1,
(
# Append 1 to the list of divisors so that we can use a fixed combination size
(np.append(u[M % u == 0],1),M) for M in Mg
)
)
# For each tuple of (divisors,M),
# generate combinations of divisors with replacement
divisor_combinations = it.chain.from_iterable(zip(
it.combinations_with_replacement(
d,
# Since we allow for replacement,
# set combination length such that a power of the smallest divisor exceeds M
int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
),
it.repeat(M)
) for d,M in divisors)
# We want to use it.takewhile, but we need info on the divisor combination that met the condition
# Use tee to duplicate the iterator and advance it by one
g = it.tee(divisor_combinations)
u_i = next(g[1])
g = zip(*g)
# Up to this point no products have been calculated, we have just set up iterators
# Now, start iterating through M, divisors, and divisior_combinations
# until the product equals M
out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
if len(out) > 0:
u_i = out[-1][1][0]
else:
u_i = list(u_i)[0]
u_i = np.array(u_i)
u_i = u_i[u_i != 1]
i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
return i%timeit i = closest_product_with_repetition(N,u)
i = closest_product_with_repetition(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())1000 loops, best of 3: 1.05 ms per loop
i: [3, 10]
u[i]: [22 36]
u[i].prod(): 792
First we define a helper function, peek_with_label to help take a look at some of the generators used.
def peek_with_label(iterable,num=10,label=''):
iterable = it.tee(iterable)
print(label)
for x in it.islice(iterable[0],num):
print(x)
return iterable[1]# Setup generator for product targets "M"
Mg = it.chain.from_iterable(
(
(N-y,N+y) if y != 0 else (N,)
for y in range(min(abs(N-u))+1)
)
)
Mg = peek_with_label(Mg,10,'M generator:')M generator:
797
796
798
795
799
794
800
793
801
792
# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
lambda f: f[0].size > 1,
(
# Append 1 to the list of divisors so that we can use a fixed combination size
(np.append(u[M % u == 0],1),M) for M in Mg
)
)
divisors = peek_with_label(divisors,10,'Divisors: ')Divisors:
(array([19, 1]), 798)
(array([5, 1]), 795)
(array([ 5, 25, 1]), 800)
(array([9, 1]), 801)
(array([ 9, 22, 33, 36, 44, 72, 1]), 792)
(array([5, 1]), 790)
(array([5, 1]), 805)
(array([31, 1]), 806)
(array([5, 1]), 785)
(array([28, 98, 1]), 784)
# For each tuple of (divisors,M),
# generate combinations of divisors with replacement
divisor_combinations = it.chain.from_iterable(zip(
it.combinations_with_replacement(
d,
# Since we allow for replacement,
# set combination length such that a power of the smallest divisor exceeds M
int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
),
it.repeat(M)
) for d,M in divisors)
divisor_combinations = peek_with_label(divisor_combinations,240,'Divisor combinations with repetition: ')Divisor combinations with repetition:
((19, 19, 19), 798)
((19, 19, 1), 798)
((19, 1, 1), 798)
((1, 1, 1), 798)
((5, 5, 5, 5, 5), 795)
((5, 5, 5, 5, 1), 795)
((5, 5, 5, 1, 1), 795)
((5, 5, 1, 1, 1), 795)
((5, 1, 1, 1, 1), 795)
((1, 1, 1, 1, 1), 795)
((5, 5, 5, 5, 5), 800)
((5, 5, 5, 5, 25), 800)
((5, 5, 5, 5, 1), 800)
((5, 5, 5, 25, 25), 800)
((5, 5, 5, 25, 1), 800)
((5, 5, 5, 1, 1), 800)
((5, 5, 25, 25, 25), 800)
((5, 5, 25, 25, 1), 800)
((5, 5, 25, 1, 1), 800)
((5, 5, 1, 1, 1), 800)
((5, 25, 25, 25, 25), 800)
((5, 25, 25, 25, 1), 800)
((5, 25, 25, 1, 1), 800)
((5, 25, 1, 1, 1), 800)
((5, 1, 1, 1, 1), 800)
((25, 25, 25, 25, 25), 800)
((25, 25, 25, 25, 1), 800)
((25, 25, 25, 1, 1), 800)
((25, 25, 1, 1, 1), 800)
((25, 1, 1, 1, 1), 800)
((1, 1, 1, 1, 1), 800)
((9, 9, 9, 9), 801)
((9, 9, 9, 1), 801)
((9, 9, 1, 1), 801)
((9, 1, 1, 1), 801)
((1, 1, 1, 1), 801)
((9, 9, 9, 9), 792)
((9, 9, 9, 22), 792)
((9, 9, 9, 33), 792)
((9, 9, 9, 36), 792)
((9, 9, 9, 44), 792)
((9, 9, 9, 72), 792)
((9, 9, 9, 1), 792)
((9, 9, 22, 22), 792)
((9, 9, 22, 33), 792)
((9, 9, 22, 36), 792)
((9, 9, 22, 44), 792)
((9, 9, 22, 72), 792)
((9, 9, 22, 1), 792)
((9, 9, 33, 33), 792)
((9, 9, 33, 36), 792)
((9, 9, 33, 44), 792)
((9, 9, 33, 72), 792)
((9, 9, 33, 1), 792)
((9, 9, 36, 36), 792)
((9, 9, 36, 44), 792)
((9, 9, 36, 72), 792)
((9, 9, 36, 1), 792)
((9, 9, 44, 44), 792)
((9, 9, 44, 72), 792)
((9, 9, 44, 1), 792)
((9, 9, 72, 72), 792)
((9, 9, 72, 1), 792)
((9, 9, 1, 1), 792)
((9, 22, 22, 22), 792)
((9, 22, 22, 33), 792)
((9, 22, 22, 36), 792)
((9, 22, 22, 44), 792)
((9, 22, 22, 72), 792)
((9, 22, 22, 1), 792)
((9, 22, 33, 33), 792)
((9, 22, 33, 36), 792)
((9, 22, 33, 44), 792)
((9, 22, 33, 72), 792)
((9, 22, 33, 1), 792)
((9, 22, 36, 36), 792)
((9, 22, 36, 44), 792)
((9, 22, 36, 72), 792)
((9, 22, 36, 1), 792)
((9, 22, 44, 44), 792)
((9, 22, 44, 72), 792)
((9, 22, 44, 1), 792)
((9, 22, 72, 72), 792)
((9, 22, 72, 1), 792)
((9, 22, 1, 1), 792)
((9, 33, 33, 33), 792)
((9, 33, 33, 36), 792)
((9, 33, 33, 44), 792)
((9, 33, 33, 72), 792)
((9, 33, 33, 1), 792)
((9, 33, 36, 36), 792)
((9, 33, 36, 44), 792)
((9, 33, 36, 72), 792)
((9, 33, 36, 1), 792)
((9, 33, 44, 44), 792)
((9, 33, 44, 72), 792)
((9, 33, 44, 1), 792)
((9, 33, 72, 72), 792)
((9, 33, 72, 1), 792)
((9, 33, 1, 1), 792)
((9, 36, 36, 36), 792)
((9, 36, 36, 44), 792)
((9, 36, 36, 72), 792)
((9, 36, 36, 1), 792)
((9, 36, 44, 44), 792)
((9, 36, 44, 72), 792)
((9, 36, 44, 1), 792)
((9, 36, 72, 72), 792)
((9, 36, 72, 1), 792)
((9, 36, 1, 1), 792)
((9, 44, 44, 44), 792)
((9, 44, 44, 72), 792)
((9, 44, 44, 1), 792)
((9, 44, 72, 72), 792)
((9, 44, 72, 1), 792)
((9, 44, 1, 1), 792)
((9, 72, 72, 72), 792)
((9, 72, 72, 1), 792)
((9, 72, 1, 1), 792)
((9, 1, 1, 1), 792)
((22, 22, 22, 22), 792)
((22, 22, 22, 33), 792)
((22, 22, 22, 36), 792)
((22, 22, 22, 44), 792)
((22, 22, 22, 72), 792)
((22, 22, 22, 1), 792)
((22, 22, 33, 33), 792)
((22, 22, 33, 36), 792)
((22, 22, 33, 44), 792)
((22, 22, 33, 72), 792)
((22, 22, 33, 1), 792)
((22, 22, 36, 36), 792)
((22, 22, 36, 44), 792)
((22, 22, 36, 72), 792)
((22, 22, 36, 1), 792)
((22, 22, 44, 44), 792)
((22, 22, 44, 72), 792)
((22, 22, 44, 1), 792)
((22, 22, 72, 72), 792)
((22, 22, 72, 1), 792)
((22, 22, 1, 1), 792)
((22, 33, 33, 33), 792)
((22, 33, 33, 36), 792)
((22, 33, 33, 44), 792)
((22, 33, 33, 72), 792)
((22, 33, 33, 1), 792)
((22, 33, 36, 36), 792)
((22, 33, 36, 44), 792)
((22, 33, 36, 72), 792)
((22, 33, 36, 1), 792)
((22, 33, 44, 44), 792)
((22, 33, 44, 72), 792)
((22, 33, 44, 1), 792)
((22, 33, 72, 72), 792)
((22, 33, 72, 1), 792)
((22, 33, 1, 1), 792)
((22, 36, 36, 36), 792)
((22, 36, 36, 44), 792)
((22, 36, 36, 72), 792)
((22, 36, 36, 1), 792)
((22, 36, 44, 44), 792)
((22, 36, 44, 72), 792)
((22, 36, 44, 1), 792)
((22, 36, 72, 72), 792)
((22, 36, 72, 1), 792)
((22, 36, 1, 1), 792)
((22, 44, 44, 44), 792)
((22, 44, 44, 72), 792)
((22, 44, 44, 1), 792)
((22, 44, 72, 72), 792)
((22, 44, 72, 1), 792)
((22, 44, 1, 1), 792)
((22, 72, 72, 72), 792)
((22, 72, 72, 1), 792)
((22, 72, 1, 1), 792)
((22, 1, 1, 1), 792)
((33, 33, 33, 33), 792)
((33, 33, 33, 36), 792)
((33, 33, 33, 44), 792)
((33, 33, 33, 72), 792)
((33, 33, 33, 1), 792)
((33, 33, 36, 36), 792)
((33, 33, 36, 44), 792)
((33, 33, 36, 72), 792)
((33, 33, 36, 1), 792)
((33, 33, 44, 44), 792)
((33, 33, 44, 72), 792)
((33, 33, 44, 1), 792)
((33, 33, 72, 72), 792)
((33, 33, 72, 1), 792)
((33, 33, 1, 1), 792)
((33, 36, 36, 36), 792)
((33, 36, 36, 44), 792)
((33, 36, 36, 72), 792)
((33, 36, 36, 1), 792)
((33, 36, 44, 44), 792)
((33, 36, 44, 72), 792)
((33, 36, 44, 1), 792)
((33, 36, 72, 72), 792)
((33, 36, 72, 1), 792)
((33, 36, 1, 1), 792)
((33, 44, 44, 44), 792)
((33, 44, 44, 72), 792)
((33, 44, 44, 1), 792)
((33, 44, 72, 72), 792)
((33, 44, 72, 1), 792)
((33, 44, 1, 1), 792)
((33, 72, 72, 72), 792)
((33, 72, 72, 1), 792)
((33, 72, 1, 1), 792)
((33, 1, 1, 1), 792)
((36, 36, 36, 36), 792)
((36, 36, 36, 44), 792)
((36, 36, 36, 72), 792)
((36, 36, 36, 1), 792)
((36, 36, 44, 44), 792)
((36, 36, 44, 72), 792)
((36, 36, 44, 1), 792)
((36, 36, 72, 72), 792)
((36, 36, 72, 1), 792)
((36, 36, 1, 1), 792)
((36, 44, 44, 44), 792)
((36, 44, 44, 72), 792)
((36, 44, 44, 1), 792)
((36, 44, 72, 72), 792)
((36, 44, 72, 1), 792)
((36, 44, 1, 1), 792)
((36, 72, 72, 72), 792)
((36, 72, 72, 1), 792)
((36, 72, 1, 1), 792)
((36, 1, 1, 1), 792)
((44, 44, 44, 44), 792)
((44, 44, 44, 72), 792)
((44, 44, 44, 1), 792)
((44, 44, 72, 72), 792)
((44, 44, 72, 1), 792)
((44, 44, 1, 1), 792)
((44, 72, 72, 72), 792)
((44, 72, 72, 1), 792)
((44, 72, 1, 1), 792)
# We want to use it.takewhile, but we need info on the divisor combination that met the condition
# Use tee to duplicate the iterator and advance it by one
g = it.tee(divisor_combinations)
u_i = next(g[1])
g = zip(*g)
g = peek_with_label(g,10,'Paired divisor combinations:')Paired divisor combinations:
(((19, 19, 19), 798), ((19, 19, 1), 798))
(((19, 19, 1), 798), ((19, 1, 1), 798))
(((19, 1, 1), 798), ((1, 1, 1), 798))
(((1, 1, 1), 798), ((5, 5, 5, 5, 5), 795))
(((5, 5, 5, 5, 5), 795), ((5, 5, 5, 5, 1), 795))
(((5, 5, 5, 5, 1), 795), ((5, 5, 5, 1, 1), 795))
(((5, 5, 5, 1, 1), 795), ((5, 5, 1, 1, 1), 795))
(((5, 5, 1, 1, 1), 795), ((5, 1, 1, 1, 1), 795))
(((5, 1, 1, 1, 1), 795), ((1, 1, 1, 1, 1), 795))
(((1, 1, 1, 1, 1), 795), ((5, 5, 5, 5, 5), 800))
# Up to this point no products have been calculated, we have just set up iterators
# Now, start iterating through M, divisors, and divisior_combinations
# until the product equals M
out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
for x in out[-10:]:
print(x)(((22, 33, 1, 1), 792), ((22, 36, 36, 36), 792))
(((22, 36, 36, 36), 792), ((22, 36, 36, 44), 792))
(((22, 36, 36, 44), 792), ((22, 36, 36, 72), 792))
(((22, 36, 36, 72), 792), ((22, 36, 36, 1), 792))
(((22, 36, 36, 1), 792), ((22, 36, 44, 44), 792))
(((22, 36, 44, 44), 792), ((22, 36, 44, 72), 792))
(((22, 36, 44, 72), 792), ((22, 36, 44, 1), 792))
(((22, 36, 44, 1), 792), ((22, 36, 72, 72), 792))
(((22, 36, 72, 72), 792), ((22, 36, 72, 1), 792))
(((22, 36, 72, 1), 792), ((22, 36, 1, 1), 792))
if len(out) > 0:
u_i = out[-1][1][0]
else:
u_i = list(u_i)[0]
u_i = np.array(u_i)
u_i = u_i[u_i != 1]
i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())i: [3, 10]
u[i]: [22 36]
u[i].prod(): 792
def closest_product(N,u,repetition=False):
# Process input
u = np.array(u)
# Setup generator for product targets "M"
Mg = it.chain.from_iterable(
(
(N-y,N+y) if y != 0 else (N,)
for y in range(min(abs(N-u))+1)
)
)
if repetition:
# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
lambda f: f[0].size > 1,
(
# Append 1 to the list of divisors so that we can use a fixed combination size
(np.append(u[M % u == 0],1),M) for M in Mg
)
)
# For each tuple of (divisors,M),
# generate combinations of divisors with replacement
divisor_combinations = it.chain.from_iterable(zip(
it.combinations_with_replacement(
d,
# Since we allow for replacement,
# set combination length such that a power of the smallest divisor exceeds M
int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
),
it.repeat(M)
) for d,M in divisors)
else:
# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
lambda f: f[0].size > 0,
(
(u[M % u == 0],M) for M in Mg
)
)
# For each tuple of (divisors,M),
# generate combinations of divisors
divisor_combinations = it.chain.from_iterable(zip(
it.chain.from_iterable(
(it.combinations(d,r+1) for r in range(d.size))
),
it.repeat(M)
) for d,M in divisors)
# We want to use it.takewhile, but we need info on the divisor combination that met the condition
# Use tee to duplicate the iterator and advance it by one
g = it.tee(divisor_combinations)
u_i = next(g[1])
g = zip(*g)
# Up to this point no products have been calculated, we have just set up iterators
# Now, start iterating through M, divisors, and divisior_combinations
# until the product equals M
out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
if len(out) > 0:
u_i = out[-1][1][0]
else:
u_i = list(u_i)[0]
u_i = np.array(u_i)
u_i = u_i[u_i != 1]
i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
return iclosest_product(8,[2,3],repetition=False)[0, 1]
closest_product(8,np.array([2,3]),repetition=True)[0, 0, 0]
closest_product(8,[8],repetition=False)[0]
closest_product(8,[8],repetition=True)[0]
closest_product(N,u,repetition=False)[3, 10]
closest_product(N,u,repetition=True)[3, 10]
If we create a random challenging case, executing time is in the milliseconds
challenging_N = np.random.randint(101,2**20)
challenging_u = np.random.choice(range(2,1010),100,False)
print('challenging_N: ',challenging_N)
print('challenging_u: ',challenging_u)challenging_N: 121048
challenging_u: [ 955 96 249 501 167 940 607 923 186 871 118 470 772 773
112 805 349 605 890 423 692 335 251 663 82 959 344 882
16 602 488 454 383 520 320 311 804 408 52 93 65 137
62 723 544 160 894 510 356 338 414 364 868 579 202 898
827 492 273 846 215 306 615 479 565 405 981 512 474 380
168 682 391 535 38 100 551 872 702 191 572 353 250 467
30 26 866 188 36 967 743 398 709 886 193 1004 584 554
731 393]
%time i = closest_product(challenging_N,challenging_u,repetition=False)
i = closest_product(challenging_N,challenging_u,repetition=False)
print('i: ',i)
print('u[i]: ',challenging_u[i])
print('u[i].prod(): ',challenging_u[i].prod())CPU times: user 6.94 ms, sys: 0 ns, total: 6.94 ms
Wall time: 6.89 ms
i: [14, 84, 88]
u[i]: [112 30 36]
u[i].prod(): 120960