How to find the maximum number of unique unit fractions that sum up to one

In my problem I have to write a program that calculates all the unique unit fractions (1/2 +1/3 + 1/6 = 1) that sum to one. I am currently at 570 fractions that sum to one in this format where the number is the denominator (this is my current solution):

(31,46,62,65,66,91,100,112,115,120,121,126,128,132,138,140,144,147,152,154,162,168,170,180,184,187,188,189,192,195,207,208,210,215,217,225,230,231,234,235,242,243,247,248,252,258,260,261,264,266,270,272,275,276,280,282,287,288,290,294,296,297,300,301,304,306,308,310,319,320,322,323,325,329,330,336,340,341,344,345,350,351,352,357,360,363,368,369,374,377,378,380,384,387,391,392,396,399,400,403,405,407,408,410,414,416,420,423,425,430,432,434,435,437,440,441,448,450,451,456,462,464,465,470,473,476,480,483,484,486,492,493,495,496,504,506,513,516,517,520,528,533,539,540,544,546,550,551,552,555,559,560,564,567,572,574,575,576,580,585,588,589,600,602,609,611,612,615,616,620,621,624,627,629,630,637,638,640,644,645,646,650,651,656,658,660,663,666,667,672,675,680,684,686,688,690,693,697,700,703,704,705,713,720,725,728,731,735,736,738,740,744,748,752,756,759,760,765,768,774,775,777,779,780,782,783,784,792,795,798,799,800,805,806,810,814,816,817,819,820,828,832,833,836,840,846,848,851,855,860,861,864,868,874,880,882,884,888,891,893,899,900,901,902,903,910,912,918,920,924,925,930,931,936,940,943,945,946,950,952,954,957,960,962,966,968,969,972,975,980,984,987,988,989,990,992,999,1000,1007,1008,1012,1015,1020,1025,1026,1029,1032,1034,1036,1040,1044,1050,1053,1054,1056,1060,1066,1073,1075,1078,1080,1081,1088,1092,1100,1102,1107,1110,1113,1116,1118,1120,1122,1125,1128,1131,1134,1144,1147,1148,1150,1152,1155,1161,1170,1175,1176,1178,1184,1188,1189,1190,1196,1200,1204,1209,1216,1218,1219,1221,1222,1224,1225,1230,1232,1242,1247,1248,1254,1258,1260,1269,1271,1275,1287,1290,1292,1295,1296,1302,1305,1312,1316,1320,1323,1325,1326,1332,1333,1344,1350,1353,1360,1363,1365,1368,1372,1375,1376,1377,1378,1380,1392,1394,1395,1400,1404,1406,1408,1410,1419,1421,1425,1426,1428,1430,1431,1435,1440,1443,1450,1452,1456,1457,1462,1470,1472,1475,1476,1480,1482,1484,1485,1488,1500,1504,1505,1508,1512,1517,1519,1520,1530,1534,1540,1548,1551,1554,1558,1560,1564,1566,1568,1581,1584,1590,1591,1593,1595,1596,1598,1599,1600,1612,1617,1620,1624,1628,1632,1634,1638,1640,1645,1650,1652,1653,1656,1664,1666,1672,1674,1677,1680,1683,1692,1696,1700,1702,1708,1710,1711,1715,1716,1720,1722,1728,1739,1740,1749,1750,1755,1760,1763,1764,1768,1769,1770,1782,1786,1792,1794,1798,1800,1802,1804,1806,1813,1820,1824,1829,1830,1836,1840,1845,1848,1850,1855,1856,1860,1862,1870,1872,1880,1886,1887,1888,1890,1891,1892,1900,1904,1908,1911,1914,1920,1927,1932,1935,1938,1944,1947,1950,1952,1960,1961,1974,1976,1978,1980,1984,1989,2000,2006,2009,2013,2014,2016,2021)

I understand that the initial input would affect the end solution but I am looking for more substantial tips on how to make more fractions.
The denominator also has to be under 2023. My method is dividing fractions up through a list using simple algorithms and dividing the ones with the least amount of fractions first and moving onto the ones with bigger fractions. I am wondering if anyone could pitch any ideas, algorithms or strategies on how to approach this and make new fractions.

Method I’m using:

Run a loop through a list of starting fractions (2,3,6) and check if they can be divided into three fractions, if not then divide into two using a*(a+b) b*(a+b) where a and b are factors that multiply to the number. It checks whether it is over the limit, already in the list or a duplicate. I then run it again and again.
I am new to coding so my code is messy, please give me tips on how to improve. Please ask me if you want me to pick out smaller pieces of code as examples or cut it down for analysis.

list = [2,3,6]
import random
lastlist = []
foul = 1
global authsuc
global possfracts
global factors

def factorise(x):
    factors1 = []
    for i in range(1, x + 1):
        if x % i == 0:
            factors1.append(i)


def make1d(x):
    global factors

    checks1 = True
    checks2 = True
    checks3 = True
    checks4 = True
    factors = []
    success = False
    for i in range(1, x + 1):
        if x % i == 0:
            factors.append(i)


    if len(factors) > 2:
        x3 = slice(3)
        factors = (factors[x3])

        xx1 = factors[0]

        xx2 = factors[1]

        xx3 = factors[2]

        total = (xx1 + xx2 + xx3)
        sum1 = (total * (x / xx1))
        sum2 = (total * (x / xx2))
        sum3 = (total * (x / xx3))

        possfracts.append(int(sum1))
        possfracts.append(int(sum2))
        possfracts.append(int(sum3))
        possfracts.remove(x)
        for i in possfracts:
            if i in lastlist:
                make1(i)


    else:
        make2(x)

def make1(x):
    global factors

    checks1 = True
    checks2 = True
    checks3 = True
    checks4 = True
    factors = []
    success = False
    for i in range(1, x + 1):
        if x % i == 0:
            factors.append(i)


    if len(factors) > 2:
        x3 = slice(3)
        factors = (factors[x3])

        xx1 = factors[0]

        xx2 = factors[1]

        xx3 = factors[2]

        total = (xx1 + xx2 + xx3)
        sum1 = (total * (x / xx1))
        sum2 = (total * (x / xx2))
        sum3 = (total * (x / xx3))

        possfracts.append(int(sum1))
        possfracts.append(int(sum2))
        possfracts.append(int(sum3))
        possfracts.remove(x)
        for i in possfracts:
            if i in lastlist or i in list:
                make1d(i)


    else:
        make2(x)

def makefracts1(x):

    global factors

    checks1 = True
    checks2 = True
    checks3 = True
    checks4 = True
    factors = []
    success = False
    for i in range(1, x + 1):
        if x % i == 0:
            factors.append(i)




    if len(factors) > 2:

        x3 = slice(3)
        factors = (factors[x3])

        xx1 = factors[0]

        xx2 = factors[1]

        xx3 = factors[2]

        total = (xx1 + xx2 + xx3)
        sum1 = (total * (x / xx1))
        sum2 = (total * (x / xx2))
        sum3 = (total * (x / xx3))

        possfracts.append(int(sum1))
        possfracts.append(int(sum2))
        possfracts.append(int(sum3))
        for i in possfracts:
            if i in lastlist or i in list:
                make1(i)

        auth(x)
    else:
        makefracts2(x)


def make2(x):
    switch = True
    while switch == True:
        f1 = random.randint(1, x)
        f2 = random.randint(1, x)
        fra1 = (f1 * (f1 + f2))
        fra2 = (f2 * (f1 + f2))
        if fra1 not in lastlist or fra2 not in lastlist or fra1 not in list or fra2 not in list: #should I remove the list checker?
            if f1 * f2 == x:

                possfracts.append(fra1)
                possfracts.append(fra2)
                possfracts.remove(x)


                switch = False
def checkdup():
    for i in possfracts:
        if possfracts.count(i) ==2:
            make1(i)
def makefracts2(x):
    switch = True

    while switch == True:
        global fra11
        global fra22
        global fra1
        global fra2

        f1 = random.randint(1,x)
        f2 = random.randint(1,x)
        if f1 * f2 == x:
            fra1 = (f1 * (f1 + f2))
            fra2 = (f2 * (f1 + f2))
            possfracts.append(fra1)
            possfracts.append(fra2)

            switch = False
    for i in possfracts:
        if i in lastlist or i in list:
            make1(i)
    auth2(x)

def auth(k):
    global possfracts
    global authsuc
    authsuc = False
    too_big = False
    alreadythere = False
    at2 = False
    at3 = False
    if len(set(possfracts)) != (len(possfracts)):
        alreadythere = True
    for i in possfracts:
        if i > 2023:
            too_big = True
    for i in possfracts:
        if i in lastlist:
            at2 = True
    for i in possfracts:
        if i in list:
            at3 = True
    if too_big == True or alreadythere == True or at2 == True or at3 == True:
        possfracts = []
        makefracts2(k)
    else:
        for i in possfracts:
            authsuc = True
            lastlist.append(i)
        print(possfracts)
        list.remove(k)
def auth2(k):
    global authsuc
    authsuc = False
    too_big = False
    alreadythere = False
    at2 = False
    at3 = False
    global foul
    foul = []
    if len(set(possfracts)) != (len(possfracts)):
        alreadythere = True
    for i in possfracts:
        if i > 2023:
            too_big = True

    for i in possfracts:
        if i in lastlist:

            at2 = True



    for i in possfracts:
        if i in list:
            at3 = True



    if alreadythere == True or too_big == True or at3 == True or at2 == True:

        lastlist.append(k)
        list.remove(k)
    else:
        for i in possfracts:

            authsuc = True
            lastlist.append(i)
        print(possfracts)
        list.remove(k)

def facts2(y):
    for i in range(1, y + 1):
        if y % i == 0:
            factors1.append(i)
def proof(x):
    x[:] = [1/i for i in x]
    print(sum(x))
def findsmall():
    for i in list:
        factors1 = []
        for x in range(1, i + 1):
            if i % x == 0:
                factors1.append(i)
        if len(factors1) < 3:
            makefracts2(i)
switch = True
while switch == True:

    factors1 = []
    possfracts = []

    if len(list) == 0:

        prooflist = []
        for i in lastlist:
            prooflist.append(i)

        print(list)
        print(len(lastlist))
        print(lastlist)
        proof(prooflist)
        count = 0
        for i in range(2024):
            amount = lastlist.count(i)

            if amount > 1:
                count = count + 1
                print(i)

        print(count)
        proof(prooflist)
        list = lastlist
        lastlist = []
        findsmall()
    pos = int(list[0])
    facts2(pos)

    if len(factors1) < 3:
        makefracts2(pos)
    else:

        makefracts1(pos)

print(len(lastlist))
print(lastlist)

count = 0
for i in range(2000):
    amount = lastlist.count(i)

    if amount > 1:
        count = count + 1
        print (i)

print(count)
proof(lastlist)

When I run the program with only the algorithm dividing the fractions into two (as mentioned previously) it generates a larger sum then using the algorithms that divide it into 3,4,5 fractions. I wonder if it is a problem with the process I am using and I am open to suggestions.

This is actually a variation of the subset sum problem (which is NP complex).

If you find the lcm between all numbers from 2 to 2020, you get a denominator (D) that you can use to represent each of the fractions as a distinct positive integer. When the sum of any subset of these integer is equal to D, you get a corresponding list of fractions that sum up to 1.

For example, if we want to represent all fractions of 1/n for value of n up to 10, we first find the lcm of 2,3,4,...10 which is 2520. Every 1/n fraction can be represented with an integer that you compute by dividing 2520 by the corresponding divisor:

1/2 = 1260 / 2520  (2520/2  = 1260)   1260   1260/2520    1/2
1/3 =  840 / 2520  (2520/3  =  840)  + 840  + 840/2520  + 1/3 
1/4 =  630 / 2520  (2520/4  =  630)
1/5 =  504 / 2520  (2520/5  =  504)
1/6 =  420 / 2520  (2520/6  =  420)  + 420  + 420/2520  + 1/6
1/7 =  360 / 2520  (2520/7  =  360)  -----  ----------  ----- 
1/8 =  315 / 2520  (2520/8  =  315)   2520   2520/2520     1
1/9 =  280 / 2520  (2520/9  =  280)
1/10=  252 / 2520  (2520/10 =  252)

So, any subset of the numbers [1260, 840, 630, 504, 420, 360, 315, 280, 252] that adds up to 2520 corresponds to a sum of the corresponding fractions that add up to 1.

Here is an example (using very naive algorithms for simplicity):

from math import gcd

def lcm(a,b): return a*b//gcd(a,b)

def subSum(A,n):
    if not A or n<1: return
    if n in A: yield [n]
    for i,a in enumerate(A):
        yield from ([a]+s for s in subSum(A[i+1:],n-a))
        
def fracSums(maxDiv):
    denom = 1
    for f in range(2,maxDiv+1):
        denom = lcm(denom,f)
    nums = [denom//d for d in range(2,maxDiv+1)]
    for ss in subSum(nums,denom):
        yield [denom//n for n in ss]

output:

for sol in fracSums(15):
    print("1 ="," + ".join(f"1/{n}" for n in sol)) 

1 = 1/2 + 1/3 + 1/6
1 = 1/2 + 1/3 + 1/10 + 1/15
1 = 1/2 + 1/4 + 1/6 + 1/12
1 = 1/2 + 1/4 + 1/10 + 1/12 + 1/15
1 = 1/3 + 1/4 + 1/6 + 1/10 + 1/12 + 1/15

You can use better subset-sum algorithms to make it more efficient, but it will remain a complex problem with an exponential time profile.

The following is more a failed attempt than an answer.

I have used the Python wrapper z3py of the Z3 solver to solve the problem.

This is my z3py model:

from z3 import *

maxDenominator = 40
merit_aim = 10
denominators = range(2, maxDenominator + 1)
solver = Solver()
a = []  #  array of unit fractions
b = []  #  array of inclusion/exclusion Booleans
for d in denominators:
    a.append(Q(1, d))  #  quotient 1/i
    b.append(Bool(f"b{d}"))

#  Constraint:
#  The sum of included unit fractions must be 1
solver.add(Sum([If(x[0], x[1], 0) for x in zip(b, a)]) == 1)

#  Constraint:
#  At least merit_aim inclusion variables must be true
solver.add(merit_aim <= Sum([If(x, 1, 0) for x in b]))

if solver.check() == sat:
    m = solver.model()

    #  count the number of inclusion variables set to true
    merit = sum([1 if is_true(m[x]) else 0 for x in b])
    print(f"fractions: {merit}")

    #  show the included unit fractions
    for d in denominators:
        if is_true(m[b[d-2]]):
            print(f"+1/{d}")
else:
    print(f"No solution with up to {merit_aim} fractions found!")

It finds the following solution in 34 seconds:

fractions: 11
+1/4
+1/8
+1/9
+1/10
+1/11
+1/12
+1/15
+1/18
+1/22
+1/24
+1/33

So, it is hopeless using this general method to beat the 570 fractions solution mentioned in the question. 2018 Boolean inclusion/exclusion variables result in a search space of 2^2018 entries (about 3 x 10^607). This is huge, unless it can be pruned by symmetry breaking or some other clever methods.

This paper "Unit fractions that sum to 1" from 2013 predicts a maximum of 62 unit fractions for denominators below 100.