Analyse https://mathoverflow.net/q/29478 by integer relation search with LLL

MO29478.sage

def integer_relation(vals_to_relate, precision = 200):
    n = len(vals_to_relate)
    m = matrix(n, n+1)
    for i in range(n): m[i,i] = 1
    for i in range(n):
        m[i,n] = int(vals_to_relate[i]*RealField(precision)(2)^(precision-20))//2^10

    L = m.LLL()
    return L.row(0)


def ordered_sequences(max_val, k):
    if k == 0:
        yield []
        return

    for last_val in range(1, max_val + 1):
        for prefix in ordered_sequences(last_val, k-1):
            yield prefix + [last_val]


def analyse_MO29478(n, k):
    def term_m(m): return sin(pi * m / (n*k + 1))
    def a(i): return 1 if i == 0 else a(i-1) * term_m(k-2+i) / term_m(n+k-1-i)

    bases = set()
    for seqlen in range(1, 8):
        for indices in ordered_sequences(n-1, seqlen):
            subterm = 1
            for idx in indices:
                subterm *= a(idx)

            # Original problem as stated in question: product/quotient of 1 - prod_{i in S} a_i for some S
            if log(1 - subterm) not in bases:
                print("-", indices, float(log(1 - subterm)))
                bases.add(log(1 - subterm))

            # Expansion implied in comment: admit 1 + prod_{i in S}
            if log(1 + subterm) not in bases:
                print("+", indices, float(log(1 + subterm)))
                bases.add(log(1 + subterm))

    print("")
    prev_Aj = dict()
    for j in range(1, n):
        Aj = log(1 - term_m(k-1) * term_m(n+k-1) / term_m(k-1+j) / term_m(n+k-1-j))
        if Aj in prev_Aj:
            print("A_" + str(j) + " = A_" + str(prev_Aj[Aj]))
            continue

        prev_Aj[Aj] = j
        print("A_" + str(j), float(Aj))

        # We seek an integer relation giving Aj in terms of the bases (or, more precisely, the logarithms of both).
        # However, the bases are not independent, so we loop removing dependent bases until we get a relation involving Aj.
        vals_to_relate = [Aj] + list(bases)
        while vals_to_relate:
            ir = integer_relation(vals_to_relate)
            print(ir)
            if ir[0] != 0: break

            killidx = max(i for i in range(1, len(ir) - 1) if ir[i])
            print("  Killing index", killidx)
            del vals_to_relate[killidx]
        print("")


#analyse_MO29478(3, 3) # Sanity-check a known good case

#analyse_MO29478(4, 2) # Finds some plausible relations
analyse_MO29478(4, 3) # Not really for A_1

#analyse_MO29478(5, 2) # Finds some plausible relations
#analyse_MO29478(5, 3) # Finds some plausible relations - TODO are they the same ones?
#analyse_MO29478(5, 4) # Not really for A_1 or A_2