def METHOD_NAME(num_vars, num_formulas, coef=(-5, 5), bounds=(-20, 20)):
    variables = [te.var("x" + str(i)) for i in range(num_vars)]
    relations = []
    for i in range(num_formulas):
        s1 = sum([v * random.randint(coef[0], coef[1]) for v in variables])
        s1 += random.randint(coef[0], coef[1])
        s2 = sum([v * random.randint(coef[0], coef[1]) for v in variables])
        s2 += random.randint(coef[0], coef[1])
        if random.random() < 0.7:
            op = tvm.tir.EQ
        else:
            # we also make sure it can correctly handle inequalities
            op = random.choice([tvm.tir.LE, tvm.tir.LT, tvm.tir.GE, tvm.tir.GT])
        relations.append(op(s1, s2))
    vranges = {v: tvm.ir.expr.Range(bounds[0], bounds[1] + 1) for v in variables}
    solution = arith.solve_linear_equations(relations, variables, vranges)
    testing.check_int_constraints_trans_consistency(solution)
    # leaving some variables as parameters should also be ok
    for k in [1, 2]:
        if len(variables) > k:
            solution = arith.solve_linear_equations(relations, variables[:-k], vranges)
            param_ranges = {v: vranges[v] for v in variables[-k:]}
            testing.check_int_constraints_trans_consistency(solution, param_ranges)