How can I generate three random integers that satisfy some condition?

import numpy as np


def generate_answer(n: int, low_limit:int, high_limit: int):
    while True:
        a = np.random.randint(low_limit, high_limit + 1, 1)[0]
        b = np.random.randint(low_limit, high_limit + 1, 1)[0]
        c = (n - 7 * a - 5 * b) / 3.0
        if int(c) == c and low_limit <= c <= high_limit:
            break

    return a, b, int(c)


if __name__ == "__main__":
    n = 30
    ans = generate_answer(low_limit=-5, high_limit=50, n=n)
    assert ans[0] * 7 + ans[1] * 5 + ans[2] * 3 == n
    print(ans)

If you select two of the numbers a, b, c, you know the third. In this case, I randomize ints for a, b, and I find c by c = (n - 7 * a - 5 * b) / 3.0.

Make sure c is an integer, and in the allowed limits, and we are done.

If it is not, randomize again.

---

If you want to generate all possibilities,

def generate_all_answers(n: int, low_limit:int, high_limit: int):
    results = []
    for a in range(low_limit, high_limit + 1):
        for b in range(low_limit, high_limit + 1):
            c = (n - 7 * a - 5 * b) / 3.0
            if int(c) == c and low_limit <= c <= high_limit:
                results.append((a, b, int(c)))

    return results

First, let me note that your task is underspecified in at least two respects:

1. The allowed range of the generated values is not specified. In particular, you don't specify whether the results may include negative integers.
2. The desired distribution of the generated values is not specified.

Normally, if not specified, one might assume that a uniform distribution on the set of possible solutions to the equation was expected (since it is, in a certain sense, the most random possible distribution on a given set). But a (discrete) uniform distribution is only possible if the solution set is finite, which it won't be if the range of results is unrestricted. (In particular, if (a, b, c) is a solution, then so is (a, b + 3k, c − 5k) for any integer k.) So if we interpret the task as asking for a uniform distribution with unlimited range, it's actually impossible!

---

On the other hand, if we're allowed to choose any distribution and range, the task becomes trivial: just make the generator always return a = −n, b = n, c = n. Clearly this is a solution to the equation (since −7n + 5n + 3n = (−7 + 5 + 3)n = 1n), and a degenerate distribution that assigns all probability mass to single point is still a valid probability distribution!

If you wanted a slightly less degenerate solution, you could pick a random integer k (using any distribution of your choice) and return a = −n, b = n + 3k, c = n − 5k. As noted above, this is also a solution to the equation for any k. Of course, this distribution is still somewhat degenerate, since the value of a is fixed.

If you want to let all return values be at least somewhat random, you could also pick a random h and return a = −n + h, b = n − 2h + 3k and c = n + h − 5k. Again, this is guaranteed to be a valid solution for any h and k, since it clearly satisfies the equation for h = k = 0, and it's also easy to see that increasing or decreasing either h or k will leave the value of the left-hand side of the equation unchanged.

In fact, it can be proved that this method can generate all possible solutions to the equation, and that each solution will correspond to a unique (h, k) pair! (One fairly intuitive way to see this is to plot the solutions in 3D space and observe that they form a regular lattice of points on a 2D plane, and that the vectors (+1, −2, +1) and (0, +3, −5) span this lattice.) If we pick h and k from some distribution that (at least in theory) assigns a non-zero probability to every integer, then we'll have a non-zero probability of returning any valid solution. So, at least for one somewhat reasonable interpretation of the task (unbounded range, any distribution with full support) the following code should solve the task efficiently:

from random import gauss

def random_solution(n):
    h = int(gauss(0, 1000))  # any distribution with full support on the integers will do
    k = int(gauss(0, 1000))
    return (-n + h, n - 2*h + 3*k, n + h - 5*k)

---

If the range of possible values is restricted, the problem becomes a bit trickier. On the positive side, if all values are bounded below (or above), then the set of possible solutions is finite, and so a uniform distribution exists on it. On the flip side, efficiently sampling this uniform distribution is not trivial.

One possible approach, which you've used yourself, is to first generate all possible solutions (assuming there's a finite number of them) and then sample from the list of solutions. We can do the solution generation fairly efficiently like this:

1. find all possible values of a for which the equation might have a solution,
2. for each such a, find all possible values of b for which there still have a solution,
3. for each such (a, b) pair, solve the equation for c and check if it's valid (i.e. an integer within the specified range), and
4. if yes, add (a, b, c) to the set of solutions.

The tricky part is step 2, where we want to calculate the range of possible b values. For this, we can make use of the observation that, for a given a, setting c to its smallest allowed value and solving the equation gives an upper bound for b (and vice versa).

In particular, solving the equation for a, b and c respectively, we get:

• a = (n − 5b − 3c) / 7
• b = (n − 7a − 3c) / 5
• c = (n − 7a − 5b) / 3

Given lower bounds on some of the values, we can use these solutions to compute corresponding upper bounds on the others. For example, the following code will generate all non-negative solutions efficiently (and can be easily modified to use a lower bound other than 0, if needed):

def all_nonnegative_solutions(n):
    a_min = b_min = c_min = 0
    a_max = (n - 5*b_min - 3*c_min) // 7
    for a in range(a_min, a_max + 1):
        b_max = (n - 7*a - 3*c_min) // 5
        for b in range(b_min, b_max + 1):
            if (n - 7*a - 5*b) % 3 == 0:
                c = (n - 7*a - 5*b) // 3
                yield (a, b, c)

We can then store the solutions in a list or a tuple and sample from that list:

from random import choice

solutions = tuple(all_nonnegative_solutions(30))
a, b, c = choice(solutions)

Ps. Apparently Python's random.choice is not smart enough to use reservoir sampling to sample from an arbitrary iterable, so we do need to store the full list of solutions even if we only want to sample from it once. Or, of course, we could always implement our own sampler:

def reservoir_choice(iterable):
    r = None
    n = 0
    for x in iterable:
        n += 1
        if randrange(n) == 0:
           r = x
    return r

a, b, c = reservoir_choice(all_nonnegative_solutions(30))

BTW, we could make the all_nonnegative_solutions function above a bit more efficient by observing that the (n - 7*a - 5*b) % 3 == 0 condition (which checks whether c = (n − 7a − 5b) / 3 is an integer, and thus a valid solution) is true for every third value of b. Thus, if we first calculated the smallest value of b that satisfies the condition for a given a (which can be done with a bit of modular arithmetic), we could iterate over b with a step size of 3 starting from that minimum value and skip the divisibility check entirely. I'll leave implementing that optimization as an exercise.