pythonnumbersdiophantine

How can I find positive integers n,k,x such that a symmetric fourth-power equation is satisfied exactly?

I'm working on a mathematical problem that involves the following identity, and I want to find positive integers k,n,x such that all four terms in the equation are perfect fourth powers and the equation holds exactly.

The equation is as follows:

latex

( 
  -1/2 * ((1/3)^(1/3) * n)^2 / ( (3n^3 - 3(n^4 - x)/n + sqrt((1/3)n^8 + 9x^2)/n)^(1/3) )
  + 1/2 * n 
  + 1/2 * (1/3)^(1/3) * (3n^3 - 3(n^4 - x)/n + sqrt((1/3)n^8 + 9x^2)/n)^(1/3) 
)^4
+
( 
  -1/2 * (1/3)^(2/3) * k^2 / ( (3k^3 - 3(k^4 - x)/k + sqrt((1/3)k^8 + 9x^2)/k)^(1/3) )
  - 1/2 * k 
  + 1/2 * (1/3)^(1/3) * (3k^3 - 3(k^4 - x)/k + sqrt((1/3)k^8 + 9x^2)/k)^(1/3)
)^4
=
(
  -1/2 * ((1/3)^(1/3) * n)^2 / ( (3n^3 - 3(n^4 - x)/n + sqrt((1/3)n^8 + 9x^2)/n)^(1/3) )
  - 1/2 * n 
  + 1/2 * (1/3)^(1/3) * (3n^3 - 3(n^4 - x)/n + sqrt((1/3)n^8 + 9x^2)/n)^(1/3)
)^4
+
(
  -1/2 * (1/3)^(2/3) * k^2 / ( (3k^3 - 3(k^4 - x)/k + sqrt((1/3)k^8 + 9x^2)/k)^(1/3) )
  + 1/2 * k 
  + 1/2 * (1/3)^(1/3) * (3k^3 - 3(k^4 - x)/k + sqrt((1/3)k^8 + 9x^2)/k)^(1/3)
)^4

This identity appears complicated but is structurally symmetric and actually resembles a taxicab-type identity for fourth powers. It seems to reduce the usual 4-variable search space for equal fourth-power sums into just 3 variables: k,n,x .

I want to write Python code that:

When:

n = 24, k = 74, x = 300783360

The equation becomes:

158^4 + 59^4 = 134^4 + 133^4

Each term's 4th root is an exact integer, so this is a valid solution.

I also tried these codes which yielded integer values for k,n,x, which while they did satisfy the equation above, they did not do so such that the fourth root of each of the four terms that make up the above equation in total,were not positive non-zero integers. Also, while the error was negligible, it was still not zero as can be seen below.

You can view the full code here:

import math
import cmath

# Define cube roots as constants
cbrt_1_3 = (1 / 3) ** (1 / 3)
cbrt_1_3_sq = (1 / 3) ** (2 / 3)


def cbrt(x):
    """Calculate cube root, handling negative numbers"""
    if x >= 0:
        return x ** (1 / 3)
    else:
        return -(-x) ** (1 / 3)


def term_inner(var, x):
    """Calculate the inner expression for the cube root"""
    return 3 * var ** 3 - 3 * (var ** 4 - x) / var + math.sqrt((1 / 3) * var ** 8 + 9 * x ** 2) / var


def lhs_expr(n, k, x):
    """Left-hand side of the equation"""
    # Calculate the inner expressions
    inner_n = term_inner(n, x)
    inner_k = term_inner(k, x)

    # Calculate cube roots
    Rn = cbrt(inner_n)
    Rk = cbrt(inner_k)

    # First term
    A = (-0.5 * (cbrt_1_3 * n) ** 2 / Rn + 0.5 * n + 0.5 * cbrt_1_3 * Rn) ** 4

    # Second term - ERROR 1 FIXED: using cbrt_1_3_sq (not cbrt_1_3) as per original equation
    B = (-0.5 * cbrt_1_3_sq * k ** 2 / Rk - 0.5 * k + 0.5 * cbrt_1_3 * Rk) ** 4

    return A + B


def rhs_expr(n, k, x):
    """Right-hand side of the equation"""
    # Calculate the inner expressions
    inner_n = term_inner(n, x)
    inner_k = term_inner(k, x)

    # Calculate cube roots
    Rn = cbrt(inner_n)
    Rk = cbrt(inner_k)

    # First term - sign change on n term
    C = (-0.5 * (cbrt_1_3 * n) ** 2 / Rn - 0.5 * n + 0.5 * cbrt_1_3 * Rn) ** 4

    # Second term - ERROR 2 FIXED: using cbrt_1_3_sq (not cbrt_1_3) as per original equation
    D = (-0.5 * cbrt_1_3_sq * k ** 2 / Rk + 0.5 * k + 0.5 * cbrt_1_3 * Rk) ** 4

    return C + D


def find_integer_solution(n_max=20, k_max=20, x_max=100, tolerance=1e-8):
    """Find integer solutions to the equation where n ≠ k ≠ x"""
    print(f"Searching for solutions with n ≤ {n_max}, k ≤ {k_max}, x ≤ {x_max}")
    print("Constraint: n ≠ k ≠ x (all three values must be different)")

    solutions = []

    for n in range(1, n_max + 1):
        for k in range(1, k_max + 1):
            for x in range(1, x_max + 1):
                # Skip if any two values are equal
                if n == k or n == x or k == x:
                    continue

                try:
                    # Calculate LHS and RHS
                    lhs = lhs_expr(n, k, x)
                    rhs = rhs_expr(n, k, x)

                    # Check if values are real and close
                    if isinstance(lhs, complex) or isinstance(rhs, complex):
                        continue

                    error = abs(lhs - rhs)

                    if error < tolerance:
                        solution = {
                            'n': n,
                            'k': k,
                            'x': x,
                            'LHS': lhs,
                            'RHS': rhs,
                            'error': error
                        }
                        solutions.append(solution)
                        print(f"✅ Solution found: n={n}, k={k}, x={x}, error={error:.2e}")

                except (ZeroDivisionError, ValueError, TypeError) as e:
                    continue
                except Exception as e:
                    # Skip any other computational errors
                    continue

    return solutions


def verify_solution(n, k, x):
    """Verify a specific solution"""
    print(f"\nVerifying solution: n={n}, k={k}, x={x}")

    try:
        lhs = lhs_expr(n, k, x)
        rhs = rhs_expr(n, k, x)

        print(f"LHS value: {lhs}")
        print(f"RHS value: {rhs}")
        print(f"Difference: {abs(lhs - rhs)}")

        return abs(lhs - rhs) < 1e-8

    except Exception as e:
        print(f"Error in verification: {e}")
        return False


# Run the search
print("Starting search for integer solutions...")
solutions = find_integer_solution(n_max=15, k_max=15, x_max=50)

if solutions:
    print(f"\n🎉 Found {len(solutions)} solution(s)!")
    for i, sol in enumerate(solutions):
        print(f"\nSolution {i + 1}:")
        for key, val in sol.items():
            if key == 'error':
                print(f"  {key} = {val:.2e}")
            else:
                print(f"  {key} = {val}")
else:
    print("\n❌ No solutions found within the given ranges.")
    print("Try expanding the search ranges or checking the equation formulation.")

# Test a few specific cases with n ≠ k ≠ x
print("\n" + "=" * 50)
print("Testing specific cases where n ≠ k ≠ x:")

# Test cases where all three values are different
test_cases = [(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1), (1, 2, 4), (2, 3, 4), (3, 4, 5)]
for n, k, x in test_cases:
    print(f"\nTesting n={n}, k={k}, x={x}:")
    try:
        lhs = lhs_expr(n, k, x)
        rhs = rhs_expr(n, k, x)
        error = abs(lhs - rhs)
        print(f"  LHS = {lhs}")
        print(f"  RHS = {rhs}")
        print(f"  Error = {error:.2e}")
        if error < 1e-6:
            print("  ✅ This might be a solution!")
    except Exception as e:
        print(f"  Error: {e}")

# Alternative approach: Check for patterns with constraint n ≠ k ≠ x
print("\n" + "=" * 50)
print("Checking for patterns with n ≠ k ≠ x:")

count = 0
for n in range(1, 100):
    for k in range(1, 100):
        for x in range(1, 100):
            if n == k or n == x or k == x:
                continue
            try:
                lhs = lhs_expr(n, k, x)
                rhs = rhs_expr(n, k, x)
                error = abs(lhs - rhs)
                if error < 1e-6:
                    count += 1
                    print(f"Pattern found: n={n}, k={k}, x={x}, error={error:.2e}")
            except:
                continue

if count == 0:
    print("No patterns found in the tested range with n ≠ k ≠ x constraint.")

What I’m asking:

  1. How can I ensure that my Python code finds integer solutions for k,n,x such that all 4 terms are perfect fourth powers?

  2. Is there a better way (e.g., using mpmath or sympy) to symbolically or precisely evaluate this equation to avoid floating-point drift?

  3. Are there better numerical or algebraic strategies to speed up or validate such searches?

Any help or feedback is greatly appreciated!

Solution

  • As I mentioned in the comment, you haven't put in the check for each of the 4 numbers to be integers too. Thus there will be many non-integral solutions to the equation as well. Adding a check for that and another tolerance value does the job.

    I modified this section -

    def lhs_expr(n, k, x):
    """Left-hand side of the equation"""
    # Calculate the inner expressions
    inner_n = term_inner(n, x)
    inner_k = term_inner(k, x)
    
    # Calculate cube roots
    Rn = cbrt(inner_n)
    Rk = cbrt(inner_k)

    # First term
    smalla = (-0.5 * cbrt_1_3_sq * n ** 2 / Rn + 0.5 * n + 0.5 * cbrt_1_3 * Rn)
    A = smalla ** 4
    # Second term - ERROR 1 FIXED: using cbrt_1_3_sq (not cbrt_1_3) as per original equation
    smallb = (-0.5 * cbrt_1_3_sq * k ** 2 / Rk - 0.5 * k + 0.5 * cbrt_1_3 * Rk)
    B = smallb ** 4
    
    return A + B,smalla,smallb


def rhs_expr(n, k, x):
    """Right-hand side of the equation"""
    # Calculate the inner expressions
    inner_n = term_inner(n, x)
    inner_k = term_inner(k, x)

    # Calculate cube roots
    Rn = cbrt(inner_n)
    Rk = cbrt(inner_k)

    # First term - sign change on n term
    smallc = (-0.5 * cbrt_1_3_sq * n ** 2 / Rn - 0.5 * n + 0.5 * cbrt_1_3 * Rn)
    C = smallc ** 4

    # Second term - ERROR 2 FIXED: using cbrt_1_3_sq (not cbrt_1_3) as per original equation
    smalld = (-0.5 * cbrt_1_3_sq * k ** 2 / Rk + 0.5 * k + 0.5 * cbrt_1_3 * Rk)
    D = smalld ** 4
    
    return C + D,smallc,smalld

    and this section:

    count = 0
COUNTMAX = 1000
for n in range(24,25,1):
    for k in range(74,75,1):
        for x in range(1, 300790000):
            
            try:
                
                lhs,sa,sb = lhs_expr(n, k, x)
                
                rhs,sc,sd = rhs_expr(n, k, x)
                
                error = abs(lhs - rhs)

                if error < 1e-6:
                    # print(sa,check_magn(sa),sb,check_magn(sb),sc,check_magn(sc),sd,check_magn(sd))
                    if check_magn(sa) and check_magn(sb) and check_magn(sc) and check_magn(sd):
                        count += 1
                        print(f"Pattern found: n={n}, k={k}, x={x}, error={error:.2e} sa:{sa} sb:{sb} sc:{sc} sd:{sd}")
                
            except:
                continue
            if count>COUNTMAX:
                break
        if count>COUNTMAX:
            break
    if count>COUNTMAX:
        break

    in the code.

    which produced the output:

    ==================================================
Pattern found: n=24, k=74, x=5374176, error=1.86e-09 sa:48.999999999999986 sb:-25.000000000000014 sc:24.99999999999999 sd:48.999999999999986
Pattern found: n=24, k=74, x=300783360, error=1.19e-07 sa:157.99999999999991 sb:58.99999999999996 sc:133.99999999999991 sd:132.99999999999994

    (I searched with n=24 and k=74. Modify the search ranges as per your requirements.)

    The first solution in my search spacewas (n,k,x) = (24, 74, 5374176) which produces the 4 numbers (49,25,-25,49). If you want these four to be positive too you can add extra constraints.