Question: I'm trying to write a Python script to search for integer solutions (x, y) with x ≠ y to the following Diophantine equation: (y + n)^4 - y^4 = (x + k)^4 - x^4 Here:  n and k are fixed small positive integers (like n = 1, k = 2), x and y range from 1 to a large number (e.g., 1 to 1,000,000), I only want solutions where x ≠ y.  Aside from n = 74 and k = 24, for which we get positive integer solutions as desired, how can this code for instance be better written to generate any positive integer solutions to n,k,x,y? Or if possible what other integer solutions exist other than n = 74 and k = 24 ? I wrote a naive brute-force solution that checks every pair (x, y) in the range, but it's way too slow: n = 1 k = 2 limit = 10**6 for x in range(1, limit): rhs = (x + k)**4 - x**4 for y in range(1, limit): if x == y: continue lhs = (y + n)**4 - y**4 if lhs == rhs: print(f"Match found: x = {x}, y = {y}") 

Answer: Efficient equivalent rhs and lhs are increasing with increasing x and y, so you can go through both in parallel, always advancing the smaller one: n = 1 k = 2 limit = 10**6 x = y = 1 while x < limit and y < limit: rhs = (x + k)**4 - x**4 lhs = (y + n)**4 - y**4 if lhs == rhs: if x != y: print(f"Match found: x = {x}, y = {y}") x += 1 y += 1 elif lhs < rhs: y += 1 else: x += 1  Multiple n and k And to cover more values of n and k, you could merge the streams for all of them: from heapq import merge limit = 10**4 def values(k): for x in range(1, limit): yield (x + k)**4 - x**4, x, k lhs = None for rhs, x, k in merge(*map(values, range(1, 101))): if lhs == rhs: print(f'{n=}, {k=}, {x=}, {y=}') lhs, y, n = rhs, x, k  That takes about 1.5 seconds and barely any memory to find the same solutions as Kelly's search in the same ranges: n=74, k=24, x=134, y=59 n=75, k=25, x=133, y=59 n=27, k=5, x=497, y=271 n=63, k=35, x=257, y=193 n=64, k=36, x=256, y=193 n=54, k=10, x=994, y=542 n=81, k=15, x=1491, y=813  Unlimited You could also get rid of the limits (the one for x and y as well as the one for k and n) and leave it running as long as you have patience. Unlimited x is as easy as using for x in itertools.count(1): in my values function. Unlimited k: Instead of the merge function, use your own heap to merge the streams for the different k-values yourself. Start with just one, and whenever the stream for the current largest k goes from its x=1 to its x=2, add the stream for k+1. from heapq import heappush, heappop from time import time end_time = time() + 3 heap = [] def push(k, x): heappush(heap, ((x + k)**4 - x**4, k, x)) push(1, 1) lhs = None while time() < end_time: rhs, k, x = heappop(heap) if lhs == rhs: print(f'{n=}, {k=}, {x=}, {y=}') push(k, x + 1) if x == 1: push(k + 1, 1) lhs, y, n = rhs, x, k  I used a three seconds time limit there and it found 84 solutions: n=24, k=74, x=59, y=134 n=25, k=75, x=59, y=133 n=12, k=150, x=7, y=227 n=5, k=27, x=271, y=497 n=82, k=220, x=7, y=157 n=35, k=63, x=193, y=257 n=36, k=64, x=193, y=256 n=48, k=148, x=118, y=268 n=50, k=150, x=118, y=266 n=24, k=300, x=14, y=454 n=28, k=256, x=103, y=514 n=72, k=222, x=177, y=402 n=75, k=225, x=177, y=399 n=10, k=54, x=542, y=994 n=164, k=440, x=14, y=314 n=70, k=126, x=386, y=514 n=72, k=128, x=386, y=512 n=36, k=450, x=21, y=681 n=204, k=226, x=271, y=298 n=73, k=281, x=222, y=558 n=183, k=411, x=103, y=359 n=96, k=296, x=236, y=536 n=100, k=300, x=236, y=532 n=128, k=336, x=222, y=503 n=48, k=600, x=28, y=908 n=27, k=577, x=76, y=1176 n=120, k=370, x=295, y=670 n=125, k=375, x=295, y=665 n=15, k=81, x=813, y=1491 n=246, k=660, x=21, y=471 n=105, k=189, x=579, y=771 n=108, k=192, x=579, y=768 n=56, k=512, x=206, y=1028 n=60, k=750, x=35, y=1135 n=144, k=444, x=354, y=804 n=39, k=119, x=878, y=1342 n=150, k=450, x=354, y=798 n=20, k=108, x=1084, y=1988 n=328, k=880, x=28, y=628 n=168, k=518, x=413, y=938 n=140, k=252, x=772, y=1028 n=175, k=525, x=413, y=931 n=144, k=256, x=772, y=1024 n=72, k=900, x=42, y=1362 n=408, k=452, x=542, y=596 n=146, k=562, x=444, y=1116 n=366, k=822, x=206, y=718 n=192, k=592, x=472, y=1072 n=200, k=600, x=472, y=1064 n=84, k=768, x=309, y=1542 n=84, k=1050, x=49, y=1589 n=256, k=672, x=444, y=1006 n=25, k=135, x=1355, y=2485 n=410, k=1100, x=35, y=785 n=175, k=315, x=965, y=1285 n=180, k=320, x=965, y=1280 n=550, k=1100, x=76, y=653 n=216, k=666, x=531, y=1206 n=225, k=675, x=531, y=1197 n=96, k=1200, x=56, y=1816 n=5, k=679, x=604, y=5048 n=384, k=464, x=878, y=997 n=54, k=1154, x=152, y=2352 n=240, k=740, x=590, y=1340 n=250, k=750, x=590, y=1330 n=30, k=162, x=1626, y=2982 n=492, k=1320, x=42, y=942 n=210, k=378, x=1158, y=1542 n=216, k=384, x=1158, y=1536 n=108, k=1350, x=63, y=2043 n=112, k=1024, x=412, y=2056 n=264, k=814, x=649, y=1474 n=612, k=678, x=813, y=894 n=275, k=825, x=649, y=1463 n=219, k=843, x=666, y=1674 n=549, k=1233, x=309, y=1077 n=35, k=189, x=1897, y=3479 n=120, k=1500, x=70, y=2270 n=288, k=888, x=708, y=1608 n=78, k=238, x=1756, y=2684 n=574, k=1540, x=49, y=1099 n=300, k=900, x=708, y=1596 n=245, k=441, x=1351, y=1799 n=192, k=460, x=1324, y=1997  Note My second and third solutions compare the current value only with the one previous value. If a value appears even three or more times, then I'm only printing the matches of neighbors. For example if (x1,k1), (x2,k2) and (x3,k3) all match, I'm not showing the match of (x1,k1) with (x3,k3).