#Calculate A072368 using integer linear programming #Based on http://list.seqfan.eu/pipermail/seqfan/2023-July/074777.html from mip import Model, xsum, minimize, BINARY, OptimizationStatus #https://www.python-mip.com/ import functools import itertools import operator def product(iterable): return functools.reduce(operator.mul, iterable, 1) for n in itertools.count(1): model = Model() #One binary variable for each triple var = [(pqr, model.add_var(var_type=BINARY)) for pqr in itertools.combinations(range(1, 3*n+1), 3)] model.objective = minimize(xsum(product(pqr) * x for pqr, x in var)) #Constraints: each number 1..3n must be used exactly once for i in range(1, 3*n+1): model += xsum(x for pqr, x in var if i in pqr) == 1 model.max_mip_gap = 0 model.verbose = 0 status = model.optimize() assert status == OptimizationStatus.OPTIMAL print(f'{n} {model.objective_value:.0f}') example = (pqr for pqr, v in var if v.x >= 1/2) #print(' + '.join(['*'.join([str(x) for x in pqr]) for pqr in example]))