2

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

我要找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))。这对每个人都一样吗？

4

2 回答 2

2

我在 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
于 2014-10-14T06:15:29.437 回答
1

在 2013a 中，这可以复制...

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

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

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

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

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

在 2013a 中，这些输入binCoef不被称为

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

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

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

可以做什么？

我可以看到两个选项；

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

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

选项一 -应该相当直截了当......

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

于 2014-10-14T10:49:16.063 回答