Computing the distribution of a Jonckheere-Terpstra statistic I have the following C++ function which I'd like to speed up, if possible. This program will have many different users, so parallelization is not really an option.

The C++ function is wrapped by C code and then called from R, using:

R CMD SHLIB -o jonckterp.so wrapper.cpp inside.cpp


from the command line, then:

dyn.load("jonckterp.so")


from within R.

inside.cpp

#include &lt;stdio.h&gt;
#include &lt;vector&gt;
#include &lt;iostream&gt;
#include &lt;R.h&gt;
#include &lt;Rmath.h&gt;

extern "C" void jonckterp(int *x,int *c, int *ngrp, double *comb);

int *finalcountdown;
static int count;

static int *tabInit(int num){
  int *w;
  w = (int *) R_alloc(num, sizeof(int));
  memset(w, '\0', num * sizeof(int));
  return(w);
}

static int mannwhit(std::vector&lt;int&gt; c, std::vector&lt;int&gt; sizes){
  int Ustat=0;
  int start=0;
  size_t Ngroups=sizes.size();
   for(int i=0; i&lt;Ngroups-1;i++){
    for(int j=start;j&lt;(start+sizes[i]);j++){
      for(int k=start+sizes[i];k&lt;c.size();k++){
    if(c[j]&lt;c[k])
      Ustat+=1;
      }
    } 
    start+=sizes[i];
  }
   return(Ustat);
    }

void allocate_to_groups_impl( std::vector&lt;char&gt;&amp; usage, std::vector&lt;int&gt;&amp; order, const           std::vector&lt;int&gt;&amp; cumsum, int group, int place, int min )
{
std::vector&lt;int&gt; store(order.size());
std::vector&lt;int&gt; sizes(cumsum.size());
sizes[0]=cumsum[0];
for(int i=1;i&lt;cumsum.size();i++)
sizes[i]=(cumsum[i]-cumsum[i-1]);
    if (place == cumsum[group]) {
    group++; min = 0;
    if (group == cumsum.size()) {
      int ij=0;
      for( std::vector&lt;int&gt;::iterator it = order.begin(); it != order.end(); ++it ){
store[ij]=*it;
    ij++;
  }
  finalcountdown[count]=mannwhit(store,sizes);
  count++;
  return;
    }
  }

  for( int v = min, max = usage.size() + place - cumsum[group]; v &lt;= max; ++v ) {
    if (!usage[v]) {
    order[place] = v;
    usage[v] = 1;
    allocate_to_groups_impl(usage, order, cumsum, group, place+1, v+1);
  usage[v] = 0;
    }
  }
}


void allocate_to_groups( int *c, int Ngroups )
{
  size_t sum_of_c = 0;
  std::vector&lt;int&gt; cumsum_of_c;
  for( int* it = c; it &lt; c + Ngroups; ++it ){
    cumsum_of_c.push_back(sum_of_c += *it);
  }
  std::vector&lt;int&gt; order(sum_of_c);
  std::vector&lt;char&gt; usage(sum_of_c);
  allocate_to_groups_impl(usage, order, cumsum_of_c, 0, 0, 0);
}

void jonckterp(int *x,int *c, int *ngrp, double *numcomb){
  count=0;  
  float p=0;
  int d[*ngrp];         
  for(int i=0; i&lt; *ngrp; i++)
    d[i]=c[i];
  finalcountdown =tabInit(*numcomb);
  allocate_to_groups(d,*ngrp);
  for(int test=0;test&lt;*numcomb;test++){
  if(finalcountdown[test]&gt;=*x){
    p+=1;   
}
}
*x=p;
}


wrapper.cpp

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt; 
#include &lt;string.h&gt;
extern "C" {
void jonckterp(int *x,int *c, int *ngrp, double *numcomb);
}


jonckterp.R

early.fn&lt;-function(stat,group.sizes){
numcomb&lt;-factorial(sum(group.sizes))/prod(factorial(group.sizes))
.C("jonckterp",p=as.integer(stat),as.integer(group.sizes),as.integer(length(group.sizes)),as.double(numcomb))$p/numcomb  
}
I didn't read in it's entirity, but change this line

static int mannwhit(std::vector&lt;int&gt; c, std::vector&lt;int&gt; sizes)


to

static int mannwhit(const std::vector&lt;int&gt;&amp; c, 
                    const std::vector&lt;int&gt;&amp; sizes)


This won't speed up your code but you should use

#include &lt;cstdio&gt; //instead of stdio.h


and memset will need to be prefixed as std::memset.

Also ensure you compiler is compiling at 02 and with out debug flags.