Fastest algorithm to take the product of all subsets

Question

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.

Answer 1

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

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

left = [0] * len(l)
right = [0] * len(l)
left[0] = l[0]
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 = [0] * len(l)
result[-1] = left[-2]
result[0] = right[1]
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.

Answer 2

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] = [1]
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.

Answer 3

Haskell, score 9974

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,[1])
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.

Answer 4

O(n log n) operations, score = 9974

Works with a binary tree.

Python

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

p = [0] * 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[1]
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.

Answer 5

Haskell, score 2994

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 [_] = [1]
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].

Answer 6

O((n-2)*n) = O(n²): 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[0]
    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).

Answer 7

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[0], [1]
    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[1] and actual_product == res[0]

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

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

print(len(nums), count)

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

Answer 8

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.

# Read input from stdin
?

# 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

Answer 9

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 :)

[Edit] 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_[0].size());
        for (size_t i = 0; i < tree_[0].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_[0].size());
        for (size_t i = 0; i < tree_[0].size(); ++i)    // steps in pairs
        {
            auto depth = tree_.size() - 1;
            auto k = i >> depth;
            auto x = ProductAtDepth2(i, depth);
            if (i + 1 < tree_[0].size())
            {
                ops_ += 2;
                result[i + 1] = tree_[0][i] * x;
                result[i] = tree_[0][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_;
    }
};