from collections import defaultdict
from itertools import combinations
from math import sqrt, floor

###############
def sigma(n):
    m = n
    facts = {}
    top = int(sqrt(m))+1
    for i in prlist:
        while m%i == 0:
            if i in facts:
                facts[i] += 1
            else:
                facts[i] = 1
            m = m/i
            top = int(sqrt(m))+1
            if i >= top:
                break
        if i >= top:
            break
    #up to sqrt not found - m is a prime too!
    if m in facts:
        facts[m] += 1
    else:
        facts[m] = 1
    
    sig = 1
    for k in facts:
        part = 1
        for i in xrange(1,facts[k]+1):
            part += k**i
        sig = sig * part
    return sig

###############
adict = defaultdict(list)   #dictionary of Dsum - [numbers]
print "sigma, x, y, n - where sigma(n) = sigma(x) = sigma(y) = 2n+x+y, n != x, n != y, x < y."
prlist = [2,3,5,7]
for i in xrange(9,500000,2  ):      #good to 250B
    top = int(sqrt(i)) + 1
    addit = True
    for j in xrange(3,top,2):
        if i % j == 0:
            addit = False
            break
    if addit:
        prlist.append(i)
n = 2
try:
    while n < 10000000000:
        if n % 50000 == 0:
            #print for progress
            print 20*"\b",
            print "{:,}".format(n),
            if n % 200000 == 0:
                #clean out the dictionary for memory and speed
                for i in adict.keys():
                    try:
                        tester = i - (2*adict[i][0]+adict[i][1])
                    except:
                        tester = i - (3*adict[i][0])
                    if tester < n:
                        del adict[i]
                    else:
                        while i - (2*adict[i][0] + adict[i][len(adict[i])-1]) < n:
                            adict[i].pop()
        x = sigma(n)
        if x > 4*n:
            adict[x].append(n)
        elif x in adict:
            for j in combinations(adict[x], 2):
                found = False
                if x == 2*j[0] + j[1] + n:
                    amin = str(min(j[1], n))
                    amax = str(max(j[1], n))
                    amiline = str(x) + " " + amin + " " + amax + " " + str(j[0])
                    found = True
                elif x == j[0] + 2*j[1] + n:
                    amin = str(min(j[0], n))
                    amax = str(max(j[0], n))
                    amiline = str(x) + " " + amin + " " + amax + " " + str(j[1]) 
                    found = True
                elif x == j[0] + j[1] + 2*n:
                    amin = str(min(j[0], j[1]))
                    amax = str(max(j[0], j[1]))
                    amiline = str(x) + " " + amin + " " + amax+ " " + str(n) 
                    found = True
                if found:
                    print 20*"\b",
                    print amiline
            adict[x].append(n)
        n += 1
        
    smallest = n
    for i in adict.keys():
        if min(adict[i]) < smallest:
            smallest = min(adict[i])
    print "smallest kept is " + str(smallest)
except:
    smallest = n
    for i in adict.keys():
        if min(adict[i]) < smallest:
            smallest = min(adict[i])
    print "smallest kept is " + str(smallest)
    raise