C++ Diophantine equation solver with constraints

diophantine.cpp

/*
    This is a quick and dirty implementation of a diophantine equation solver.
    Given a vector of coefficient, a vector of constraints and a results value, the solver returns a vector of solutions.
    
    Let's take this equation as an example:
        40*x1 + 296*x2 + 945*x3 + 2048*x4 + 4500*x5 + 8640*x6 = 616103

    We have the following constraints: x1, x2, x3, x4, x5 and x6 shall be between 29 and 95.

    The following instructions will compute all the possible solutions:
        auto res = solve({40,296,945,2048,4500,8640}, {std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95)}, 616103);

    To compile, it should be linked with lgmp and boost needs to be installed:
        g++ -std=c++17 -o diophantine diophantine.cpp -lgmp
*/
#include <iostream>
#include <vector>
#include <algorithm>
#include <tuple>
#include <utility>

#include <boost/multiprecision/gmp.hpp>

typedef boost::multiprecision::mpz_int bigint;

/*
    Return the gcd of n1 and n2.
    The gcd is always positive.
    The gcd of 0 and 0 is 0.
    The gcd of 0 and n is abs(n).
*/
bigint gcd(bigint n1, bigint n2) {
    if(n1 == 0 && n2 == 0)
        return 0;
    else if(n1 == 0)
        return abs(n2);
    else if(n2 == 0)
        return abs(n1);

    n1 = abs(n1);
    n2 = abs(n2);

    bigint reminder = std::min(n1, n2);
    bigint quotient = std::max(n1, n2);

    while(reminder != 0) {
        bigint tmp_reminder = quotient % reminder;
        quotient = reminder;
        reminder = tmp_reminder;
    }

    return quotient;
}

/*
    Return the gcd of a list of integers.
*/
bigint gcd(std::vector<bigint> v) {
    v.erase(std::remove_if(v.begin(), v.end(), [] (auto n) { return n == 0; }));

    if(v.empty())
        return 0;

    if(v.size() == 1)
        return v[0];

    bigint current_gcd = gcd(v[0], v[1]);

    for(std::size_t i = 2; i < v.size(); i++)
        current_gcd = gcd(current_gcd, v[i]);

    return current_gcd;
}

/*
    Return the sign of x.
    If x is positive (or zero), return 1.
    Otherwise return -1.
*/
bigint sign(bigint x) {
    if(x >= 0)
        return 1;

    return -1;
}

/*
    Solve the diophantine equation: x1*a + x2*b = y
    The return value is a tuple of 4 numbers:
     1. a_0: particular solution for a
     2. b_0: particular solution for b
     3. a_n: step to compute the general solution for a
     4. b_n: step to compute the general solution for b

    The general solution can compute that way, with k an integer:
        a = a_0 + a_n * k
        b = b_0 + b_n * k
        
*/
std::tuple<bigint, bigint, bigint, bigint> solve(bigint x1, bigint x2, bigint y) {
    bigint gcd_x = gcd(x1, x2);

    if(gcd_x == 0)
        throw("Can't solve 0a + 0b = c equations");

    if(y % gcd_x != 0)
        throw(std::string() + "No solutions for " + std::string(x1) + "a + " + std::string(x2) + "b = " + std::string(y));

    x1 /= gcd_x;
    x2 /= gcd_x;
    y /= gcd_x;

    if(std::min(abs(x1), abs(x2)) == 0) {
        if(abs(x1) == 0)
            return std::make_tuple(0, y, 0, 0);
        else
            return std::make_tuple(y, 0, 0, 0);
    }

    if(std::min(abs(x1), abs(x2)) == 1) {
        if(abs(x1) == 1)
            return std::make_tuple(-x1*(x2-1)*y, y, x1*x2, -1);
        else
            return std::make_tuple(y, -x2*(x1-1)*y, -1, x1*x2);
    }

    std::vector<bigint> quotients;
    bigint numerator = std::max(abs(x1), abs(x2));
    bigint reminder = std::min(abs(x1), abs(x2));

    while(reminder != 1) {
        quotients.push_back(numerator / reminder);
        bigint tmp_reminder = numerator % reminder;
        numerator = reminder;
        reminder = tmp_reminder;
    }

    bigint sol_max = 0;
    bigint sol_min = 1;

    for(auto it = quotients.rbegin(); it != quotients.rend(); it++) {
        bigint tmp_sol_min = sol_min;
        sol_min = (sol_min * (-(*it))) + sol_max;
        sol_max = tmp_sol_min;
    }

    if(x1 > x2)
        return std::make_tuple(sign(x1) * sol_max * y, sign(x2) * sol_min * y, sign(x1)*x2, -sign(x2)*x1);
    else
        return std::make_tuple(sign(x1) * sol_min * y, sign(x2) * sol_max * y, sign(x1)*x2, -sign(x2)*x1);
    
}

/*
    Find all the solutions for the given parameters, constraints and results.
    The return value is a vector with a size that corresponds to the number of solution satisfaying the constraints.
*/
std::vector<std::vector<bigint>> solve(std::vector<bigint> x_vec, std::vector<std::pair<bigint, bigint>> constraints, bigint y) {
    if(x_vec.size() < 2)
        throw("Not enough parameters");

    bigint x_gcd = gcd(x_vec);

    if(x_gcd == 0)
        throw("Unable to solve");

    if(y % x_gcd != 0)
        throw("No solutions");
    
    if(x_vec.size() == 2) {
        auto results = solve(x_vec[0], x_vec[1], y);
        bigint x_0 = std::get<0>(results);
        bigint y_0 = std::get<1>(results);
        bigint xn = std::get<2>(results);
        bigint yn = std::get<3>(results);

        std::vector<std::vector<bigint>> return_value;

        for(int n = 0; x_0 - n * xn * sign(xn) >= constraints[0].first; n++) {
            bigint current_x = x_0 - n * xn * sign(xn);
            bigint current_y = y_0 - n * yn * sign(xn);

            if(current_x <= constraints[0].second && current_y >= constraints[1].first && current_y <= constraints[1].second)
                return_value.push_back(std::vector<bigint>({current_x, current_y}));
        }

        for(int n = 1; x_0 + n * xn * sign(xn) <= constraints[0].second; n++) {
            bigint current_x = x_0 - n * xn * sign(xn);
            bigint current_y = y_0 - n * yn * sign(xn);

            if(current_x >= constraints[0].first && current_y >= constraints[1].first && current_y <= constraints[1].second)
                return_value.push_back(std::vector<bigint>({current_x, current_y}));
        }

        return return_value;
    } else {
        bigint x2 = x_vec.back();
        x_vec.pop_back();
        bigint x1_gcd = gcd(x_vec);

        auto results = solve(x1_gcd, x2, y);
        bigint w_0 = std::get<0>(results);
        bigint x_0 = std::get<1>(results);
        bigint wn = std::get<2>(results);
        bigint xn = std::get<3>(results);

        std::vector<std::vector<bigint>> return_value;

        auto current_constraint = constraints.back();
        constraints.pop_back();

        for(int n = 0; x_0 - n * xn * sign(xn) >= current_constraint.first; n++) {
            bigint current_x = x_0 - n * xn * sign(xn);
            bigint current_w = w_0 - n * wn * sign(xn);

            if(current_x <= current_constraint.second) {
                auto current_result = solve(x_vec, constraints, current_w * x1_gcd);

                for(auto r : current_result) {
                    r.push_back(current_x);
                    return_value.push_back(r);
                }
            }
        }

        for(int n = 1; x_0 + n * xn * sign(xn) <= current_constraint.second; n++) {
            bigint current_x = x_0 - n * xn * sign(xn);
            bigint current_w = w_0 - n * wn * sign(xn);

            if(current_x >= current_constraint.first) {
                auto current_result = solve(x_vec, constraints, current_w * x1_gcd);

                for(auto r : current_result) {
                    r.push_back(current_x);
                    return_value.push_back(r);
                }
            }
        }

        return return_value;
    }
}

int main() {
    auto res = solve({40,296,945,2048,4500,8640}, {std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95),std::make_pair(29,95)}, 616103);
    
    for(auto e : res) {
        for(auto v : e) {
            std::cout << v << " ";
        }
        std::cout << "\n";
    }

    return 0;
}