#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]))