这个有前辈已经实现了，部分代码如下——

当然了，提取时不要忘记提取底层函数和一些个全局变量定义。
p.s.:怎么找相关函数？

在R下用UltraEdit搜一下，轻松搞定！

//////////////////////////RNG.h/////////////////////////////////////

#defined(RNG_EXPORTS)
#define RNG_API __declspec(dllexport)
#else
#define RNG_API __declspec(dllimport)
#endif

#ifdef __cplusplus
extern "C" {
#endif

RNG_API void __stdcall rmultinom(long* n, double* prob, long* size, long* K, long* randomNum);
#ifdef __cplusplus
}
#endif

//////////////////////////RNG.cpp/////////////////////////////////////
#include "RNG.h"
…

void __stdcall rmultinom(long* n, double* prob, long* size, long* K, long* randomNum)
{
  
  long* rN = randomNum;
  for(long i = 0; i < n[0]; i++)
  {
    __rmultinom(size[0], prob, K[0], rN);
    for(long j = 0; j < K[0]; j++){
        randomNum+j] = rN[j];
    }
  }
  return;

}

///以下来自R package：

void __rmultinom(long n, double* prob, long K, long* rN)
{
  long k;
  double pp, p_tot = 0.;

  for(k = 0; k < K; k++) {
    pp = prob[k];
    p_tot += pp;
    rN[k] = 0;
  }
  if (n == 0) return;
  if (K == 1 && p_tot == 0.) return;

  for(k = 0; k < K-1; k++) {
    if(prob[k]) {
        pp = prob[k] / p_tot;
        rN[k] = ((pp < 1.) ? (long) __rbinom((double) n, pp) :   n);
        n -= rN[k];
    }
    else rN[k] = 0;
    if(n <= 0) return;
    p_tot -= prob[k];
  }
  rN[K-1] = n;
  return;
}

///以下代码来自R
long __rbinom(double nin, double pp)
{
  static double c, fm, npq, p1, p2, p3, p4, qn;
  static double xl, xll, xlr, xm, xr;

  static double psave = -1.0;
  static long nsave = -1;
  static long m;

  double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
  double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
  long i,ix,k,n;

  n = (long)floor(nin + 0.5);
  
  if (n == 0 || pp == 0.) return 0;
  if (pp == 1.) return n;

  p = fmin2(pp, 1. – pp);
  q = 1. – p;
  np = n * p;
  r = p / q;
  g = r * (n + 1);

  /* Setup, perform only when parameters change [using static (globals): */

  /* FIXING: Want this thread safe
    — use as little (thread globals) as possible
  */
  if (pp != psave || n != nsave) {
  psave = pp;
  nsave = n;
  if (np < 30.0) {
    /* inverse cdf logic for mean less than 30 */
    qn = pow(q, (double) n);
    goto L_np_small;
  } else {
    ffm = np + p;
    m = (long)ffm;
    fm = m;
    npq = np * q;
    p1 = (long)(2.195 * sqrt(npq) – 4.6 * q) + 0.5;
    xm = fm + 0.5;
    xl = xm – p1;
    xr = xm + p1;
    c = 0.134 + 20.5 / (15.3 + fm);
    al = (ffm – xl) / (ffm – xl * p);
    xll = al * (1.0 + 0.5 * al);
    al = (xr – ffm) / (xr * q);
    xlr = al * (1.0 + 0.5 * al);
    p2 = p1 * (1.0 + c + c);
    p3 = p2 + c / xll;
    p4 = p3 + c / xlr;
  }
  } else if (n == nsave) {
  if (np < 30.0)
    goto L_np_small;
  }

  /*————————– np = n*p >= 30 : ——————- */
  repeat {
    u = unif_rand() * p4;
    v = unif_rand();
    /* triangular region */
    if (u <= p1) {
    ix = (long)(xm – p1 * v + u);
    goto finis;
    }
    /* parallelogram region */
    if (u <= p2) {
    x = xl + (u – p1) / c;
    v = v * c + 1.0 – fabs(xm – x) / p1;
    if (v > 1.0 || v <= 0.)
      continue;
    ix = (long)x;
    } else {
    if (u > p3) {   /* right tail */
      ix = (long)(xr – log(v) / xlr);
      if (ix > n)
      continue;
      v = v * (u – p3) * xlr;
    } else {/* left tail */
      ix = (long)(xl + log(v) / xll);
      if (ix < 0)
      continue;
      v = v * (u – p2) * xll;
    }
    }
    /* determine appropriate way to perform accept/reject test */
    k = abs(ix – m);
    if (k <= 20 || k >= npq / 2 – 1) {
    /* explicit evaluation */
    f = 1.0;
    if (m < ix) {
      for (i = m + 1; i <= ix; i++)
      f *= (g / i – r);
    } else if (m != ix) {
      for (i = ix + 1; i <= m; i++)
      f /= (g / i – r);
    }
    if (v <= f)
      goto finis;
    } else {
    /* squeezing using upper and lower bounds on log(f(x)) */
    amaxp = (k / npq) * ((k * (k / 3. + 0.625) + 0.1666666666666) / npq + 0.5);
    ynorm = -k * k / (2.0 * npq);
    alv = log(v);
    if (alv < ynorm – amaxp)
      goto finis;
    if (alv <= ynorm + amaxp) {
      /* stirling's formula to machine accuracy */
      /* for the final acceptance/rejection test */
      x1 = ix + 1;
      f1 = fm + 1.0;
      z = n + 1 – fm;
      w = n – ix + 1.0;
      z2 = z * z;
      x2 = x1 * x1;
      f2 = f1 * f1;
      w2 = w * w;
      if (alv <= xm * log(f1 / x1) + (n – m + 0.5) * log(z / w) + (ix – m) * log(w * p / (x1 * q)) + (13860.0 – (462.0 – (132.0 – (99.0 – 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 – (462.0 – (132.0 – (99.0 – 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 – (462.0 – (132.0 – (99.0 – 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 – (462.0 – (132.0 – (99.0 – 140.0 / w2) / w2) / w2) / w2) / w / 166320.)
      goto finis;
    }
    }
}

L_np_small:
  /*———————- np = n*p < 30 : ————————- */

repeat {
  ix = 0;
  f = qn;
  u = unif_rand();
  repeat {
  if (u < f)
      goto finis;
  if (ix > 110)
      break;
  u -= f;
  ix++;
  f *= (g / ix – r);
  }
}

finis:
  if (psave > 0.5)
  ix = n – ix;
return ix;
}