Solving a Linear Diophantine Equation(see description for examples)

Question

Let me start off by clarifying that(before you guys dismiss me), this is not a homework problem and I'm not a university student. :)

EDIT Thanks to @Klas and others, my question now boils down to a mathematical equation which needs to be solved programmatically.

I'm looking for an algorithm/code which solves Linear Diophantine Equation. For lesser mortals like me, here's how such an equation looks like:

Example 1: 3x + 4y + 5z = 25 (find all possible values of x,y,z)

Example 2: 10p + 5q + 6r + 11s = 224 (find all possible values of p,q,r,s)

Example 3: 8p + 9q + 10r + 11s + 12t = 1012 (find all possible values of p,q,r,s,t)

I tried googling to no avail. I would have thought that some code would already be written to solve this. Do let me know if you guys have come across some kind of library which has already implemented this. And if the solution is in Java, nothing can be cooler!. Algorithm/pseudo code will also do. Thanks much.

Vieira

I'm sorry about my bad math terminologies, haven't done from long. I'm trying generate a question paper randomly, based on certain constraints(which are complicated and unnecessary for others to know). I've tried to make this problem independent and simplified as much as possible. – Vieira

Voted to close; not programming related. Should be on something like math.stackexchange.com – Sleep

I'm looking to programmatically solve this problem. And after Klas' answer, I'm looking for code which solves Diophantine Equations. It's definitely programming related IMHO. – Vieira

Answer 1

Brute-force recursion is an option, depending on how large you will allow the value or number of values to become.

Assumptions: The user inputs (the multiplicands) are always distinct positive integers. The coefficients to be found must be non-negative integers.

Algorithm:

Of the multiplicands, let M be the largest.
Calculate C=floor(F/M).
If F=M*C, output solution of the form (0,0,...,C) and decrement C
If M is the only multiplicand, terminate processing
Loop from C down to 0 (call that value i)
  Let F' = F - i*M
  Recursively invoke this algorithm:
    The multiplicands will be the current set minus M
    The goal value will be F'
  For each solution output by the recursive call:
     append i to the coefficient list and output the solution

Statvolt

Also please check the following link for the implementation of the above algorithm - #5513629 – Vieira

This is worse. for any Algorithm scaling should and must be taken into account. Imagine how many programmers and software developers with the same problem that apply this exponential algorithm to their projects. Its worse. for anyone reading this please first do your research on linear diophantine, implementation although complicated is worth it. its under Integer Programming Algebra. – Sjambok

Its crazy that in Math stack exchange there is a far better answer that require unimodular row operations on the 'A' matrix or something. but the asker has being downvoted, so it might be hard to find the question. you only need to implement after understanding the Math. – Sjambok

Answer 2

This is a mathematical question rather than a programming one. Once you have a suitable algorithm, impelementing it shouldn't be too hard.

I suggest you google on Diophantine Equations.

I found an explanation for you.

Cassiecassil

I think the question is about implementing a solution in Java, so this is a programming question. Implementing a solution for any number of possible coefficients seem to be an interesting and relevant coding problem. A link would be more useful than saying "go google it". E.g. mathworld.wolfram.com/DiophantineEquation.html – Junkie

@Steve: This is fundamentally a math question; there's nothing Java-specific about the question (the algorithm used, whatever it may be, is the same regardless of the language). – Sleep

From the FAQ "if your question generally covers … a software algorithm … then you’re in the right place". It might not be java specific, but Q. is askign for a software algorithm so it is on topic. – Junkie

Guys, would there be an already existing code to solve Linear Diophantine Equations in some library? I tried googling to no avail. – Vieira

Answer 3

I happened to write Java code for this. Please help yourself. The solutions are not extensively tested, but it seems to work well so far.

package expt.qp;

import java.util.HashMap;
import java.util.LinkedHashMap;
import java.util.Map;

public class LinearDiophantine {

    private Map<Integer, Integer> sol = new LinkedHashMap<Integer, Integer>();
    private Map<Integer, Integer> coeff = new HashMap<Integer, Integer>();

    /**
     * @param args
     */
    public static void main(String[] args) {
        // Fill up the data
        // 3x + 4y + 5z + 3a = 25
        LinearDiophantine ld = new LinearDiophantine();
        ld.coeff.put(1, 1);ld.coeff.put(2, 2);ld.coeff.put(3, 3);ld.coeff.put(4, 4);
        Map<Integer, Integer> coeffCopy = new HashMap<Integer, Integer>(ld.coeff);
        int total=30;

        // Real algo begins here
        ld.findPossibleSolutions(total, coeffCopy);

    }

    private void findPossibleSolutions(int total, Map<Integer, Integer> coeff) {
        int index=returnLargestIndex(coeff);
        int range = (int) Math.floor(total/coeff.get(index));
        if(range*coeff.get(index) == total) {
            sol.put(index, range);
            displaySolution();
            //System.out.println();
            range--;
        }
        if(coeff.size() == 1) {
            return;
        }
        while(range>=0) {
            int remTotal = total - range*coeff.get(index);
            Map<Integer, Integer> coeffCopy = new HashMap<Integer, Integer>(coeff);
            coeffCopy.remove(index);
            sol.put(index, range);
            findPossibleSolutions(remTotal, coeffCopy);
            range--;
        }
    }

    private void displaySolution() {
        int total = 0;
        for(int i : sol.keySet()) {
            //System.out.print(coeff.get(i)+"("+sol.get(i)+"), ");
            total = total + (coeff.get(i)*sol.get(i));
        }
        if(total != 30)
            System.out.print(total+",");
    }

    /**
     * @param coeff
     */
    private int returnLargestIndex(Map<Integer, Integer> coeff) {
        int largestKey = coeff.keySet().iterator().next();
        for(int i : coeff.keySet()) {
            if(coeff.get(i)>coeff.get(largestKey)) {
                largestKey=i;
            }
        }
        return largestKey;
    }

}

Vieira

The solution is based on @Dave Costa's (accepted) answer for the question. – Vieira

Just running the above example seems to show some issues with the code/algorithm. Many totals != 30 are printed, implying some kind of issue. – Traylor

Yes, I did make some minor changes to the above code. Unfortunately can't remember where, nor do I have access to that code now. – Vieira

Answer 4

Adding on to Klas' very accurate answer:

Hilbert's 10th problem asked if an algorithm existed for determining whether an arbitrary Diophantine equation has a solution. Such an algorithm does exist for the solution of first-order Diophantine equations. However, the impossibility of obtaining a general solution was proven by Yuri Matiyasevich in 1970

taken from: Wolfram MathWorld

Angeliqueangelis

I think he is talking about non-linear Diophantine equations. For linear Diophantine equations it should be possible. – Vieira

Answer 5

A brute force algorithm is as follows (3 variable case):

int sum = 25;
int a1 = 3;
int a2 = 4;
int a3 = 5;
for (int i = 0; i * a1 <= sum; i++) {
    for (int j = 0; i * a1 + j * a2 <= sum; j++) {
        for (int k = 0; i * a1 + j * a2 + k * a3 <= sum; k++) {
            if (i * a1 + j * a2 + k * a3 == sum) {
                System.out.println(i + "," + j + "," + k);
            }
        }
    }
}

To generalize this for the N variable case, you need to convert into a recursive form.

This algorithm is O(f(size, a)^N) for some function f.

• We can place bounds on f as follows: size / MaxValue(a) <= f(size, a) <= size / MinValue(a).
• In the worst case (where all of the a[i]s are 1) f(size, a) is size.

Either way, this is pretty horrendous for large values of N. So while the recursive N variable algorithm would be more elegant, it is probably not very practical.

---

If you are willing to fork out 34 Euro's to Springer Verlag, here's a link to a paper which (according to the abstract) includes an algorithm for solving the N variable case.

Guevara

Thanks for solution @Stephen, I will try to see if there some libraries existing already. Do let me know if you find any. – Vieira

Why can't i, j and/or k be negative? – Stichous

If they can be negative, there can be an infinite number of solutions. (Indeed, if N > 1 and there is any solution, there will be an infinite number of solutions. There's a simple construction that proves this.) – Guevara

@Vieira - Does your program included a survey of previous literature on the subject, a formal proof of correctness, a formal complexity analysis and benchmarks? – Guevara

Answer 6

There are either no, or infinitely many solutions. It is often the case that you have an extra constraint that the solution must match. Is this the case in your problem?

Let's start with the most simple situation where there are two unkowns a*x + b*y = c:

The first step is using the Euclidean algorithm to find the GCD of a and b, let's call itd. As a bonus, the algorithm provides x' and y' such that a*x' + b*y' = d. If d doesn't divide c, then there is no solution. Otherwise, a solution is:

x = x' * (c/d)
y = y' * (c/d)

The second step is to find all solutions. This means we must find all p and q such that a*p + b*q = 0. For if both (x,y) and (X, Y) are solutions, then

a * (X-x) + b * (Y-y) = 0

The answer to this is p = b/d and q = -a/d where d = GCD(a,b) and is already calculated in step 1. The general solution now is:

x = x' * (c/d) + n * (b/d)
y = y' * (c/d) - n * (a/d)

where n is an integer.

The first step is easy to extend to multiple variables. I am not sure about generalizing the second step. My first guess would be to find a solution for all pairs of coefficients and combine these solutions.

Pinsky

Answer 7

This is a solution in Perl. rather a hack by using Regex.

Following this blog post to solve algebraic equations using regex.

we can use the following script for 3x + 2y + 5z = 40

#!/usr/bin/perl
$_ = 'o' x 40;
$a = 'o' x 3;
$b = 'o' x 2;
$c = 'o' x 5;
$_ =~ /^((?:$a)+)((?:$b)+)((?:$c)+)$/;
print "x = ", length($1)/length($a), "\n";
print "y = ", length($2)/length($b), "\n";
print "z = ", length($3)/length($c), "\n";

output: x=11, y = 1, z = 1

the famous Oldest plays the piano puzzle ends up as a 3 variable equation

This method applies for a condition that the variables are actually positive and the constant is positive.

Gurl

And how does this scale? I believe large values will be very slow – Cowgirl

Answer 8

The reason why there is no library that does this is similar to why you won't find a library to do permutations - you generate so much data that this is probably the wrong thing to do.

More precisely, if you have n variables whose sum totals X, then you will have O(Xn-1) answers. X and n don't have to be very large for this to become an issue.

That said, here is some Python that fairly efficiently figures out all of the necessary information to encode the answer:

def solve_linear_diophantine (*coeff_tuple):
    coeff = list(coeff_tuple)
    constant = coeff.pop()

    cached = []
    for i in range(len(coeff)):
        # Add another cache.
        cached.append({})

    def solve_final (i, remaining_constant):
        if remaining_constant in cached[i]:
            return cached[i][remaining_constant]
        answer = []
        if i+1 == len(coeff):
            if 0 == remaining_constant%coeff[i]:
                answer = [(remaining_constant/coeff[i], [])]
        else:
            for j in range(remaining_constant/coeff[i] + 1):
                finish = solve_final(i+1, remaining_constant - j*coeff[i])
                if finish is not None:
                    answer.append((j, finish))
        if 0 == len(answer):
            cached[i][remaining_constant] = None
            return None
        else:
            cached[i][remaining_constant] = answer
            return answer

    return solve_final(0, constant)

When I say "encode", the data structure looks like this. For each possible coefficient, we'll get an array of (coefficient, [subanswers]). Whenever possible it reuses subanswers to make the data structure as small as possible. If you cant you can recursively expand this back out into an array of answers, and in so doing you'll be pretty efficient. (In fact it is O(nX).) You'll do very little repeating of logic to discover the same facts over and over again. (However the returned list may get very large simply because there is a large list of answers to be returned.)

Hungary