Matlab nchoosek 使用 int64 和 sym 得到不同的答案

Question

这是一个关于Matlab中的函数nchoosek的问题。

我要找nchoosek(54,25)，和54C25一样。由于答案大约是 10^15，所以我最初使用int64。然而，就象征性而言，答案是错误的。

输入：

nchoosek(int64(54),int64(25))
nchoosek(sym(54),sym(25))

输出：

1683191473897753
1683191473897752

您可以看到它们相差一个。这不是一个真正紧迫的问题，因为我现在使用sym。但是有人能告诉我为什么会这样吗？

---

编辑：

我正在使用 R2013a。

我看了一下nchoosek.m，发现如果输入在int64，代码可以简化成

function c = nchoosek2(v,k)

    n = v;  % rename v to be n. the algorithm is more readable this way.

    classOut = 'int64';
    nd = double(n);
    kd = double(k);
    nums = (nd-kd+1):nd;
    dens = 1:kd;
    nums = nums./dens;      %%
    c = round(prod(nums));
    c = cast(c,classOut);
end

但是，int64(prod(nums./dens)) 的结果对我来说与prod(sym(nums)./sym(dens)) 不同。每个人都一样吗？

Answer 1

我在 R2014a 上没有这个问题：

数字

>> n = int64(54);
>> k = int64(25);
>> nchoosek(n,k)
ans =
     1683191473897752    % class(ans) == int64

象征性

>> nn = sym(n);
>> kk = sym(k);
>> nchoosek(nn,kk)
ans =
1683191473897752         % class(ans) == sym

% N!/((N-K)! K!)
>> factorial(nn) / (factorial(nn-kk) * factorial(kk))
ans =
1683191473897752         % class(ans) == sym

---

如果您查看函数 edit nchoosek.m 的源代码，您会发现它使用单独的算法专门处理 64 位整数的情况。我不会在这里复制代码，但这里是重点：

function c = nchoosek(v,k)
    ...

    if int64type
        % For 64-bit integers, use an algorithm that avoids
        % converting to doubles
        c = binCoef(n,k,classOut);
    else
        % Do the computation in doubles.
        ...
    end

    ....
end

function c = binCoef(n,k,classOut)
    % For integers, compute N!/((N-K)! K!) using prime factor cancellations
    ...
end

Answer 2

在 2013a 中可以复制...

@Amro 在nchoosek 中显示了一个特殊情况，用于 int64 或 unit64 的 classOut，
但是在 2013a 中，这仅适用于答案介于两者之间的情况

• flintmax（没有参数）和
• double(intmax(classOut)) + 2*eps(double(intmax(classOut)))

int64 给出 9007199254740992 和 9223372036854775808，解决方案不在...

---

如果解决方案落在这些值之间，它将使用子函数 binCoef 重新计算 帮助说明：对于整数，使用素因数取消计算 N!/((N-K)!M!)

binCoef 函数会为给定的 int64 输入生成正确答案

在 2013a 中，这些输入 binCoef 未被调用

取而代之的是“默认”帕斯卡三角法，其中：

• 输入转换为双精度
• 取向量((n-k+1):n)./(1:k)的乘积
• 此向量包含k 分数的双重表示。

所以我们几乎可以肯定浮点错误。

---

可以做什么？

我可以看到两个选项；

1. 根据binCoef中的代码制作自己的函数，
2. 修改 nchoosek 并从第 81 行删除 && c >= flintmax

删除此表达式将强制 Matlab 对 int64 和 uint64 的输入使用更精确的基于整数的计算，以获取其精度范围内的任何值。这会稍微慢一些，但可以避免浮点错误，这在处理整数类型时是理所当然的。

选项一 - 应该相当直接...

选项二 - 我建议保留原始函数的未更改备份，或者制作带有修改的函数的副本并改用它。