1

I want to solve the following problem with python, if possible with sympy.

Let n be a fixed positive number. Let p=(p_1,...p_n) be a fixed known vector of positive integers. Let d be a fixed, known positive integer. Let q=(q_1,...,q_n) be a vector of unknown nonnegative integers.

How can I get all the solutions of p.q=d?

Where . means dot product.

Actually I can solve this for each individual n. But I want to create a function

def F(n,p,d):
...
return result

Such that result is a, e.g., list of all solutions. Note that from the restrictions made above, there is a finite number of solutions for each triplet of data (n,p,d).

I can't figure a way to do this, so any suggestion will be appreciated.


Added.

Example: suppose n=3 (the case n=2 is trivial), p=(2,1,3), d=3. Then I would do something like

res=[]
for i in range (d):
    for j in range (d):
        k=d-p[0]*i-p[2]*j
        if k>=0:
            res.append([i,k,j])

Then res=[[0, 3, 0], [0, 0, 1], [1, 1, 0]] which is correct.

As you can imagine, the bigger n is, the more for loops I need if I want to follow the same idea. So I do not think this is a good way to do it for arbitrary n, say n=57 or whatever big enough...

PepeToro
  • 559
  • 4
  • 14
  • 1
    Check out `itertools.product`. – Joel Cornett Sep 11 '14 at 13:25
  • Minor cleanup: You don't need to pass `F` the value of `n`. `n = len(p)`. – Jonathan Eunice Sep 11 '14 at 13:26
  • I don't fully understand the math problem but generally if the amount of nested loops is dependent on an unknown variable, **[recursion](http://en.wikipedia.org/wiki/Recursion)** is the way to go. So you have a function with a loop and within the loop the function either calls itself with slightly altered parameters (`n-1` most likely here) or returns under some condition (often something like `n==0`) – asontu Sep 11 '14 at 13:26
  • @funkwurm: The amount of nested loops is dependent on the length of `p`, which can be determined (using `len(p)`) prior to entering any loops. – Joel Cornett Sep 11 '14 at 13:30
  • @JoelCornett sure, but you don't know `len(p)` when writing the python, do you? Either you or me are not understanding this problem, I'm wondering why neither OP's code nor your code is even remotely mentioning `q` when the question was all solution of `p·q=d` – asontu Sep 11 '14 at 13:33
  • @funkwurm: According to OP: "Let `p=(p_1, ..., p_n)` be a fixed *known* vector." – Joel Cornett Sep 11 '14 at 13:34
  • @funkwurm: Also, see edit. I have mentioned `q` – Joel Cornett Sep 11 '14 at 13:35
  • @funkwurm Well I have substituted 2 of the q's by i and j as they range over some integers and I check if the condition I want to satisfy is met. So, for the specific example (q1,q2,q3)=(i,k,j). – PepeToro Sep 11 '14 at 13:35
  • possible duplicate of [Find all combinations of coins when given some dollar value](http://stackoverflow.com/questions/1106929/find-all-combinations-of-coins-when-given-some-dollar-value) – David Eisenstat Sep 11 '14 at 13:41
  • @DavidEisenstat I'd say this is a bit more general as the "denominations" here are not fixed. If I understood correctly, there n=4, and that's it. But indeed it is very related, thanks for the suggestion. – PepeToro Sep 11 '14 at 14:12
  • @user58533 Most of the answers solve the general problem. – David Eisenstat Sep 11 '14 at 16:39
  • 1
    This is a constrained linear Diophantine equation. I don't know a general algorithm and wasn't able to find one in a web search. If you really need to solve this in generality (e.g., for large values of `d` or `n`), maybe your best bet is to treat it as an integer programming problem and call an IP solver from python and parse the output. – Robert Dodier Sep 11 '14 at 20:38
  • @RobertDodier. Unfortunately, my knowledge of algorithms is not as good as to fully comprehend your comment. If some context of my problem helps here it goes: I am computing normal forms of vector fields. For such a task one expands the vector field into its taylor series and constructs a function that kills some perturbative terms. In principle, I need to do this for each degree d and for arbitrary variables n. However, a closed solution is in general very difficult if not impossible and we are content with knowing stuff for small enough n and d, principally small d. – PepeToro Sep 12 '14 at 08:58
  • @user58533 since `d` is small, a direct search will be OK. Glad to hear you got the problem solved. – Robert Dodier Sep 12 '14 at 17:26

1 Answers1

3

Following the algorithm you provided:

from itertools import product
dot = lambda X, Y: sum(x * y for x, y in zip(X, Y))
p = [1, 2, 3, ...] # Whatever fixed value you have for `p`
d = 100 # Fixed d
results = []
for q in product(range(0, d+1), repeat=len(p)):
    if dot(p, q) == d:
        results.append(q)

However this is slightly inefficient since it is possible to determine prior to computing the entire dot product, whether k will be positive. So let's define the dot product like this:

def dot(X, Y, d):
    total = 0
    for x, y in zip(X, Y):
        total += x * y
        if total > d:
            return -1
    return total

Now, as soon as the total exceeds d, the calculation exits. You can also express this as a list comprehension:

results = [q for q in product(range(0, d+1), repeat=len(p)) if dot(p, q, d) == d]
Joel Cornett
  • 24,192
  • 9
  • 66
  • 88
  • Thanks a lot. In your first program, I think a minor adjustment is required `for q in product(range(0, d+1), repeat=len(p)):` and `if k==0`. The first is because I made a mistake in the post and q must be nonnegative and second because k, well k is the difference between the dot product and d. So whenever k=0 we have `p.q=d` – PepeToro Sep 11 '14 at 13:52
  • 2
    NP. See edits. No need for k, since you can just test `p.q == d` – Joel Cornett Sep 11 '14 at 13:58
  • 1
    This brute-force search requires `d^n` tests. That's OK for small values of `d` and `n`. If that's all OP needs, I guess all is well. Otherwise .... – Robert Dodier Sep 11 '14 at 20:33