Walnut 7 performance: Frobenius Coin Problem - Walnut-Theorem-Prover/Walnut GitHub Wiki

(For a quick refresher on this problem, see the Wikipedia entry)

At the bottom of this page is a Python script that generates a set of Frobenius Coin problems for Walnut, with a given number of variables and a given numerical bound on its constants. An example execution (below) generates 3 problems, each with 2 variables, and with an upper bound of 350.

python3 generate_fcp.py 2 3 350
def test1 "~ E x1 E x2 (137*x1 + 303*x2 = P) & (A R ((~ E x1 E x2 (137*x1 + 303*x2 = R)) =>(R <= P)))":
def test2 "~ E x1 E x2 (217*x1 + 333*x2 = P) & (A R ((~ E x1 E x2 (217*x1 + 333*x2 = R)) =>(R <= P)))":
def test3 "~ E x1 E x2 (59*x1 + 244*x2 = P) & (A R ((~ E x1 E x2 (59*x1 + 244*x2 = R)) =>(R <= P)))":

Running 100 trials, we find Walnut 7 takes a median of 1 second, with a max of 12 seconds. This compares very favorably to other, similar tools. A 2024 paper on Amaya and benchmarks indicate that similar calculations (2 variables, constants increasing up to 350 ?) take a median of 3-4 seconds in Amaya (with 10% timing out, i.e., taking over 60 seconds), and 6 seconds in Lash (with 20% timing out, i.e., taking over 60 seconds). The SMT-based tools (Z3, cvc5, PRINCESS) showed much worse performance.

As we increase the number of variables and the bounds, these Frobenius problems become significantly slower. 3 variables and a bound of 1000 yields

def test1 "~ E x1 E x2 E x3 (202*x1 + 264*x2 + 327*x3 = P) & (A R ((~ E x1 E x2 E x3 (202*x1 + 264*x2 + 327*x3 = R)) =>(R <= P)))":

which takes closer to 7 seconds. 4 variables and a bound of 1000 yields

def test1 "~ E x1 E x2 E x3 E x4 (4*x1 + 27*x2 + 189*x3 + 712*x4 = P) & (A R ((~ E x1 E x2 E x3 E x4 (4*x1 + 27*x2 + 189*x3 + 712*x4 = R)) =>(R <= P)))":

which takes closer to 22 seconds.

These results do not compare favorably to special-purpose Frobenius solvers, such as within Wolfram Alpha. But for automata-based tools, this appears quite good.

#!/usr/bin/env python3
import random
import sys
import math

def generate_coefficients(num_vars, lower_bound=1, upper_bound=100):
    # Ensure there is a large enough range to pick distinct numbers.
    if upper_bound - lower_bound + 1 < num_vars:
        raise ValueError("Range too small to generate the required number of distinct integers.")
   
    while True:
        # Pick distinct integers and sort them.
        coeffs = random.sample(range(lower_bound, upper_bound + 1), num_vars)
        coeffs.sort()  # smallest integer corresponds to x1, etc.
        # Compute the GCD of all coefficients.
        current_gcd = coeffs[0]
        for c in coeffs[1:]:
            current_gcd = math.gcd(current_gcd, c)
        if current_gcd == 1:
            return coeffs

def generate_formula(num_vars, test_index, lower_bound=1, upper_bound=100):
    # Generate coefficients with the desired properties.
    coefficients = generate_coefficients(num_vars, lower_bound, upper_bound)
    # Create variable names: x1, x2, ..., xN.
    var_names = [f"x{i+1}" for i in range(num_vars)]
    # Build the summation string, e.g., "2*x1 + 3*x2 + 5*x3"
    sum_str = " + ".join(f"{coef}*{var}" for coef, var in zip(coefficients, var_names))
   
    # Create the existential quantifier string as "E x1 E x2 ...".
    existential = " ".join(f"E {var}" for var in var_names) + f" ({sum_str} = P)"
    existential_R = " ".join(f"E {var}" for var in var_names) + f" ({sum_str} = R)"
   
    # Build the final formula string.
    formula = (
        f'def test{test_index} "~ {existential} & (A R ((~ {existential_R}) =>(R <= P)))":'
    )
    return formula

def main():
    if len(sys.argv) < 4:
        print("Usage: python3 script.py num_variables num_tests upper_bound")
        sys.exit(1)
   
    try:
        num_vars = int(sys.argv[1])
        num_tests = int(sys.argv[2])
        upper_bound = int(sys.argv[3])
    except ValueError:
        print("Error: Both arguments must be integers.")
        sys.exit(1)
   
    for i in range(1, num_tests + 1):
        print(generate_formula(num_vars, i, upper_bound=upper_bound))

if __name__ == "__main__":
    main()