How can I represent a fraction as a finite sum of reciprocal squares?

Question

I've found this result : $$ \frac{1}2 = \frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\frac1{5^2}+\frac1{6^2}+\frac1{15^2}+\frac1{18^2}+\frac1{36^2}+\frac1{60^2}+\frac1{180^2}$$ I've tried my best to make a program to find that result but i can't find a solution. I found this solution to an easier problem : the representation of a fraction as the sum of distinct reciprocals. I tried to expend this algorithm to solve my problem but when i execute my program the numbers given get too large:

def pgcd(a,b):
    return a if b==0 else pgcd(b,a%b)
def subFrac(p1,q1,p2,q2):
    u,v=p1*q2-p2*q1,q1*q2
    d=pgcd(u,v)
    return (u//d,v//d)
def dec(p,q):
    print(p,q)
    frac=[]
    k=2
    while p>0:
        while 1/(k**2)>p/q:
            k+=1
        p,q=subFrac(p,q,1,k**2)
        print(k,p,q)
        frac.append(k)
        k+=1
    return frac

Gives : $$ \frac{1}2 = \frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\frac1{5^2}+\frac1{6^2}+\frac1{11^2}+\frac1{54^2}+\frac1{519^2}+\frac1{59429^2}+...$$ In general, I'd like to make a program that gives a representation with reasonable number like the one given or this one from proofWiki: $$ \frac{1}2 = \frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\frac1{5^2}+\frac1{7^2}+\frac1{12^2}+\frac1{15^2}+\frac1{20^2}+\frac1{28^2}+\frac1{35^2}$$ Thanks in advance!

EDIT:

Thanks for your answers. I looked into the problem 152 of Euler project and I saw that they mostly brute force to find solutions and try to reduce the initial list size using math. I will probably try this later. Here is a link to a solution (not mine) written in C++ and here in python.

I've tried James Arathoon solution and it works well. Here is the full python code if you are interested ( to find combinations i used a method which works but i think using a recursive function like they do in most project Euler 152 solutions is better ):

import numpy.polynomial.polynomial as nppol # To make the variable you called list
# Useful function, i could use already existing module
def Divisors(k):
    l=[]
    for ii in range(1,k+1):
        if k%ii==0:
            l.append(ii)
    return l
def pgcd(a,b):
    return a if b==0 else pgcd(b,a%b)
def addFrac(p1,q1,p2,q2):
    assert q1!=0 and q2!=0
    u,v=p1*q2+p2*q1,q1*q2
    d=pgcd(u,v)
    return (u//d,v//d)
def subFrac(p1,q1,p2,q2):
    u,v=p1*q2-p2*q1,q1*q2
    d=pgcd(u,v)
    return (u//d,v//d)
# The interesting part
def test(alist, frac):
    for i in alist:
        testDivisors=Divisors(i)
        testDivisorsSq=[k*k for k in testDivisors]

        prod=[1] # Its a polynomial
        for j in testDivisorsSq:
            prod=nppol.polymul(prod,[1]+(j-1)*[0]+[1])
        llist=[k for k in range(len(prod)) if prod[k]!=0]

        totalInvTestDivisorsSq=(0,1)
        for k in testDivisorsSq[1:]:    # Remove 1 which is always the first element of TestDivisorsSq
            totalInvTestDivisorsSq=addFrac(*totalInvTestDivisorsSq,1,k)
        totalInvTestDivisorsSq=subFrac(*totalInvTestDivisorsSq,*frac)
        
        testnum=testDivisorsSq[-1]/totalInvTestDivisorsSq[1]*totalInvTestDivisorsSq[0]
        if testnum in llist:
            print("Possible Solutions : ")
            print(i,"  ",testDivisorsSq,"  ",testnum) # Max(Divisors(i))=i always   and no need to reverse the list
            # Find all the correct combination that add up to testnum
            l=len(testDivisorsSq)
            for j in range(1,2**l):
                auth=("0"*l+bin(j)[2:])[-l:]
                s=0
                k=0
                while s<=testnum and k<l:
                    if auth[k]=="1":
                        s+=testDivisorsSq[k]
                    k+=1
                if s==testnum:
                    print("Find a combination : ",[testDivisorsSq[k] for k in range(1,l) if auth[l-1-k]=="0"])
# its l-1-k because its reversed
test([1,2,6,12,18,30,36,60,90,120,150,180],(1,2))
# It took a few seconds to run this program

It gives :

Possible Solutions : 
120    [1, 4, 9, 16, 25, 36, 64, 100, 144, 225, 400, 576, 900, 1600, 3600, 14400]    500.0
Find a combination :  [4, 9, 16, 25, 64, 100, 225, 400, 576, 900, 1600, 3600, 14400]
Find a combination :  [4, 9, 16, 25, 64, 100, 144, 576, 900, 1600, 3600, 14400]

Possible Solutions : 
180    [1, 4, 9, 16, 25, 36, 81, 100, 144, 225, 324, 400, 900, 1296, 2025, 3600, 8100, 32400]    1086.0
Find a combination :  [4, 9, 16, 25, 81, 100, 144, 225, 400, 8100, 32400]
Find a combination :  [4, 9, 16, 25, 36, 225, 324, 1296, 3600, 32400]
Find a combination :  [4, 9, 16, 25, 36, 225, 400, 1296, 2025, 3600, 8100]
Find a combination :  [4, 9, 16, 25, 36, 144, 1296, 2025, 3600, 8100]
Find a combination :  [4, 9, 16, 25, 81, 100, 144, 324, 400, 900, 3600, 8100]
Find a combination :  [4, 9, 16, 25, 81, 100, 144, 225, 900, 1296, 2025, 3600]

Answer 1

This is a concept for a basic algorithm (without proof and optimization). Its quick and dirty and incomplete; the final stage involves a manual search operation, which can be automated with lots of extra effort.

Start with a list of the Heinz Numbers: that is $\{1,2,6,12,18,30,36,60,90,120,150,180,210,270,300,360,420,450,540,600,630,750,840,900\}$

For the desired summation of an Egyptian Fraction, $1/2$ say, the algorithm searches within the list of Heinz Numbers provided. Solutions for a $1/2$ are only found in the above list, for the numbers $120,180,300,360,420,540,600,840$ and $900$

The algorithm spits out the Heinz Number used e.g $120$, the list of divisors of the squared Heinz Number (e.g. $120^2=14400$) in reverse order and lastly a number which is equal the sum of one or more subsets of the divisors squared. For the Heinz Number of $120$ it is $500$. This is testnum in the Mathematica code.

The code basically searches for testnum in a list of all the possible subset sums of the divisors squared for that Heinz Number. If it finds a match, there exists one or more solutions. The number of elements in the list is always smaller than the Heinz Number squared, but nevertheless will slow the algorithm as the size size of Heinz number grows larger.

$120\text{ }\{14400,3600,1600,900,576,400,225,144,100,64,36,25,16,9,4,1\}\text{ }500$

Obviously one way summing to $500$ is $400+100$. That is elements $6$ and $9$ in the above list. Then strike out elements $6$ and $9$ form the list below (i.e. $1/36$ and $1/144$), plus the first one (being 1), and we have a solution.

$\left\{1,\frac{1}{4},\frac{1}{9},\frac{1}{16},\frac{1}{25},\frac{1}{36},\frac{1}{64},\frac{1}{100},\frac{1}{144},\frac{1}{225},\frac{1}{400},\frac{1}{576},\frac{1}{900},\frac{1}{1600},\frac{1}{3600},\frac{1}{14400}\right\}$

For the Heinz Number of $180$ the testnum sum of squares is $1086$, with some familiar solutions resulting.

$180\text{ }\{32400,8100,3600,2025,1296,900,400,324,225,144,100,81,36,25,16,9,4,1\}\text{ }1086$

$\left\{1,\frac{1}{4},\frac{1}{9},\frac{1}{16},\frac{1}{25},\frac{1}{36},\frac{1}{81},\frac{1}{100},\frac{1}{144},\frac{1}{225},\frac{1}{324},\frac{1}{400},\frac{1}{900},\frac{1}{1296},\frac{1}{2025},\frac{1}{3600},\frac{1}{8100},\frac{1}{32400}\right\}$

Mathematica Code used:

Inputs:

alist - is the list of Heinz Numbers

frac - is the desired fraction e.g. $1/2$

test[alist_, frac_] := (
  Do[testDivisors = Divisors[i]; testDivisorsSq = Divisors[i]^2;
   testDivisorsSqTotal = Total[testDivisorsSq];
   Clear[prod];
   prod = 1;
   Do[prod = prod*(1 + y^j);, {j, testDivisorsSq}];
   list = Exponent[Expand[prod], y, List];
   invtestDivisorsSq = 1/testDivisorsSq;
   totalinvtestDivisorsSq = Total[invtestDivisorsSq] - 1;
   testnum = (Last[testDivisorsSq]/
       Denominator[totalinvtestDivisorsSq - frac])*
     Numerator[totalinvtestDivisorsSq - frac];
   If[MemberQ[list, testnum] == True, 
    Print[Max[testDivisors], "  ", Reverse[testDivisorsSq], "  ", 
     testnum];
    Print[invtestDivisorsSq];
    ];
   , {i, alist}])