Alternatives to brute force algebra with additive/sigma functions

Question

I have an equation

       7405
210000 = Σ (x * (1 - 1/y)^n)
       n = 0

My goal is to find various x and y values that make it true. I also require for x and y to be whole numbers.

The way I have been doing it so far takes an incredible amount of time and is not very efficient... and right now I get no certain answers as they start to end in .999999999999

I am not sure if the numbers are not exactly precise which you get with floats on occasion or they actually have that many trailing 9s.

def closest_to_21000(a, b):
    diff_a = abs(a - 21000)
    diff_b = abs(b - 21000)

    if diff_a < diff_b:
        return a
    else:
        return b


def func():
    n = 7405 # for sigma
    maxResult = 0 # to keep track of our closest winners
    maxX = 0
    for x in range(1, 1001): # looping through X vals
        for y in range(1, 10001): # looping through Y vals
            result = sum(x * (1 - 1/y)**i for i in range(1, n+1)) # summization func
            if result > 21000: # if it's over 21,000 we can move to the next Y val
                break

            resu = closest_to_21000(maxResult, result) # is maxResult or new result is closer
            if maxResult == resu and x == maxX: # if that and x are the same as last one, new Y
                break
            maxResult = resu # max result = whichever was closer
            if maxResult == result: # if the new result was
                maxX = x # update maxX
                print(result, x, y)
            

func()
# closest answer I got was "20999.999999999996 700 31"

Is there a more efficient way I can do it? And a way to guarantee it's a whole number so I don't have to guess if it just has a trailing 0.99999999 after it? It takes a super long time if I want it to equal 21000000 instead of 21000. As of right now I am fine with it only checking x and y values between 1000 and 10000.

Answer 1

You can first rewrite your equation as:

             7405
210000 = x *  Σ (1 - 1/y)^n
             n = 0

You can then look up the formula for "finite sums of geometric sequences" to realize that the sum can be put into closed form. Let r = (1 - 1/y) and it is

(1 - r^(7406)) / (1 - r).  

So

r = (1 - 1/y)
result = x * (1 - r^7406) / (1 - r)

should be a lot faster to calculate.

---

Update:

I realized in the middle of the night that this can be made even simpler. Once you pull out the x from the summation, it's trivial to write x in terms of y.

You don't have to look through every possibility of x and y. Just look at a particular value of y, calculate x, and see if it is an integer.

It may also turn out that if you expand out the equation for x in terms of y, you'll find something interesting. I don't have pencil and paper nearby, but this may turn out to be a math problem rather than a programming problem.