Calling from Matlab a recusive function written in C

Hi,

I have a written a recursive function in C to replace nchoosek() in Matlab. When i call this function in Matlab, it gives error if i use output parameter. If i don't use output parameter then it does not return the required answer. I have tried to debug it but of no use.

Code without output parameter is:

#include "math.h"

#include "mex.h"

double test(double n, double k)
{

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

void mexFunction(int nlhs, mxArray *plhs[], int nrhs,
        const mxArray *prhs[])
{
    double y, n, k;

      /* Create matrix for the return argument. */
      plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
      /* Assign pointers to each input and output. */
      n = mxGetScalar(prhs[0]);
      k = mxGetScalar(prhs[1]);
       /* Call the test subroutine. */
      test(n,k);
  }

Code with output parameter is:

#include "math.h"

#include "mex.h"

double test(double *y, double n, double k)
{

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

void mexFunction(int nlhs, mxArray *plhs[], int nrhs,const mxArray *prhs[])
{

      double *y, n, k;

      /* Create matrix for the return argument. */
      plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
      /* Assign pointers to each input and output. */
      n = mxGetScalar(prhs[0]);
      k = mxGetScalar(prhs[1]);
      y = mxGetPr(plhs[0]);
      /* Call the test subroutine. */
      test(y,n,k);
  }

I have to implement this function since nchoosek(n,k) gives warning if i compute it for higher values of n & k.

I will be highly obliged if any one is kind enough to help me out.

Ameer

Answer 1

[James Tursa]

In addition to what Geoff has already written, in your "code without output parameter" you are not setting the value of the return variable. E.g., these lines:

      plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
      /* Assign pointers to each input and output. */
      n = mxGetScalar(prhs[0]);
      k = mxGetScalar(prhs[1]);
       /* Call the test subroutine. */
      test(n,k); // <-- Calculates a value then throws it away, doesn't get into plhs[0]

should look something like this instead:

      /* Assign pointers to each input and output. */
      n = mxGetScalar(prhs[0]);
      k = mxGetScalar(prhs[1]);
       /* Call the test subroutine. */
      plhs[0] = mxCreateDoubleScalar(test(n,k)); // <-- Calculates value and gets it into plhs[0]

As a side note, it would be much better to code this up in a loop inside mexFunction rather than chewing up a lot of stack space and overhead on all of those recursive function calls. The loop will be cleaner, use less resources, and certainly run faster.

Also, if you are getting warnings or errors from nchoosek because the numbers you are using are too large, then a mex routine isn't likely to fix that. You may need to rethink how you are solving your problem.

[Ameer Ahmed]

 I have done it like

plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL); /* Assign pointers to each input and output. */

      n = mxGetScalar(prhs[0]);

      k = mxGetScalar(prhs[1]);

      y = mxGetPr(plhs[0]);

      /* Call the test subroutine. */

      test(y,n,k);

But it still gives error. I don't know how to fix it.

[Ameer Ahmed]

Thanks James

I just added

plhs[0] = mxCreateDoubleScalar(test(n,k)); into my code and i got the result. Thanks for your advice.

[James Tursa]

E.g., here is a bare bones loop implementation. This is for example purposes only. At the very least, much more argument checking would have to be inserted in order for it to be "production quality" code.

/* File:  nchoosekx.c                                                    */
/* Programmer: James Tursa                                               */
/* Bare bones C-mex replacement of nchoosek (for example purposes)       */
/* Includes ------------------------------------------------------------ */

#include "mex.h"

/* Gateway ------------------------------------------------------------- */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double i, n, k, m, result;
    // Need to add more input checking here (complex, size, etc)
    if( nrhs == 2 && mxIsNumeric(prhs[0]) && mxIsNumeric(prhs[1]) ) {
        n = mxGetScalar(prhs[0]);
        k = mxGetScalar(prhs[1]);
        // Add checks here for valid n and k values (integer, k<=n, etc)
        result = 1;
        m = k < n-k ? k : n-k;
        for( i=1; i<=m; i++ ) {
            result = result * (n--) / i;
        }
        // Might need to add rounding here
        plhs[0] = mxCreateDoubleScalar(result);
    } else {
        mexErrMsgTxt("nchoosekx: Invalid inputs");
    }
}

Just spot checking results against a vpa implementation, this loop implementation does appear to give slightly more accurate results than the MATLAB nchoosek function when the inputs are too large to generate an exact double result. The mex loop implementation also runs faster. Maybe you can morph this code into what you are really wanting to do.

Answer 2

[Geoff Hayes]

Ameer - you are getting an error with the output parameter code because of your recursive call to test

 *y=(n*test(n-1,k-1))/k;

So two parameters are being passed in to test but three are required

 double test(double *y, double n, double k)

so you will need to either pass in the third parameter to test or re-work this method in order to return a value instead of updating the input y

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

Note also that you will want to re-think how you are calculating n choose k. Remember, that this is really

 n!/((n-k)!*(k!))

rather than the

 n!/k!

that you seem to have.

[James Tursa]

test(n,k) =

n! / ( (n-k)! * (k!) ) =

n * (n-1)! / ( (n-k)! * k * (k-1)! ) =

n * (n-1)! / ( ( (n-1) - (k-1) )! * (k-1)! ) / k =

n * test(n-1,k-1) / k

So the recursive call in the above code appears to be correct.

[Geoff Hayes]

hmmm...I was getting a different answer when comparing it with the MATLAB nchoosek. But now I'm not - weird!

Answer 3

[Ameer Ahmed]

Thanks James & Geoff.

I have tried doing it with two parameters as

double test(double n, double k) {

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

But it does not show the result. Just returns 0. I have implemented it in C without Mex Function. This works fine and returns the required answer. I don't know where is the mistake.

Actually, i want to implement a complex expression in C using Mex function which could replace three nchoosek(n,k) functions

nchoosek(n,k)*nchoosek(n,k)/nchoosek(n,k)

where each nchoosek() has different n and k.

I will be really helping if some one help me in implementing it.