如何将递归函数转换为 mex 代码？

Question

我在MATLAB代码中有一个递归函数选择如下：

   function nk=choose(n, k)
        if (k == 0)
            nk=1;
        else
            nk=(n * choose(n - 1, k - 1)) / k;
        end
    end

该代码用于计算 n 和 k 之间的组合。我想通过使用 mex 代码来加速它。我尝试将 mex 代码编写为

double choose(double* n, double* k)
{
   if (k==0) 
        return 1;
   else
        return (n * choose(n - 1, k - 1)) / k;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *n, *k, *nk;
    int mrows, ncols;
    plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
    /* Assign pointers to each input and output. */
    n = mxGetPr(prhs[0]);    
    k = mxGetPr(prhs[1]);
    nk = mxGetPr(plhs[0]);
    /* Call the recursive function. */
    nk=choose(n,k);
}

但是，它不起作用。你能帮我修改一下可以实现上述MATLAB代码的mex代码吗？谢谢

Answer 1

不需要mex二项式系数可以在Matlab中实现：

function nk=nchoosek2(n, k)
    if n-k > k
        nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))));
    else
        nk = prod((n-k+1:n) .* prod((1:k).^ (-1/k)) ) ;
    end
end

Answer 2

以下代码修复了您的 C mex 实现。
问题当然不是递归...
您的代码使用指针而不是值（在 C 中，仅在正确的位置使用指针很重要）。

您可以使用 Matlab 内置函数：nchoosek
参见：http://www.mathworks.com/help/matlab/ref/nchoosek.html

以下代码有效：

//choose.c

#include "mex.h"

double choose(double n, double k)
{
    if (k==0) 
    {
        return 1;
    }
    else
    {
        return (n * choose(n - 1, k - 1)) / k;
    }
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *n, *k, *nk;
    int mrows, ncols;
    plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
    /* Assign pointers to each input and output. */
    n = mxGetPr(prhs[0]);    
    k = mxGetPr(prhs[1]);
    nk = mxGetPr(plhs[0]);

    /* Call the recursive function. */
    //nk=choose(n,k);
    *nk = choose(*n, *k);
}

在 Matlab 中编译它： mex choose.c

执行：
choose(10,5)
ans =
</code><br> <code>252

这不是低效的实施...
我正在帮助修复您的实现，用作 "inefficient example".

测量 rahnema1 实施的执行情况：
tic;n = 1000000;k = 500000;nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))));toc
经过的时间是 0.022855 秒。

测量 choose.mexw64 实现的执行情况：
tic;n = 1000000;k = 500000;nk = 选择 (1000000, 500000);toc
经过的时间是 0.007952 秒。
（比 prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k)))) 花费的时间少一点）。

测量 Matlab 递归，出现错误（即使对于 n=700 和 k=500）：
ic;n = 700;k = 500;nk = RecursiveFunctionTest(n, k);toc
达到最大递归限制 500。使用 set(0,'RecursionLimit',N) 更改限制。请注意，超过您的可用堆栈 space 可以使 MATLAB and/or 计算机崩溃。

tic;n = 700;k = 400;nk = RecursiveFunctionTest(n, k);toc
经过的时间是 0.005635 秒。效率很低...

测量 Matlab 内置函数 nchoosek:
tic;nchoosek(1000000, 500000);toc
警告：结果可能不准确。系数大于9.007199e+15且只精确到15位在 nchoosek 92 经过的时间是 0.005081 秒。

结论：
需要在不使用递归的情况下实现C mex文件，采取措施

无递归测量：

static double factorial(double number) 
{
    int x;
    double fac = 1;

    if (number == 0)
    {
        return 1.0;
    }    

    for (x = 2; x <= (int)number; x++)
    {
        fac = fac * x;
    }

    return fac;
}



double choose(double n, double k)
{
    if (k == 0) 
    {
        return 1.0;
    }
    else
    {
        //n!/((n–k)! k!) 
        return factorial(n)/(factorial(n-k)*factorial(k));
    }
}

tic;choose(1000000, 500000);toc

Elapsed time is 0.003079 seconds.

更快...

如何将递归函数转换为 mex 代码？

How to convert a recursive function to mex code?

matlab

mex