Erdős-Straus-conjecture using polynomials in Python

Question

I am trying to write a code in Python to do the following. We can express the Erdős-Straus-conjecture in function of some polynomials $x(k), y(k), z(k) \in \mathbb{Q}[k]$ such that $\frac4k = \frac{1}{x(k)} + \frac{1}{y(k)} + \frac{1}{z(k)}$. Notice that the coefficients are rational numbers, but they have to be integers to be able to use them for this conjecture. Following formulas are given:

1. $ \frac4k = \frac1{uk} + \frac1{vk} + \frac{4uv - u - v}{uvk} \quad(P1)$
2. $ \frac4k = \frac1{uk} + \frac{4u-1}{uk+v} + \frac{(4u-1)v}{uk(uk+v)}\quad (P2) $
3. $ \frac4k = \frac{4u}{uk+v+1} + \frac4{k(uk+v+1)} + \frac{4v}{k(uk+v+1)}\quad (P3)$
4. $ \frac4k = \frac4{k+v} + \frac{4v}{k(k+4u+v)} + \frac{16uv}{k(k+v)(k+4u+v)} \quad (P4) $

The following theorem is used to determine when polynomials yield integer values: Consider $f \in \mathbb{Q}[k]$ en choose $m \in \mathbb{Z}$ such that $m \cdot f$ only has integer coefficients. Now consider $n = am+b \in \mathbb{Z}$ with $0 \le b < m$. Then $f(n) \in \mathbb{Z}$ iff $(m\cdot f)(b) \equiv 0\,\,(mod\,\,m)$.

This theorem implies that only looking at the residue classes modulo $m$ is sufficient to find out for which integers $n$ also $f(n)$ will be an integer.

For the Erdös-Straus-conjucture I have three polynomials $f_1(n), f_2(n), f_3(n)$. I want to find out for which $n$ all three polynomials have integer coefficients. Consider the corresponding $m$-values $m_1,m_2,m_3$ and let $m$ be equal to $lcm(m_1,m_2,m_3)$. Applying the above theorem (with the $m$ found), it's enough to check for which residue classes modulo $m$ all three polynomials $f_1,f_2,f_3$ will be integers.

So, I'd like to write a function in Python that has concrete values of $u$ and $v$ as input. I want to calculate the correct value for $m$ and to check for which residue classes modulo $m$ the polynomials have integer values. Therefore, the output should be a tuple ($m$,L) with L the list of residue classes modulo $m$ (that satisfy the conditions). I want to be able to use this function for (P1), (P2), (P3) and (P4).

This is what I've got so far:

P.<k> = PolynomialRing(QQ)
var('x,y,z')
def P1(u,v):
    x = u*k
    y = v*k
    z = (u*v*k)/(4*u*v-u-v)
    return (x,y,z)
def P2(u,v):
    x = u*k
    y = (u*k+v)/(4*u-1)
    z = (u*k*(u*k+v))/((4*u-1)*v)
    return (x,y,z)
def P3(u,v):
    x = (u*k+v+1)/(4*u)
    y = (k*(u*k+v+1))/4
    z = (k*(u*k+v+1))/(4*v)
    return (x,y,z)
def P4(u,v):
    x = (k+v)/4
    y = (k*(k+4*u+v))/(4*v)
    z = (k*(k+v)*(k+4*u+v))/(16*u*v)
    return (x,y,z)

var('m,k')
def findM(f):
    m = (f(k=1)^-1).denominator() 
    return m  
var('u,v,P,m1,m2,m3,m,k')
def mP(u,v,P):
    D = P(u,v)
    m1 = findM(D[0])
    m2 = findM(D[1])
    m3 = findM(D[2])
    m = lcm([m1,m2,m3])
    L = []
    for i in xrange(1,m-1):
         if((D[0](k=i)^-1).denominator() == 1 and (D[1](k=i)^-1).denominator() == 1 and (D[2](k=i)^-1).denominator() == 1 and Mod(m*((D[0]^-1) + (D[1]^-1) + (D[2]^-1))(k=b),m) == 0):
             L.append(i)
    return (m,L)

The problem is a keep getting an empty list. I'm using Sage Math to compute this (which is based on Python - kinda). Can someone please help me.

Infinitely grateful!

Answer 1

The following is a possible implementation, which you can test here.

R. = PolynomialRing(QQ)
D. = PolynomialRing(ZZ)
def P1(u,v):
    x = u*k
    y = v*k
    z = (u*v*k)/(4*u*v-u-v)
    return (x,y,z)
def P2(u,v):
    x = u*k
    y = (u*k+v)/(4*u-1)
    z = (u*k*(u*k+v))/((4*u-1)*v)
    return (x,y,z)
def P3(u,v):
    x = (u*k+v+1)/(4*u)
    y = (k*(u*k+v+1))/4
    z = (k*(u*k+v+1))/(4*v)
    return (x,y,z)
def P4(u,v):
    x = (k+v)/4
    y = (k*(k+4*u+v))/(4*v)
    z = (k*(k+v)*(k+4*u+v))/(16*u*v)
    return (x,y,z)
def findM(f):
    m = f.denominator() 
    return m
def mP(u,v,P):
    f = P(u,v)
    m1 = findM(f[0])
    m2 = findM(f[1])
    m3 = findM(f[2])
    m = lcm([m1,m2,m3])
    F = [D(m*f[0]),D(m*f[1]),D(m*f[2])]
    L = []
    for i in xrange(1,m+1):
         if(Mod(F[0](i),m)==0 and Mod(F[1](i),m)==0 and Mod(F[2](i),m)==0):
             L.append(i)
    return (m,L)
#Example:
[u,v,P] = [3/5,5,P4]
[m,S] = mP(u,v,P)
[x,y,z] = P(u,v)

print("[u,v,m] = ",[u,v,m])
print("[k,x(k),y(k),z(k),check]:")
for k in S:
    [xk,yk,zk] = [x(k),y(k),z(k)]
    chk = 1/xk + 1/yk + 1/zk
    print [k,xk,yk,zk,chk==4/k]

The example, using $(u,v)=(3/5,5)$ and $P4$, will return

[75, 20, 309, 10300, True]
[175, 45, 1596, 119700, True]
[375, 95, 7170, 1135250, True]
[475, 120, 11457, 2291400, True]
[675, 170, 23031, 6525450, True]
[775, 195, 30318, 9853350, True]
[975, 245, 47892, 19555900, True]
[1075, 270, 58179, 26180550, True]

The first line would mean $$ \frac{4}{75} = \frac{1}{20} + \frac{1}{309} + \frac{1}{10300} $$

---

You declared

P.<k> = PolynomialRing(QQ)

but later this removes the declaration:

var('u,v,P,m1,m2,m3,m,k')

You might want to assign a universal ring $R:=\mathbb Q[k]$ that is untouched throughout instead. Might also be a good idea to have an integral ring $D:=\mathbb Z[k]$ since you will work in it later, but I guess it's not necessary. i.e.

R.<k> = PolynomialRing(QQ)  
D.<k> = PolynomialRing(ZZ)

For finding $m$, it suffices to use

def findM(f):
    m = f.denominator() 
    return m  

The current command

(f(k=1)^-1).denominator()

will find the numerator instead, and evaluated at $f(1)$ so not the true denominator.

One last point is xrange(1,m-1) gives you the range $[1,2,\dots,m-2]$ but you want $[1,2,\dots,m]$.