What I want, big picture: I want to know how to mimic Mathematica's ExtendedGCD[...] functionality in Java. Info about that function can be found here, but I'll describe it briefly for completeness.
For example, in Mathematica: ExtendedGCD[550,420] returns {10,{13,-17}} because the GCD of 550 and 420 is 10, and the "Bezout coefficients" (stemming from Bezout's Theorem) 13 and 17 satisfy the identity 13(550)-17(420)=10.
This also works for n>2 integers: ExtendedGCD[550,420,3515] returns {5,{-4563, 5967, 1}} because the GCD is 5 and "Bezout coefficients" -4563, 5967, and 1 satisfy -4563(550)+5967(420)+1(3515)=5.
What I can do currently: I can compute the GCD (sans Bezout coefficients) of two integers:
public static int getGcd(int x, int y)
{
if (x == 0)
return y;
return gcd(y % x, x);
}
I can use this to compute the GCD (sans Bezout coefficients) of an array integers:
public static int findGCD(int arr[])
{
int gcd = arr[0];
for (int i = 1; i < arr.length; i++) {
gcd = getGcd(arr[i], gcd);
}
return gcd;
}
I can also (independently) find the GCD + Bezout Coefficients of two integers:
public static int[] gcdWithBezout(int p, int q) {
if (q == 0)
return new int[] { p, 1, 0 };
int[] vals = gcdWithBezout(q, p % q);
int d = vals[0];
int a = vals[2];
int b = vals[1] - (p / q) * vals[2];
return new int[] { d, a, b };
}
What I want, more precisely: So, to summarize my goal: I would really like to extend the method gcdWithBezout to accept an integer array of arbitrary length and to output the GCD of the array and the Bezout coefficients of the array.
What I did before asking this question: I spent lots of time attempting to modify my algorithms. I spent lots of other time Googling things + perusing StackOverflow in attempts to find things I could use/modify. I even went so far as to download + read online lecture notes, etc. from University courses in computational mathematics + number theory + whatever.
I definitely did my due diligence.
What I know: I know a little bit about a little bit: I know that GCD(a,b,c)=GCD(a,GCD(b,c)); I know Bezout's theorem; I also know how to (messily) do the process of finding Bezout coefficients by-hand, and know about the Chinese Remainder Theorem, etc.
Edit: I also know about this post on math.stackexchange, as well as the very (mathematically-)insightful answer upvoted beneath it. This helps me understand how to visualize the problem recursively, but it doesn't help clear up the fog I seem to have about making my data structures line up, etc.
Conclusion + Precise Question: Can anyone help me get from what I have to what I want?
Edit: The methods I've given are in Java and my optimal end result will be as well; however, I am open to + appreciative of non-Java solutions as well. I can adapt other languages to my needs!
I was able to find an algorithm that given the input { 3515, 550, 420 } produces the result:
-5*3515 + 6*550 + 34*420 = 5
The technique is similar to the implementation of gcdWithBezout shown in the question, which computes remainders on the way down the recursive call stack, and then computes the coefficients on the way up.
That technique can be extended to an array of inputs with the following algorithm:
find the smallest non-zero element of the array
compute the remainders for all other elements using the smallest element as the divisor
recursively call the function until only one non-zero element remains
that element is the gcd
the coefficient for that element starts at 1
all other coefficients start at 0
on the way back up the call stack, update the coefficient of the smallest element as follows
for each element of the array that was reduced
compute the difference between the original value and the reduced value
divide by the smallest element (the difference is guaranteed to be a multiple of the smallest value)
multiply by the coefficient for the reduced element
subtract the result from the coefficient of the smallest element
For example, let's say that we start with {3515, 550, 420}. On the way down, the array contents are as shown below. Each line has the remainders after dividing the elements by the smallest. (The smallest remains unchanged.)
3515 550 420 // top level
155 130 420 // third level
25 130 30 // second level
25 5 5 // first level
0 5 0 // base case, the remaining non-zero element is the gcd
On the way back up, the coefficients are
0 1 0 // base case
0 1 0 // first level up
-5 1 0 // second level up
-5 6 0 // third level up
-5 6 34 // top level
At the deepest level of recursion, the coefficient of the only non-zero element is set to 1. If we compute the sum of products (0*0 + 1*5 + 0*0), the result is the gcd. That's the goal, and we want to maintain that result.
At the first level up, the coefficients don't change. Although the first and last elements in the array increase in value, their coefficients are zero, so the sum doesn't change.
At the second level, 5 changes to 130. The difference is 125. To compensate, the coefficient in the first column is set to -5. The sum of products is now -5*25 + 1*130 + 0*30 = 5.
At the third level, 25 changes to 155, with a coefficient of -5. To compensate, we need five more 130s. So the coefficient in the center column changes from 1 to 6.
Finally, 155 changes to 3515 with a coefficient of -5, and at the same time 130 changes to 550 with a coefficient of 6. The difference 3515-155 is 8*420, so we need 8*420*5 or 40*420 to compensate. The difference 550-130 is 420 so we need -6*420 to compensate. In total we need 34*420 to compensate.
Bottom line: the algorithm is guaranteed to find a solution, but not necessarily the optimal solution. The optimal solution for this example is { -1, -2, 11 }, which can be found with a brute force algorithm.
Here's a sample implementation in C:
int bezout3(int oldarray[], int coeff[], int length)
{
// the array must have at least 2 elements
assert(length >= 2);
// make a copy of the array
int array[length];
memcpy(array, oldarray, sizeof(array));
// find the smallest non-zero element of the array
int smallest = 0;
for (int i = 1; i < length; i++)
{
if (array[i] != 0 && (array[smallest] == 0 || array[i] < array[smallest]))
smallest = i;
}
// all elements must be non-negative, and at least one must be non-zero
assert(array[smallest] > 0);
// for every element of the array, compute the remainder when dividing by the smallest
int changed = false;
for (int i = 0; i < length; i++)
if (i != smallest && array[i] != 0)
{
array[i] = array[i] % array[smallest];
changed = true;
}
// base case: the array has only one non-zero element, which is the gcd
if (!changed)
{
coeff[smallest] = 1;
return array[smallest];
}
// recursively reduce the elements of the array till there's only one
int result = bezout3(array, coeff, length);
// update the coefficient of the smallest element
int total = coeff[smallest];
for (int i = 0; i < length; i++)
if (oldarray[i] > array[i])
{
int q = oldarray[i] / array[smallest];
total -= q * coeff[i];
}
coeff[smallest] = total;
// return the gcd
return result;
}
int main(void)
{
int array[3] = { 3515, 550, 420 };
int coeff[3] = { 0, 0, 0 };
int length = 3;
int result = bezout3(array, coeff, length);
printf("gcd=%d\n", result);
int total = 0;
for (int i = 0; i < length; i++)
{
printf("%4d * %4d\n", array[i], coeff[i]);
total += array[i] * coeff[i];
}
printf("total=%d\n", total);
}
Note: I occurs to me that choosing the base case as 0 0 5 instead of 0 5 0 results in the optimal solution of { -1, -2, 11 }. So it seems that choosing the divisor in the proper order can result in an optimal solution.
The numbers look like this:
3515 550 420
155 130 420
25 130 30
25 5 5
0 0 5
0 0 1
0 0 1
-1 0 1
-1 -2 1
-1 -2 11