通过扩展欧几里德算法在任意长度数组上查找 Bezout 系数
Finding Bezout coefficients via Extended Euclidean Algorithm on Array of Arbitrary Length
我想要什么，大图： 我想知道如何在 Java 中模仿 Mathematica 的 ExtendedGCD[...] 功能。可以找到有关该功能的信息 here，但为了完整起见，我将对其进行简要描述。
例如，在 Mathematica 中：ExtendedGCD[550,420] returns {10,{13,-17}} 因为 550 和 420 的 GCD 是 10，而 "Bezout coefficients" （源于 Bezout's Theorem) 13和17满足恒等式13(550)-17(420)=10.
这也适用于 n>2 整数：ExtendedGCD[550,420,3515] returns {5,{-4563, 5967, 1}} 因为 GCD 为 5 并且 "Bezout coefficients" -4563、5967 和 1 满足-4563(550)+5967(420)+1(3515)=5.
我目前能做的：我能计算两个整数的 GCD（无 Bezout 系数）：
public static int getGcd(int x, int y)
{
if (x == 0)
return y;
return gcd(y % x, x);
}
我可以用它来计算整数数组的 GCD（无 Bezout 系数）：
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;
}
我可以也（独立地）找到两个整数的 GCD + Bezout 系数：
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 };
}
我想要的，更准确地说：所以，总结一下我的目标：我真的想扩展方法gcdWithBezout 接受任意长度的整数数组并输出数组的 GCD 和 数组的 Bezout 系数。
问这个问题之前我做了什么：我花了很多时间试图修改我的算法。我花了很多其他时间谷歌搜索东西 + 细读 Whosebug 试图找到我能找到的东西 use/modify。我什至从计算数学+数论+任何大学课程中下载+阅读在线讲义等。
我确实做了尽职调查。
我所知道的：我知道一点点：我知道GCD(a,b,c)=GCD(a,GCD(b,c))；我知道 Bezout 的定理；我也知道如何（乱七八糟地）手工求 Bezout 系数，知道中国剩余定理等等
编辑： 我也知道 this post on math.stackexchange, as well as the very (mathematically-)insightful answer upvoted beneath it。这有助于我理解如何递归地可视化问题，但它无助于消除我似乎对使我的数据结构排列等的迷雾。
结论+精准提问：谁能帮我从我所拥有的变成我想要的？
编辑： 我给出的方法在 Java 中，我的最佳最终结果也是如此；但是，我也愿意 + 欣赏非 Java 解决方案。我可以根据自己的需要调整其他语言！
我找到了一个算法，给定输入 { 3515, 550, 420 } 产生结果：
-5*3515 + 6*550 + 34*420 = 5
该技术类似于问题中显示的 gcdWithBezout 的实现，它在递归调用堆栈向下计算余数，然后在向上计算系数。
可以使用以下算法将该技术扩展到输入数组：
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
例如，假设我们从 {3515, 550, 420} 开始。往下走，数组内容如下图。每行都有元素除以最小值后的余数。 （最小的保持不变。）
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
返回途中，系数为
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
在递归的最深层，唯一的非零元素的系数设置为1。如果我们计算乘积之和（0*0 + 1*5 + 0*0），结果是gcd。这就是目标，我们希望保持这一结果。
在第一层，系数不变。虽然数组中第一个和最后一个元素的值增加了，但它们的系数为零，所以总和没有改变。
在第二级，5变为130。差值是125。为了补偿，第一列的系数设置为-5。乘积之和现在是 -5*25 + 1*130 + 0*30 = 5.
第三层25变为155，系数-5。为了补偿，我们还需要五个 130。所以中心列中的系数从 1 变为 6。
最后155变为3515，系数为-5，同时130变为550，系数为6。3515-155相差8*420，所以我们需要8*420* 5 或 40*420 来补偿。 550-130 的差异是 420，所以我们需要 -6*420 来补偿。一共需要34*420来补偿
底线：算法保证找到解，但不一定是最优解。这个例子的最优解是{-1, -2, 11}，可以用蛮力算法找到。
这是 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);
}
注意：我想到将基本情况选择为 0 0 5 而不是 0 5 0 会得到 { -1, -2, 11 } 的最优解。因此，似乎以适当的顺序选择除数可以得到最优解。
数字是这样的：
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
BigInteger 的实现：
public static void main(String[] args) {
BigInteger[] gcdArgs = new BigInteger[] {
BigInteger.valueOf(550),
BigInteger.valueOf(420),
BigInteger.valueOf(3515) };
BigInteger[] bezoutCoefficients = new BigInteger[3];
BigInteger gcd = ExtendedGCD.extendedGCD(gcdArgs, bezoutCoefficients);
System.out.println("GCD: "+gcd.toString());
System.out.println("Bezout Coefficients: ");
for (int i = 0; i < bezoutCoefficients.length; i++) {
System.out.print(" "+bezoutCoefficients[i].toString());
}
}
/**
* Calculate the extended GCD
*
* @param gcdArgs
* an array of positive BigInteger numbers
* @param bezoutsCoefficients
* returns the Bezout Coefficients
* @return
*/
public static BigInteger extendedGCD(final BigInteger[] gcdArgs, BigInteger[] bezoutsCoefficients) {
BigInteger factor;
BigInteger gcd = gcdArgs[0];
Object[] stepResult = extendedGCD(gcdArgs[1], gcd);
gcd = (BigInteger) stepResult[0];
bezoutsCoefficients[0] = ((BigInteger[]) stepResult[1])[0];
bezoutsCoefficients[1] = ((BigInteger[]) stepResult[1])[1];
for (int i = 2; i < gcdArgs.length; i++) {
stepResult = extendedGCD(gcdArgs[i], gcd);
gcd = (BigInteger) stepResult[0];
factor = ((BigInteger[]) stepResult[1])[0];
for (int j = 0; j < i; j++) {
bezoutsCoefficients[j] = bezoutsCoefficients[j].multiply(factor);
}
bezoutsCoefficients[i] = ((BigInteger[]) stepResult[1])[1];
}
return gcd;
}
/**
* Returns the gcd of two positive numbers plus the bezout relation
*/
public static Object[] extendedGCD(BigInteger numberOne, BigInteger numberTwo) throws ArithmeticException {
Object[] results = new Object[2];
BigInteger dividend;
BigInteger divisor;
BigInteger quotient;
BigInteger remainder;
BigInteger xValue;
BigInteger yValue;
BigInteger tempValue;
BigInteger lastxValue;
BigInteger lastyValue;
BigInteger gcd = BigInteger.ONE;
BigInteger mValue = BigInteger.ONE;
BigInteger nValue = BigInteger.ONE;
boolean exchange;
remainder = BigInteger.ONE;
xValue = BigInteger.ZERO;
lastxValue = BigInteger.ONE;
yValue = BigInteger.ONE;
lastyValue = BigInteger.ZERO;
if ((!((numberOne.compareTo(BigInteger.ZERO) == 0) || (numberTwo.compareTo(BigInteger.ZERO) == 0)))
&& (((numberOne.compareTo(BigInteger.ZERO) == 1) && (numberTwo.compareTo(BigInteger.ZERO) == 1)))) {
if (numberOne.compareTo(numberTwo) == 1) {
exchange = false;
dividend = numberOne;
divisor = numberTwo;
} else {
exchange = true;
dividend = numberTwo;
divisor = numberOne;
}
BigInteger[] divisionResult = null;
while (remainder.compareTo(BigInteger.ZERO) != 0) {
divisionResult = dividend.divideAndRemainder(divisor);
quotient = divisionResult[0];
remainder = divisionResult[1];
dividend = divisor;
divisor = remainder;
tempValue = xValue;
xValue = lastxValue.subtract(quotient.multiply(xValue));
lastxValue = tempValue;
tempValue = yValue;
yValue = lastyValue.subtract(quotient.multiply(yValue));
lastyValue = tempValue;
}
gcd = dividend;
if (exchange) {
mValue = lastyValue;
nValue = lastxValue;
} else {
mValue = lastxValue;
nValue = lastyValue;
}
} else {
throw new ArithmeticException("ExtendedGCD contains wrong arguments");
}
results[0] = gcd;
BigInteger[] values = new BigInteger[2];
values[0] = nValue;
values[1] = mValue;
results[1] = values;
return results;
}
我想要什么，大图： 我想知道如何在 Java 中模仿 Mathematica 的 ExtendedGCD[...] 功能。可以找到有关该功能的信息 here，但为了完整起见，我将对其进行简要描述。
例如，在 Mathematica 中：ExtendedGCD[550,420] returns {10,{13,-17}} 因为 550 和 420 的 GCD 是 10，而 "Bezout coefficients" （源于 Bezout's Theorem) 13和17满足恒等式13(550)-17(420)=10.
这也适用于 n>2 整数：ExtendedGCD[550,420,3515] returns {5,{-4563, 5967, 1}} 因为 GCD 为 5 并且 "Bezout coefficients" -4563、5967 和 1 满足-4563(550)+5967(420)+1(3515)=5.
我目前能做的：我能计算两个整数的 GCD（无 Bezout 系数）：
public static int getGcd(int x, int y)
{
if (x == 0)
return y;
return gcd(y % x, x);
}
我可以用它来计算整数数组的 GCD（无 Bezout 系数）：
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;
}
我可以也（独立地）找到两个整数的 GCD + Bezout 系数：
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 };
}
我想要的，更准确地说：所以，总结一下我的目标：我真的想扩展方法gcdWithBezout 接受任意长度的整数数组并输出数组的 GCD 和 数组的 Bezout 系数。
问这个问题之前我做了什么：我花了很多时间试图修改我的算法。我花了很多其他时间谷歌搜索东西 + 细读 Whosebug 试图找到我能找到的东西 use/modify。我什至从计算数学+数论+任何大学课程中下载+阅读在线讲义等。
我确实做了尽职调查。
我所知道的：我知道一点点：我知道GCD(a,b,c)=GCD(a,GCD(b,c))；我知道 Bezout 的定理；我也知道如何（乱七八糟地）手工求 Bezout 系数，知道中国剩余定理等等
编辑： 我也知道 this post on math.stackexchange, as well as the very (mathematically-)insightful answer upvoted beneath it。这有助于我理解如何递归地可视化问题，但它无助于消除我似乎对使我的数据结构排列等的迷雾。
结论+精准提问：谁能帮我从我所拥有的变成我想要的？
编辑： 我给出的方法在 Java 中，我的最佳最终结果也是如此；但是，我也愿意 + 欣赏非 Java 解决方案。我可以根据自己的需要调整其他语言！
我找到了一个算法，给定输入 { 3515, 550, 420 } 产生结果：
-5*3515 + 6*550 + 34*420 = 5
该技术类似于问题中显示的 gcdWithBezout 的实现，它在递归调用堆栈向下计算余数，然后在向上计算系数。
可以使用以下算法将该技术扩展到输入数组：
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
例如，假设我们从 {3515, 550, 420} 开始。往下走，数组内容如下图。每行都有元素除以最小值后的余数。 （最小的保持不变。）
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
返回途中，系数为
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
在递归的最深层，唯一的非零元素的系数设置为1。如果我们计算乘积之和（0*0 + 1*5 + 0*0），结果是gcd。这就是目标，我们希望保持这一结果。
在第一层，系数不变。虽然数组中第一个和最后一个元素的值增加了，但它们的系数为零，所以总和没有改变。
在第二级，5变为130。差值是125。为了补偿，第一列的系数设置为-5。乘积之和现在是 -5*25 + 1*130 + 0*30 = 5.
第三层25变为155，系数-5。为了补偿，我们还需要五个 130。所以中心列中的系数从 1 变为 6。
最后155变为3515，系数为-5，同时130变为550，系数为6。3515-155相差8*420，所以我们需要8*420* 5 或 40*420 来补偿。 550-130 的差异是 420，所以我们需要 -6*420 来补偿。一共需要34*420来补偿
底线：算法保证找到解，但不一定是最优解。这个例子的最优解是{-1, -2, 11}，可以用蛮力算法找到。
这是 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);
}
注意：我想到将基本情况选择为 0 0 5 而不是 0 5 0 会得到 { -1, -2, 11 } 的最优解。因此，似乎以适当的顺序选择除数可以得到最优解。
数字是这样的：
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
BigInteger 的实现：
public static void main(String[] args) {
BigInteger[] gcdArgs = new BigInteger[] {
BigInteger.valueOf(550),
BigInteger.valueOf(420),
BigInteger.valueOf(3515) };
BigInteger[] bezoutCoefficients = new BigInteger[3];
BigInteger gcd = ExtendedGCD.extendedGCD(gcdArgs, bezoutCoefficients);
System.out.println("GCD: "+gcd.toString());
System.out.println("Bezout Coefficients: ");
for (int i = 0; i < bezoutCoefficients.length; i++) {
System.out.print(" "+bezoutCoefficients[i].toString());
}
}
/**
* Calculate the extended GCD
*
* @param gcdArgs
* an array of positive BigInteger numbers
* @param bezoutsCoefficients
* returns the Bezout Coefficients
* @return
*/
public static BigInteger extendedGCD(final BigInteger[] gcdArgs, BigInteger[] bezoutsCoefficients) {
BigInteger factor;
BigInteger gcd = gcdArgs[0];
Object[] stepResult = extendedGCD(gcdArgs[1], gcd);
gcd = (BigInteger) stepResult[0];
bezoutsCoefficients[0] = ((BigInteger[]) stepResult[1])[0];
bezoutsCoefficients[1] = ((BigInteger[]) stepResult[1])[1];
for (int i = 2; i < gcdArgs.length; i++) {
stepResult = extendedGCD(gcdArgs[i], gcd);
gcd = (BigInteger) stepResult[0];
factor = ((BigInteger[]) stepResult[1])[0];
for (int j = 0; j < i; j++) {
bezoutsCoefficients[j] = bezoutsCoefficients[j].multiply(factor);
}
bezoutsCoefficients[i] = ((BigInteger[]) stepResult[1])[1];
}
return gcd;
}
/**
* Returns the gcd of two positive numbers plus the bezout relation
*/
public static Object[] extendedGCD(BigInteger numberOne, BigInteger numberTwo) throws ArithmeticException {
Object[] results = new Object[2];
BigInteger dividend;
BigInteger divisor;
BigInteger quotient;
BigInteger remainder;
BigInteger xValue;
BigInteger yValue;
BigInteger tempValue;
BigInteger lastxValue;
BigInteger lastyValue;
BigInteger gcd = BigInteger.ONE;
BigInteger mValue = BigInteger.ONE;
BigInteger nValue = BigInteger.ONE;
boolean exchange;
remainder = BigInteger.ONE;
xValue = BigInteger.ZERO;
lastxValue = BigInteger.ONE;
yValue = BigInteger.ONE;
lastyValue = BigInteger.ZERO;
if ((!((numberOne.compareTo(BigInteger.ZERO) == 0) || (numberTwo.compareTo(BigInteger.ZERO) == 0)))
&& (((numberOne.compareTo(BigInteger.ZERO) == 1) && (numberTwo.compareTo(BigInteger.ZERO) == 1)))) {
if (numberOne.compareTo(numberTwo) == 1) {
exchange = false;
dividend = numberOne;
divisor = numberTwo;
} else {
exchange = true;
dividend = numberTwo;
divisor = numberOne;
}
BigInteger[] divisionResult = null;
while (remainder.compareTo(BigInteger.ZERO) != 0) {
divisionResult = dividend.divideAndRemainder(divisor);
quotient = divisionResult[0];
remainder = divisionResult[1];
dividend = divisor;
divisor = remainder;
tempValue = xValue;
xValue = lastxValue.subtract(quotient.multiply(xValue));
lastxValue = tempValue;
tempValue = yValue;
yValue = lastyValue.subtract(quotient.multiply(yValue));
lastyValue = tempValue;
}
gcd = dividend;
if (exchange) {
mValue = lastyValue;
nValue = lastxValue;
} else {
mValue = lastxValue;
nValue = lastyValue;
}
} else {
throw new ArithmeticException("ExtendedGCD contains wrong arguments");
}
results[0] = gcd;
BigInteger[] values = new BigInteger[2];
values[0] = nValue;
values[1] = mValue;
results[1] = values;
return results;
}