pythonnumbersdiophantine

How can I efficiently find integer solutions (x ≠ y) to a Diophantine equation using Python?


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:

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}")

Solution

  • 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).