algorithm - 用功能语言了解此多项式除法算法

Question

**这是一种以相当古老的语言标准ML实现的算法。我很难理解此算法：

fun polyquotremd ts ((n,b) :: us) = 
    let fun quo [] qs = (rev qs, [])
        | quo ((m,a) :: ts) qs = 
            if m < n then (rev qs, (m,a) :: ts)
            else
                quo (polysum ts
                    (map (termproduct(m-n, ~a/b)) us))
                    ((m-n, a/b) :: qs)
    in
        quo ts []
    end;

这是目标。我们通过使用元组列表来密集表示多项式，每个元组包含索引为0的幂和索引为1的系数。例如，$ x ^ 2 + 1 $表示为 [(2, 1.0), (0, 1.0)]。这由类型 (int*real)list给出，如下面的辅助函数所示。我们想要获得一个输出 (quo, rem)，每个输出都是分别表示商和余数的元组列表。 ts是要除的多项式，而 us是除数。

这是我想象的算法的工作原理：


首先，如果 ts = []意味着我们没有多项式可除，那么我们返回qs，即我们的商数结果作为答案。由于我们一直在qs列表的开头添加 ((m-n, a/b) :: qs)，因此我们需要反转列表以使列表以最大的指数开头。
第一个if条件处理除数n的指数大于我们正在处理的当前幂m的情况。在这种情况下，我们只是不添加元素并返回 (rev qs, (m,a) :: ts)，即 (quo, rem)
其他条件是我感到困惑的部分。我知道第二个括号 ((m-n, a/b) :: qs)是我们要添加的累积结果，但是 polysum ts (map (termproduct(m-n, ~a/b)) us))实际在做什么？这背后的数学是什么？


通常，我们会将多项式除数立即除以一起，即，将xx + 1一起除。但是此算法首先使用$ x $进行划分，然后使用$ 1 $进行划分。这怎么工作？为什么会有 ~a/b（供您参考，表示 -a/b）。



辅助功能：

fun take ([], _) = []
| take (x::xs, i) = if i>0
then x :: take(xs, i-1)
else [];

fun drop ([], _) = []
| drop (x::xs, i) = if i>0 then drop(xs,i-1)
else x::xs;

fun termproduct (m,a) (n,b) = (m+n, a*b) : (int*real);

fun polyproduct [] us = []
    | polyproduct [(m,a)] us = map (termproduct(m,a)) us
    | polyproduct ts us =
    let 
        val k = length ts div 2
    in 
        polysum (polyproduct (take(ts,k)) us)
                (polyproduct (drop(ts,k)) us)
    end;

有人可以向我解释第3点吗？我对此感到非常困惑，尤其是对该算法的数学运算感到困惑。但是它可以工作。仅供参考，您可以编译此代码，并可以在标准ML中对其进行测试。

Answer 1

我认为您当前的理解大部分是正确的，而且我不确定您在else分支中到底不了解什么。

基本思想和不变性

为了避免混淆名称阴影，让ts0为传递给ts的polyquotremd的值，并且us0 = ((n,b) :: us)。然后对于递归的每个步骤中从外部上下文捕获的quo ts qs的us0，除了最后一个，即由返回quo [] qs = (rev qs, [])和if m < n then (rev qs, (m,a) :: ts)分支处理的qs，以下不变式成立

ts0 = us0 * rev qs + ts

换句话说，我们迭代地（通过递归）在答案（ ts）中产生下一项，而新的 ts具有较低的最高幂，表示该除法运算的其余部分。显然，如果我们可以通过此过程降低 m的最高功率，使其最高功率 us0小于 n最高功率 qs，那么我们在 ts和 ts中有完整的答案。因此，现在我们要做的就是降低 ((m,a) :: ts)的最大功耗。我们将得到的只是标准的 long division算法。

长除法算法细节

长除法的第一步如下：


对于股息中的最高项（当前余数） ((n,b) :: us)和除数 (m,a)，即对于 (n,b)和 m < n，检查是否为 (m,a)。如果是，请完成算法。
将 (n,b)除以 (m-n, -a/b)并找到它们的商。显然是 (map (termproduct(m-n, a/b)) ((n,b) :: us)))
将整个除数乘以该商，即 map。请注意，这不是您的代码所在的行。请稍等，我们也会到达那里。


旁注或如何在代码中完成。万一这还不清楚： (termproduct(m-n, a/b))是一个经典的高阶函数，它接受另一个函数和一个集合，并返回相同大小的新集合，其中每个元素都是对该函数应用相应结果的结果源元素。现在 termproduct称为 partial application：它接受两个参数 (m-n, a/b)的函数，并产生一个新参数的函数，其中一个参数（原始函数中的第二个）的第一个参数固定为 ((n,b) :: us)的值。因此，整个构造如下：将 us的每一项（或代码中的 (m-n, a/b)）乘以 -1。显然，这只是将多项式乘以一个项的一种方法。


从当前剩余部分中减去该产品。我们可以通过在上一步中乘以 (m-n, a/b)来用总和代替该减法。所以现在是

(polysum ((m,a) :: ts)
                (map (termproduct(m-n, ~a/b)) ((n,b) :: us)))

现在我们可以注意到，最高条款实际上会相互抵消。不足为奇，因为我们将乘数计算为其商，即正好使它​​们彼此抵消。现在，我们可以删除它们，并将其替换为代码中包含的那一行：

(polysum ts
                (map (termproduct(m-n, ~a/b)) us))

最后，我们需要将 rev qs的商/乘数添加到我们的结果多项式中。唯一的次要问题是我们以相反的顺序填写它，即我们首先找到最高的术语，然后找到下一个，依此类推。在纸上可以正常工作，但ML中的列表以相反的顺序填充。这就是为什么我们的两个“出口”分支都在其中包含 quo的原因。


如果您想举一个例子-您可以在我上面引用的 Polynomial long division上的Wiki文章中关注一个例子。只是重申一下：该算法的每个步骤都是所有 termproduct (m-n, ~a/b) (n,b)的另一递归。

希望这可以帮助。让我知道是否仍然不清楚。

更新（评论的答案）


我还没有完全理解〜a / b在数学上如何相互抵消。


首先，您是否尝试遵循我上面引用的 example in the wiki？如果是这样，您在那里还不清楚什么呢？

回到您的问题，有几种方法可以解决这个问题。从纯机械角度来看，这里的问题是： (m-n + n, ~a/b*b)的结果是什么？显然，它是 (m, -a) = ts。因此，这与 (m-n, ~a/b)中的最高术语完全相反。

从更高级别的算法出发，它们将完全抵消，因为我们计算乘数 2*x^2 + 3*x + 4的唯一目的是抵消了它们。看来您没有得到我在“基本概念和不变式”部分中描述的内容，而且我不确定是否可以用其他方式重述，但我会尽力而为。除法可以看作是多次相减，而我们却做了几次。如果我们这样看，很显然，我们应该尝试减去除数的一些乘数，以抵消除数的最高项。

示例1：假设我们将 x^2 + x + 1除以 2*x^2。如果要删除 2项，则必须将除数乘以 2，然后减去。因此结果是商= 2*x^2 + 3*x + 4 - 2*(x^2 + x + 1)，其余为 x + 2 = 2*x^2 + 3*x + 4。

示例2：假设我们将 x^2 + 2*x + 3除以 2。请注意，尽管除数不同于之前的示例，但我们仍然必须将其乘以 2才能去除股息中的最高项。这是因为确定乘数的唯一因素是最高期限的比率。其他不合理的术语会影响最终结果（其余），但不会影响当前结果。请注意，结果仍然是商= 2*x^2 + 3*x + 4 - 2*(x^2 + 2*x + 3)，但其余部分是 -x - 2 = 2*x^3+3*x^2

例3：让我们将 x+1除以 x+2和 x^3来观察差异。请注意，要取消 2*x^3/x，我们的第一个乘数（即商中的最高项）在两种情况下均应相同 2*x^2 = x+1
所以对于 x+2

x+1: 2*x^3+3*x^2 = 2*x^2*(x+1) + (2*x^3+3*x^2) - 2*x^2*(x+1) = 2*x^2*(x+1) + x^2

对于 +x^2同样：

x+2: 2*x^3+3*x^2 = 2*x^2*(x+2) + (2*x^3+3*x^2) - 2*x^2*(x+2) = 2*x^2*(x+2) - x^2

从现在开始，结果将有所不同，因为我们继续对 x+1使用不同的中间休息 -x^2，对于 x+2使用不同的中间休息。这意味着下一个乘数将是 +x^2/x = +x和 -x^2/x = -x

x+1: 2*x^3+3*x^2 = (2*x^2+x)*(x+1) + x^2 - x*(x+1) = (2*x^2+x)*(x+1) - x
x+2: 2*x^3+3*x^2 = (2*x^2-x)*(x+2) - x^2 - (-x)*(x+2) = (2*x^2-x)*(x+2) + 2*x

现在，对 -x =>乘数其余为 x+1的最后一步为 -1，对于 2*x =>乘数为 x+2的 2

x+1: 2*x^3+3*x^2 = (2*x^2+x-1)*(x+1) - x - (-1)(x+1) = (2*x^2+x-1)*(x+1) + 1
x+2: 2*x^3+3*x^2 = (2*x^2-x+2)*(x+2) + 2*x - (2)*(x+2) = (2*x^2-x+2)*(x+2) - 4

同样，您可以看到除数的非最高项的更改会影响结果，但不会影响结果的最高项。