import numpy
import scipy.signal

# to test, try the following
# > import random
# > xs = list(set(random.sample(xrange(10000), 500)))
# > subset_sums_up_to_u(xs,50000) == naive_subset_sums_up_to_u(xs,50000)

def naive_subset_sums_up_to_u(xs,u):
    sums = set([0])
    for x in xs:
        sums2 = set(sums)
        for y in sums:
            if x+y<=u :
                sums2.add(x+y)
        sums = sums2
    return sums

def subset_sums_up_to_u(xs,u):
    epsilon = 0.0001 # account for floating error
    # taking Minkowski sum of two sets, such that
    # the sets are contained in the union of intervals
    # [i*lc, i*uc] for all 0<=i<=k
    def smart_minkowski(A,B,k,lc,uc):
        n = min(1+u//lc,k)
        m = n*(uc-lc)
        gap = lc

        # in this case, we can be dumb and do normal minkowski
        # which is when the first coordinate is 0.
        if (lc<=m) or (m>=u):
            n = 0
            m = min(max(A)+max(B)+1,u)
            gap = u+1 # this make sure first coordinate is always 0

        AA = numpy.zeros((n+1,m+1))
        BB = numpy.zeros((n+1,m+1))

        for x in A:
            (i,j) = divmod(x, gap)
            AA[i,j] = 1
        for x in B:
            (i,j) = divmod(x, gap)
            BB[i,j] = 1
        CC = scipy.signal.fftconvolve(AA, BB)
        C = []

        # Find all the non-zeros
        (na,nb) = CC.shape
        for i in xrange(n+1):
            for j in xrange(m+1):
                if CC[i,j]>epsilon and i*lc+j<=u:
                    C.append(i*lc+j)
        return set(C)

    # taking minkowski sum of two sets with some extra information
    def minkowski((A,ka,la,ua),(B,kb,lb,ub)):
        lc = min(la,lb)
        uc = max(ua,ub)
        kc = ka+kb
        C = smart_minkowski(A,B,kc,lc,uc)
        return (C,kc,lc,uc)

    # combine a list, where each element of the list is
    # (set of subset sums, number of generators, lower bound of generators, upper bound of generators)
    def combine(xs):
        if len(xs)==1:
            return xs[0]
        evens = xs[::2]
        odds  = xs[1::2]
        extra = []
        if len(xs)%2 != 0:
            extra = [xs[-1]]

        ys = [minkowski(A,B) for (A,B) in zip(evens,odds)] + extra
        return combine(ys)

    # The elements in xs are called generators
    # Assume the input is actually a set
    ys = [(set([0,x]),1,x,x) for x in sorted(xs) if x<=u]
    # output is list of subset sums
    return combine(ys)[0]