在任意长度的数组上使用扩展欧几里得算法找到贝祖系数

Question

在任意长度的数组上使用扩展欧几里得算法找到贝祖系数

arraysalgorithmwolfram-mathematicanumber-theory

4

我希望的目标: 我想知道如何在Java中模仿Mathematica的ExtendedGCD[...]功能。有关该函数的信息可以在此处找到here，但为了完整起见，我将简要描述它。

例如，在Mathematica中：ExtendedGCD[550,420]返回{10,{13,-17}}，因为550和420的最大公约数是10，并且源自Bézout's定理的"Bezout系数"13和17满足恒等式13(550)-17(420)=10。

这对于n>2的整数也适用：ExtendedGCD[550,420,3515]返回{5,{-4563, 5967, 1}}，因为GCD为5，而"Bezout系数" -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); 
}

我可以使用这个方法计算一个整数数组的最大公约数（不包括 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;
}

我可以独立地找到两个整数的最大公约数和贝祖等式系数：

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，使其能够接受任意长度的整数数组，并输出该数组的最大公约数以及 Bezout 系数。

在提出问题之前我做了什么： 我花了很多时间尝试修改算法。我还花了很多时间在 Google 和 StackOverflow 上搜索和浏览课程讲义等东西，以寻找可以使用/修改的内容。我甚至下载并阅读了计算数学、数字论等课程的在线讲义等等。

我肯定已经尽了自己的努力。

我知道的： 我了解一点点：我知道 GCD(a,b,c)=GCD(a,GCD(b,c))；我知道贝祖定理；我也知道如何（笨拙地）手动求解贝祖系数，以及中国剩余定理等。

编辑： 我还知道关于这篇文章以及下面得到赞同的非常（数学上）有洞见的答案。这些帮助我理解如何递归地可视化问题，但它并没有帮助我澄清我所遇到的数据结构不匹配等问题。

结论和具体问题： 有人能帮我从现有的代码到达我想要的结果吗？

编辑： 我提供的方法是使用Java编写的，我的最终目标也是如此；然而，我也接受并欣赏非Java语言的解决方案。我可以将其他语言适应我的需求！

- cstover

请注意，{5，{-4563，5967，1}}是一种解决方案，但它不是最佳解决方案。例如，更好的解决方案是{5，{-2，11，-1}}。递归应用ExtendedGCD的问题在于贝祖定理系数变得非常大。 - user3386109

@user3386109 - 感谢您的评论。有两个回复：（a）对于我的初始目的，系数的最优性并不重要。（b）是否有一种方法[数学或编程]来隔离“最优”贝祖系数？ - cstover

我找到了一个相当有效的解决方案，提供了一个不错但不是最优的系数集。但我不懂Java，所以代码是用C编写的。如果您想看它，请删除Java标签，以便我不会因在Java问题上发布C答案而受到谴责。 - user3386109

@user3386109 - 我已经移除了Java标签，并对问题进行了明确的编辑，这为任何/所有语言/伪代码解决方案打开了大门。感谢您的考虑和分享意愿！ - cstover

我已经发布了一个答案。我会离开一段时间，但稍后会回来看看你是否有任何问题。 - user3386109

@user3386109 - 非常感谢您详细的帖子！我也会离开一段时间（为我的手术后随访+我妻子的化疗），但是只要我能够深入研究您的答案，我将尝试通过提出问题来跟进。再次感谢！ - cstover

2个回答

1

一个关于 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;
    }

- axelclk

非常感谢您的回答！就使用BigInteger而言，它非常有见地，只需进行最小程度的调整即可直接奏效。然而，我最终投票支持@user3386109的答案，因为它给了我更多关于如何思考问题，如何努力优化比宇系数等方面的直觉。再次感谢您的工作！ - cstover

不用谢。这只是我项目的一小部分，可以在 https://github.com/axkr/symja_android_library 找到。 - axelclk

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user3386109 · Accepted Answer

我能找到一个算法，输入为 { 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以-5的系数变成了3515，同时130以6的系数变成了550。差值为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