# Fastest algorithm to take the product of all subsets

Given n numbers in an array (you can't assume they are integers), I would like to compute the product of all subsets of size n-1.

You can do this by multiplying all the numbers together and then dividing by each one in turn, as long as none of the numbers is zero. However, how quickly can you do this with no division?

If you don't allow division, what is the minimum number of arithmetic operations (e.g. multiplication and addition) needed to compute the product of all subsets of size n-1?

Clearly you can do it in (n-1)*n multiplications.

To clarify, the output is n different products and the only operations apart from reading and writing to memory allowed are multiplication, addition and subtraction.

Example

If the input has three numbers 2,3,5, then the output is three numbers 15 = 3*5, 10 = 2*5 and 6 = 2*3.

Winning criterion

Answers should give an exact formula for the number of arithmetic operations their code will use in terms of n. To make life simple, I will just plug n = 1000 into your formula to judge its score. The lower the better.

If it's too hard to produce an exact formula for your code, you can just run it for n = 1000 and count the arithmetic operations in the code. An exact formula would be best however.

You should add your score for n=1000 to your answer for easy comparison.

• May we count multiplying by 1 as free? Otherwise I'd define a custom-multiplication function that does this.
– xnor
Jun 18, 2017 at 18:22
• Would it be against the rules to do a whole bunch of multiplications in parallel by concatenating numbers together with sufficiently many spacer 0 digits?
– xnor
Jun 18, 2017 at 18:37
• Do the operations such as + on indices count? If this is the case, does array indexing count as well? (since it is after all syntactic sugar for addition and dereferencing).
– nore
Jun 18, 2017 at 18:56
• @nore OK I give in :) Just count arithmetic operations that involve the input in some way.
– user7467
Jun 18, 2017 at 20:12
• Clearly you can do it in (n-1)*n multiplications You mean (n-2)*n, right? Jun 18, 2017 at 21:16

# Python, 3(n-2) operations, score = 2994

l = list(map(float, input().split()))
n = len(l)

left =  * len(l)
right =  * len(l)
left = l
right[-1] = l[-1]
for i in range(1,len(l)-1):
left[i] = l[i] * left[i - 1]
right[-i-1] = l[-i-1] * right[-i]

result =  * len(l)
result[-1] = left[-2]
result = right
for i in range(1, len(l) - 1):
result[i] = left[i - 1] * right[i+1]

print(result)


The arrays left and right contain the cumulated products of the array from the left and from the right, respectively.

## EDIT: Proof that 3(n-2) is the optimal number of operations needed for n >= 2, if we only use multiplication.

We will do that by induction; by the above algorithm, we just have to prove that for n >= 2, 3(n-2) is a lower bound on the number of multiplications needed.

For n = 2, we need at least 0 = 3(2-2) multiplications, so the result is trivial.

Let n > 2, and suppose for n - 1 elements, we need at least 3(n-3) multiplications. Consider a solution for n elements with k multiplications. Now, we remove the last of these elements as if it was 1, and simplify all multiplications directly by it. We also remove the multiplication leading to the product of all the other elements, since that one is not needed as it can never be used as an intermediate value to get the product of n-2 of the other elements, since division is not allowed. This leaves us with l multiplications, and a solution for n - 1 elements.

By induction hypothesis, we have l >= 3(n-3).

Now, let's have a look at how many multiplications were removed. One of them was the one leading to the product of all elements except the last. Moreover, the last element was used directly in at least two multiplications: if it was used in only one, then it was used when multiplying by an intermediate result consisting in some product of the other elements; let's say, without loss of generality, this this intermediate result included the first element in the product. Then, there is no way to get the product of all the elements but the first, since every product that contains the last element is either the last element, or contains the first element.

We thus have k >= l+3 >= 3(n-2), proving the claimed theorem.

• This turns out very clean in Haskell: f l = zipWith (*) (scanl (*) 1 l) (scanr (*) 1 $tail l). – xnor Jun 18, 2017 at 22:25 • Comments are not for extended discussion; this conversation has been moved to chat. Jun 21, 2017 at 3:19 # Haskell, score 2994 group :: Num a => [a] -> [[a]] group (a:b:t) = [a,b] : group t group [a] = [[a]] group [] = [] (%) :: (Num a, Eq a) => a -> a -> a a % 1 = a 1 % b = b a % b = a * b prod_one_or_two :: (Num a, Eq a) => [a] -> a prod_one_or_two [a, b] = a % b prod_one_or_two [x] = x insert_new_value :: (Num a, Eq a) => ([a], a) -> [a] insert_new_value ([a, b], c) = [c % b, c % a] insert_new_value ([x], c) = [c] products_but_one :: (Num a, Eq a) => [a] -> [a] products_but_one [a] =  products_but_one l = do combination <- combinations ; insert_new_value combination where pairs = group l subresults = products_but_one$ map prod_one_or_two pairs
combinations = zip pairs subresults


Try it online!

Say we're given the list [a,b,c,d,e,f,g,h]. We first group it into pairs [[a,b],[c,d],[e,f],[g,h]]. Then, we recurse on the half-size list pairs of their products to get subresults

[a*b, c*d, e*f, g*h] -> [(c*d)*(e*f)*(g*h), (a*b)*(e*f)*(g*h), (a*b)*(c*d)*(g*h), (a*b)*(c*d)*(e*f)]


If we take the first element (c*d)*(e*f)*(g*h), and multiply it by b and a respectively, we get the product of all but a and all but b. Doing this for every pair and recursive result with that pair missing, we get out final result. The odd-length case is specially handled by having the odd element passed unpaired to the recursive step, and the product of the remaining elements returned is the product without it.

The number of multiplications t(n) is n/2 for the pairing product, t(n/2) for the recursive call, and another n for the products with individual elements. This gives t(n) = 1.5 * n + t(n/2) for odd n. Using a more precise counts for odd n and ignoring multiplying with 1 for the base case gives score 2997 for n=1000.

• This is very nice.
– user7467
Jun 18, 2017 at 20:14
• I think the reason why the score is 2995 and not 2994 as in my answer is that it computes the product of all numbers as well in the non-power of two case, which is later truncated. Maybe a careful handling in products_but_one' could avoid that by returning something of the correct length.
– nore
Jun 18, 2017 at 20:19
• @nore I found I had an extra multiplication in my count because I forgot the base case has a 1 that is free to multiply. I think the padding 1 didn't affect things, but I cleaned up my algorithm to not use them.
– xnor
Jun 18, 2017 at 20:38
• Does this code assume the input is integers?
– user9206
Jun 19, 2017 at 3:52
• @Lembik It does, but only in the optional type annotations. I'll change them all to float.
– xnor
Jun 19, 2017 at 3:57

partition :: [Float] -> ([Float], [Float])
partition = foldr (\a (l1,l2) -> (l2, a:l1)) ([],[])

(%) :: Float -> Float -> Float
a % 1 = a
1 % b = b
a % b = a*b

merge :: (Float, [Float]) -> (Float, [Float]) -> (Float, [Float])
merge (p1,r1) (p2, r2) = (p1%p2, map(%p1)r2 ++ map(%p2)r1)

missing_products' :: [Float] -> (Float, [Float])
missing_products' [a] = (a,)
missing_products' l = merge res1 res2
where
(l1, l2) = partition l
res1 = missing_products' l1
res2 = missing_products' l2

missing_products :: [Float] -> [Float]
missing_products = snd . missing_products'


Try it online!

A divide-and-conquer strategy, very reminiscent of merge sort. Doesn't do any indexing.

The function partition splits the list into as-equal-as-possible halves by putting alternating elements on opposite sides of the partition. We recursively merge the results (p,r) for each of the halves, with r the list of products-with-one-missing, and p the overall product.

For the output for the full list, the missing element must be in one of the halves. The product with that element missing is a one-missing-product for the half it's in, multiplied by the full product for the other half. So, we multiply each product-with-one-missing by the full product of the other half and make a list of the results, as map(*p1)r2 ++ map(*p2)r1). This takes n multiplications, where n is the length. We also need to make a new full-product p1*p2 for future use, costing 1 more multiplication.

This gives the general recursion for for the number of operations t(n) with n even: t(n) = n + 1 + 2 * t(n/2). The odd one is similar, but one of the sublists is 1 larger. Doing out the recursion, we get n*(log_2(n) + 1) multiplications, though the odd/even distinction affects that exact value. The values up to t(3) are improved by not multiplying by 1 by defining a variant (%) of (*) that shortcuts the _*1 or 1*_ cases.

This gives 9975 multiplications for n=1000. I believe Haskell's laziness means the unused overall product in the outer layer is not computed for 9974; if I'm mistaken, I could omit it explicitly.

• You beat me by the timestamp one minute earlier.
– nore
Jun 18, 2017 at 18:15
• If it's hard to work out the formula exactly, feel free to just run it for n = 1000 and count the arithmetic operations in the code.
– user7467
Jun 18, 2017 at 18:18
• Since our code is basically the same, I don't understand how you got to 9974 and not 9975 multiplications for n = 1000 (in the case of computing the overall product in the outer layer). Did you include a 1 in the input you used to test it?
– nore
Jun 18, 2017 at 18:39
• @nore You're right, I was off by one. I wrote code to do the recursion for the number of multiplication function calls. Counting calls directly would be more reliable -- does anyone know how I'd do that in Haskell?
– xnor
Jun 18, 2017 at 18:43
• @xnor you can use trace from Debug.Trace with a catch-all | trace "call!" False = undefined guard, I think. But this uses unsafePerformIO under the hood, so it's not really that much of an improvement. Jun 19, 2017 at 15:11

## O(n log n) operations, score = 9974

Works with a binary tree.

## Python

l = list(map(int, input().split()))
n = len(l)

p =  * n + l
for i in range(n - 1, 1, -1):
p[i] = p[i + i] * p[i + i+1]

def mul(x, y):
if y == None:
return x
return x * y

r = [None] * n + [[None]] * n
for i in range(n - 1, 0, -1):
r[i] = [mul(p[i + i + 1], x) for x in r[i + i]] + [mul(p[i + i], x) for x in r[i + i + 1]]

u = r
j = 1
while j <= n:
j += j
print(u[n+n-j:] + u[:n+n-j])


This also requires list addition operations, and some arithmetic on numbers that are not the input values; not sure if that counts. The mul function is there to save n operations for the base case, to avoid wasting them by multiplying by 1. In any case, this is O(n log n) operations. The exact formula is, if only counting arithmetic operations on input numbers, with j = floor(log_2(n)): j * (2^(j + 1) - n) + (j + 1) * (2 * n - 2^(j + 1)) - 2.

Thanks to @xnor for saving one operation with the idea of not computing the outer product!

The last part is to output the products in order of the missing term.

• If it's hard to work out the formula exactly, feel free to just run it for n = 1000 and count the arithmetic operations in the code.
– user7467
Jun 18, 2017 at 18:17
• I counted 10975 operations...? Jun 18, 2017 at 18:29
• p[i] = p[i + i] * p[i + i+1] isn't counted Jun 18, 2017 at 18:29
• This is about n log2 n + n operations (which is O(nlogn) btw Jun 18, 2017 at 18:30
• @HyperNeutrino the operations in p[i] = p[i + i] * p[i + i + 1] should be saved by the multiplication optimization. I might have counted one too many, however.
– nore
Jun 18, 2017 at 18:34

group :: [a] -> Either [(a, a)] (a, [(a, a)])
group [] = Left []
group (a : l) = case group l of
Left pairs -> Right (a, pairs)
Right (b, pairs) -> Left ((a, b) : pairs)

products_but_one :: Num a => [a] -> [a]
products_but_one [_] = 
products_but_one [a, b] = [b, a]
products_but_one l = case group l of
Left pairs ->
let subresults =
products_but_one [a * b | (a, b) <- pairs]
in do ((a, b), c) <- zip pairs subresults; [c * b, c * a]
Right (extra, pairs) ->
let subresult : subresults =
products_but_one (extra : [a * b | (a, b) <- pairs])
in subresult : do ((a, b), c) <- zip pairs subresults; [c * b, c * a]


Try it online!

### How it works

This is a cleaned up version of xnor’s algorithm that deals with the odd case in a more straightforward way (edit: it looks like xnor has cleaned it up the same way):

[a, b, c, d, e, f, g] ↦
[a, bc, de, fg] ↦
[(bc)(de)(fg), a(de)(fg), a(bc)(fg), a(bc)(de)] by recursion ↦
[(bc)(de)(fg), a(de)(fg)c, a(de)(fg)b, a(bc)(fg)e, a(bc)(fg)d, a(bc)(de)g, a(bc)(de)f]

[a, b, c, d, e, f, g, h] ↦
[ab, cd, ef, gh] ↦
[(cd)(ef)(gh), (ab)(ef)(gh), (ab)(cd)(gh), (ab)(cd)(ef)] by recursion ↦
[(cd)(ef)(gh)b, (cd)(ef)(gh)a, (ab)(ef)(gh)d, (ab)(ef)(gh)c, (ab)(cd)(gh)f, (ab)(cd)(gh)e, (ab)(cd)(ef)h, (ab)(cd)(ef)g].

• "Given n numbers in an array (you can't assume they are integers), " We can't assume they are integers
– user9206
Jun 19, 2017 at 5:15

# O((n-2)*n) = O(n2): Trivial Solution

This is just the trivial solution that multiplies together each of the subsets:

# Python

def product(array): # Requires len(array) - 1 multiplication operations
if not array: return 1
result = array
for value in array[1:]:
result *= value
return result

def getSubsetProducts(array):
products = []
for index in range(len(array)): # calls product len(array) times, each time calling on an array of size len(array) - 1, which means len(array) - 2 multiplication operations called len(array) times
products.append(product(array[:index] + array[index + 1:]))
return products


Note that this also requires n list-addition operations; not sure if that counts. If that is not allowed, then product(array[:index] + array[index + 1:]) can be replaced to product(array[:index]) * product(array[index + 1:]), which changes the formula to O((n-1)*n).

– user7467
Jun 18, 2017 at 18:03
• Doesn't need your product function O(n) operations? one for each element in the array (althought this can easily be changed to O(n-1)) Jun 18, 2017 at 18:34
• @RomanGräf True. I'll change it to O(n-1) but thanks for pointing that out. Jun 18, 2017 at 18:42
• This has been changed to atomic-code-golf... Jun 18, 2017 at 19:33
• @EriktheOutgolfer What does that make my score now? Unless I am being blatantly stupid, don't the tag and the specs contradict each other now? Jun 18, 2017 at 19:44

# Python, 7540

A tripartite merge strategy. I think I can do even better than this, with a yet large merge. It's O(n log n).

EDIT: Fixed a miscount.

count = 0
def prod(a, b):
if a == 1: return b
if b == 1: return a
global count
count += 1
return a * b

def tri_merge(subs1, subs2, subs3):
total1, missing1 = subs1
total2, missing2 = subs2
total3, missing3 = subs3

prod12 = prod(total1, total2)
prod13 = prod(total1, total3)
prod23 = prod(total2, total3)

new_missing1 = [prod(m1, prod23) for m1 in missing1]
new_missing2 = [prod(m2, prod13) for m2 in missing2]
new_missing3 = [prod(m3, prod12) for m3 in missing3]

return prod(prod12, total3), new_missing1 + new_missing2 + new_missing3

def tri_partition(nums):
split_size = len(nums) // 3
a = nums[:split_size]
second_split_length = split_size + (len(nums) % 3 == 2)
b = nums[split_size:split_size + second_split_length]
c = nums[split_size + second_split_length:]
return a, b, c

def missing_products(nums):
if len(nums) == 1: return nums, 
if len(nums) == 0: return 1, []
subs = [missing_products(part) for part in tri_partition(nums)]
return tri_merge(*subs)

def verify(nums, res):
actual_product = 1
for num in nums:
actual_product *= num
actual_missing = [actual_product // num for num in nums]
return actual_missing == res and actual_product == res

nums = range(2, int(input()) + 2)
res = missing_products(nums)

print("Verified?", verify(nums, res))
if max(res) <= 10**10: print(res)

print(len(nums), count)


The relevant function is missing_products, which gives the overall product and all of the ones with a missing element.

• did you count the multiplications in tri_merge? Also you can replace the 2 * split_size + ... in tri_partition with split_size + split_size + .... Jun 18, 2017 at 19:51
• @RomanGräf I restructured it as per your suggestion. Jun 19, 2017 at 0:03

# dc, score 2994

#!/usr/bin/dc -f

# How it works:
# The required products are
#
#   (b × c × d × e × ... × x × y × z)
# (a) × (c × d × e × ... × x × y × z)
# (a × b) × (d × e × ... × x × y × z)
# ...
# (a × b × c × d × e × ... × x) × (z)
# (a × b × c × d × e × ... × x × y)
#
# We calculate each parenthesised term by
# multiplying the one above (on the left) or below
# (on the right), for 2(n-2) calculations, followed
# by the n-2 non-parenthesised multiplications
# giving a total of 3(n-2) operations.

?

# We will store input values into stack 'a' and
# accumulated product into stack 'b'.  Initialise
# stack b with the last value read.
sb

# Turnaround function at limit of recursion: print
# accumulated 'b' value (containing b..z above).
[Lbn[ ]nq]sG

# Recursive function - on the way in, we stack up
# 'a' values and multiply up the 'b' values.  On
# the way out, we multiply up the 'a' values and
# multiply each by the corresponding 'b' value.
[dSalb*Sb
z1=G
lFx
dLb*n[ ]n
La*]dsFx

# Do the last a*b multiplication
dLb*n[ ]n

# And we have one final 'a' value that doesn't have a
# corresponding 'b':
La*n


I'm assuming that the integer comparison z1= (which terminates the recursion when we reach the last value) is free. This is equivalent to the likes of foreach in other languages.

## Demonstrations

for i in '2 3 5' '2 3 5 7' '0 2 3 5' '0 0 1 2 3 4'
do printf '%s => ' "$i"; ./127147.dc <<<"$i"; echo
done
2 3 5 => 15 10 6
2 3 5 7 => 105 70 42 30
0 2 3 5 => 30 0 0 0
0 0 1 2 3 4 => 0 0 0 0 0 0


A demo with large and small inputs:

./127147.dc <<<'.0000000000000000000542101086242752217003726400434970855712890625 1 18446744073709551616'
18446744073709551616 1.0000000000000000000000000000000000000000000000000000000000000000 .0000000000000000000542101086242752217003726400434970855712890625


# C++, score: 5990, O([2NlogN]/3)

This implementation uses a binary tree look up table. My first implementation was O(NlogN), but a last minute optimization, which looks up the product of all array elements minus a pair, + 2 multiplications saved the day. I think this could still be optimized a bit further, maybe another 16%...

I've left some debugging traces in, only because it's easier to delete them than to rewrite them :)

 the actual complexity is measured at O([2NlogN]/3) for 100. It is actually a bit worse than O(NlogN) for small sets, but tends towards O([NlogN]/2) as the array grows very large O(0.57.NlogN) for a a set of 1 million elements.

#include "stdafx.h"
#include <vector>
#include <iostream>
#include <random>
#include <cstdlib>

using DataType = long double;

using DataVector = std::vector<DataType>;

struct ProductTree
{
std::vector<DataVector> tree_;
size_t ops_{ 0 };

ProductTree(const DataVector& v) : ProductTree(v.begin(), v.end()) {}
ProductTree(DataVector::const_iterator first, DataVector::const_iterator last)
{
Build(first, last);
}

void Build(DataVector::const_iterator first, DataVector::const_iterator last)
{
tree_.emplace_back(DataVector(first, last));

auto size = std::distance(first, last);
for (auto n = size; n >= 2; n >>= 1)
{
first = tree_.back().begin();
last = tree_.back().end();

DataVector v;
v.reserve(n);
while (first != last) // steps in pairs
{
auto x = *(first++);
if (first != last)
{
++ops_;
x *= *(first++); // could optimize this out,small gain
}
v.push_back(x);
}
tree_.emplace_back(v);
}
}

// O(NlogN) implementation...
DataVector Prod()
{
DataVector result(tree_.size());
for (size_t i = 0; i < tree_.size(); ++i)
{
auto depth = tree_.size() - 1;
auto k = i >> depth;
result[i] = ProductAtDepth(i, depth);
}
return result;
}

DataType ProductAtDepth(size_t index, size_t depth)
{
if (depth == 0)
{
return ((index ^ 1) < tree_[depth].size())
? tree_[depth][index ^ 1]
: 1;
}
auto k = (index >> depth) ^ 1;

if ((k < tree_[depth].size()))
{
++ops_;
return tree_[depth][k] * ProductAtDepth(index, depth - 1);
}
return ProductAtDepth(index, depth - 1);
}

// O([3NlogN]/2) implementation...
DataVector Prod2()
{
DataVector result(tree_.size());
for (size_t i = 0; i < tree_.size(); ++i)    // steps in pairs
{
auto depth = tree_.size() - 1;
auto k = i >> depth;
auto x = ProductAtDepth2(i, depth);
if (i + 1 < tree_.size())
{
ops_ += 2;
result[i + 1] = tree_[i] * x;
result[i] = tree_[i + 1] * x;
++i;
}
else
{
result[i] = x;
}
}
return result;
}

DataType ProductAtDepth2(size_t index, size_t depth)
{
if (depth == 1)
{
index = (index >> 1) ^ 1;
return (index < tree_[depth].size())
? tree_[depth][index]
: 1;
}
auto k = (index >> depth) ^ 1;

if ((k < tree_[depth].size()))
{
++ops_;
return tree_[depth][k] * ProductAtDepth2(index, depth - 1);
}
return ProductAtDepth2(index, depth - 1);
}

};

int main()
{
//srand(time());

DataVector data;
for (int i = 0; i < 1000; ++i)
{
auto x = rand() & 0x3;          // avoiding overflow and zero vaolues for testing
data.push_back((x) ? x : 1);
}

//for (int i = 0; i < 6; ++i)
//{
//  data.push_back(i + 1);
//}

//std::cout << "data:[";
//for (auto val : data)
//{
//  std::cout << val << ",";
//}
//std::cout << "]\n";

ProductTree pt(data);
DataVector result = pt.Prod2();

//std::cout << "result:[";
//for (auto val : result)
//{
//  std::cout << val << ",";
//}
//std::cout << "]\n";
std::cout << "N = " << data.size() << " Operations :" << pt.ops_ << '\n';

pt.ops_ = 0;
result = pt.Prod();

//std::cout << "result:[";
//for (auto val : result)
//{
//  std::cout << val << ",";
//}
//std::cout << "]\n";

std::cout << "N = " << data.size() << " Operations :" << pt.ops_ << '\n';

return 0;
}


I'm adding @nore's algorithm, for completeness. It's really nice, and is the fastest.

class ProductFlat
{
private:
size_t ops_{ 0 };

void InitTables(const DataVector& v, DataVector& left, DataVector& right)
{
if (v.size() < 2)
{
return;
}

left.resize(v.size() - 1);
right.resize(v.size() - 1);

auto l = left.begin();
auto r = right.rbegin();
auto ol = v.begin();
auto or = v.rbegin();

*l = *ol++;
*r = *or++;
if (ol == v.end())
{
return;
}

while (ol + 1 != v.end())
{
ops_ += 2;
*l = *l++ * *ol++;
*r = *r++ * *or++;
}
}

public:
DataVector Prod(const DataVector& v)
{
if (v.size() < 2)
{
return v;
}

DataVector result, left, right;
InitTables(v, left, right);

auto l = left.begin();
auto r = right.begin();
result.push_back(*r++);
while (r != right.end())
{
++ops_;
result.push_back(*l++ * *r++);
}
result.push_back(*l++);
return result;
}

auto Ops() const
{
return ops_;
}
};