Rewritten mkl_random.multinomial to use MKL's own generator

```
In [1]: import mkl_random, numpy as np

In [2]: %time mkl_random.multinomial(527, np.array([1/16, 3/16, 5/16, 7/16], 'd'), size=10**6)
CPU times: user 820 ms, sys: 0 ns, total: 820 ms
Wall time: 819 ms
Out[2]:
array([[ 29,  90, 176, 232],
       [ 32,  91, 185, 219],
       [ 35,  90, 169, 233],
       ...,
       [ 42,  93, 157, 235],
       [ 28, 105, 152, 242],
       [ 32,  79, 179, 237]], dtype=int32)

In [3]: %timeit mkl_random.multinomial(527, np.array([1/16, 3/16, 5/16, 7/16], 'd'), size=10**6)
803 ms ± 9.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [4]: %timeit mkl_random.multinomial(527, np.array([1/16, 3/16, 5/16, 7/16], 'd'), size=10**5)
80.5 ms ± 874 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

Previously,

```
In [1]: import mkl_random, numpy as np

In [2]: %timeit mkl_random.multinomial(527, np.array([1/16, 3/16, 5/16, 7/16], 'd'), size=10**6)
1.35 s ± 2.35 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit mkl_random.multinomial(527, np.array([1/16, 3/16, 5/16, 7/16], 'd'), size=10**5)
135 ms ± 613 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

Author: oleksandr-pavlyk

mkl_random/mklrand.pyx

@@ -131,6 +131,7 @@ cdef extern from "mkl_distributions.h":
     void irk_geometric_vec(irk_state *state, npy_intp len, int *res, double p) nogil
     void irk_negbinomial_vec(irk_state *state, npy_intp len, int *res, double a, double p) nogil
     void irk_binomial_vec(irk_state *state, npy_intp len, int *res, int n, double p) nogil
+    void irk_multinomial_vec(irk_state *state, npy_intp len, int *res, int n, int d, double *pvec) nogil
     void irk_hypergeometric_vec(irk_state *state, npy_intp len, int *res, int ls, int ss, int ms) nogil
 
     void irk_poisson_vec_PTPE(irk_state *state, npy_intp len, int *res, double lam) nogil
@@ -5544,27 +5545,13 @@ cdef class RandomState:
             raise ValueError("sum(pvals[:-1]) > 1.0")
 
         shape = _shape_from_size(size, d)
-
         multin = np.zeros(shape, np.int32)
+
         mnarr = <ndarray>multin
         mnix = <int*>PyArray_DATA(mnarr)
         sz = PyArray_SIZE(mnarr)
-        with self.lock, nogil, cython.cdivision(True):
-            i = 0
-            while i < sz:
-                Sum = 1.0
-                dn = n
-                for j from 0 <= j < d-1:
-#                    mnix[i+j] = irk_binomial(self.internal_state, dn, pix[j]/Sum)
-                    irk_binomial_vec(self.internal_state, 1, mnix + (i+j), dn, pix[j]/Sum)
-                    dn = dn - mnix[i+j]
-                    if dn <= 0:
-                        break
-                    Sum = Sum - pix[j]
-                if dn > 0:
-                    mnix[i+d-1] = dn
-
-                i = i + d
+
+        irk_multinomial_vec(self.internal_state, sz // d, mnix, n, d, pix)
 
         return multin
 

mkl_random/src/mkl_distributions.cpp

@@ -1136,6 +1136,31 @@ irk_binomial_vec(irk_state *state, npy_intp len, int *res, const int n, const do
     }
 }
 
+void
+irk_multinomial_vec(irk_state *state, npy_intp len, int *res, const int n, const int k, const double *pvec)
+{
+    int err;
+
+    if(len < 1)
+        return;
+
+    if (n==0) {
+        memset(res, 0, len*k*sizeof(int));
+    }
+    else {
+        while (len > MKL_INT_MAX) {
+            err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, state->stream, MKL_INT_MAX, res, n, k, pvec);
+            assert(err == VSL_STATUS_OK);
+            res += k*MKL_INT_MAX;
+            len -= k*MKL_INT_MAX;
+        }
+
+        err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, state->stream, len, res, n, k, pvec);
+        assert(err == VSL_STATUS_OK);
+    }
+}
+
+
 void
 irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p)
 {
@@ -1953,184 +1978,3 @@ static double irk_double(irk_state *state) {
 
     return res;
 }
-
-#ifdef COMPILE_UNUSED_CODE
-
-#ifndef min
-#define min(x,y) (((x)<(y))?(x):(y))
-#endif
-
-static long irk_binomial_btpe(irk_state *state, long n, double p)
-{
-    double r,q,fm,p1,xm,xl,xr,c,laml,lamr,p2,p3,p4;
-    double a,u,v,s,F,rho,t,A,nrq,x1,x2,f1,f2,z,z2,w,w2,x;
-    long m,y,k,i;
-
-    r = min(p, 1.0-p);
-    q = 1.0 - r;
-    fm = n*r+r;
-    m = (long)floor(fm);
-    p1 = floor(2.195*sqrt(n*r*q)-4.6*q) + 0.5;
-    xm = m + 0.5;
-    xl = xm - p1;
-    xr = xm + p1;
-    c = 0.134 + 20.5/(15.3 + m);
-    a = (fm - xl)/(fm-xl*r);
-    laml = a*(1.0 + a/2.0);
-    a = (xr - fm)/(xr*q);
-    lamr = a*(1.0 + a/2.0);
-    p2 = p1*(1.0 + 2.0*c);
-    p3 = p2 + c/laml;
-    p4 = p3 + c/lamr;
-
-  /* sigh ... */
-  Step10:
-    nrq = n*r*q;
-    u = irk_double(state)*p4;
-    v = irk_double(state);
-    if (u > p1) goto Step20;
-    y = (long)floor(xm - p1*v + u);
-    goto Step60;
-
-  Step20:
-    if (u > p2) goto Step30;
-    x = xl + (u - p1)/c;
-    v = v*c + 1.0 - fabs(m - x + 0.5)/p1;
-    if (v > 1.0) goto Step10;
-    y = (long)floor(x);
-    goto Step50;
-
-  Step30:
-    if (u > p3) goto Step40;
-    y = (long)floor(xl + log(v)/laml);
-    if (y < 0) goto Step10;
-    v = v*(u-p2)*laml;
-    goto Step50;
-
-  Step40:
-    y = (long)floor(xr - log(v)/lamr);
-    if (y > n) goto Step10;
-    v = v*(u-p3)*lamr;
-
-  Step50:
-    k = labs(y - m);
-    if ((k > 20) && (k < ((nrq)/2.0 - 1))) goto Step52;
-
-    s = r/q;
-    a = s*(n+1);
-    F = 1.0;
-    if (m < y)
-    {
-        for (i=m+1; i<=y; i++)
-        {
-            F *= (a/i - s);
-        }
-    }
-    else if (m > y)
-    {
-        for (i=y+1; i<=m; i++)
-        {
-            F /= (a/i - s);
-        }
-    }
-    if (v > F) goto Step10;
-    goto Step60;
-
-    Step52:
-    rho = (k/(nrq))*((k*(k/3.0 + 0.625) + 0.16666666666666666)/nrq + 0.5);
-    t = -k*k/(2*nrq);
-    A = log(v);
-    if (A < (t - rho)) goto Step60;
-    if (A > (t + rho)) goto Step10;
-
-    x1 = y+1;
-    f1 = m+1;
-    z = n+1-m;
-    w = n-y+1;
-    x2 = x1*x1;
-    f2 = f1*f1;
-    z2 = z*z;
-    w2 = w*w;
-    if (A > (xm*log(f1/x1)
-           + (n-m+0.5)*log(z/w)
-           + (y-m)*log(w*r/(x1*q))
-           + (13680.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320.
-           + (13680.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320.
-           + (13680.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320.
-           + (13680.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.))
-    {
-        goto Step10;
-    }
-
-  Step60:
-    if (p > 0.5)
-    {
-        y = n - y;
-    }
-
-    return y;
-}
-
-static long irk_binomial_inversion(irk_state *state, long n, double p)
-{
-    double q, qn, np, px, U;
-    long X, bound;
-
-    q = 1.0 - p;
-    qn = exp(n * log(q));
-    np = n*p;
-    bound = min(n, np + 10.0*sqrt(np*q + 1));
-
-    X = 0;
-    px = qn;
-    U = irk_double(state);
-    while (U > px)
-    {
-        X++;
-        if (X > bound)
-        {
-            X = 0;
-            px = qn;
-            U = irk_double(state);
-        } else
-        {
-            U -= px;
-            px  = ((n-X+1) * p * px)/(X*q);
-        }
-    }
-    return X;
-}
-
-#undef min
-
-static long irk_binomial(irk_state *state, long n, double p)
-{
-    double q;
-
-    if (p <= 0.5)
-    {
-        if (p*n <= 30.0)
-        {
-            return irk_binomial_inversion(state, n, p);
-        }
-        else
-        {
-            return irk_binomial_btpe(state, n, p);
-        }
-    }
-    else
-    {
-        q = 1.0-p;
-        if (q*n <= 30.0)
-        {
-            return n - irk_binomial_inversion(state, n, q);
-        }
-        else
-        {
-            return n - irk_binomial_btpe(state, n, q);
-        }
-    }
-
-}
-
-#endif

mkl_random/src/mkl_distributions.h

@@ -88,6 +88,9 @@ extern void irk_noncentral_f_vec(irk_state *state, npy_intp len, double *res, do
 extern void irk_triangular_vec(irk_state *state, npy_intp len, double *res, double left, double mode, double right);
 
 extern void irk_binomial_vec(irk_state *state, npy_intp len, int *res, const int n, const double p);
+
+extern void irk_multinomial_vec(irk_state *state, npy_intp len, int *res, const int n, const int k, const double* pvec);
+
 extern void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p);
 extern void irk_negbinomial_vec(irk_state *state, npy_intp len, int *res, const double a, const double p);
 extern void irk_hypergeometric_vec(irk_state *state, npy_intp len, int *res, const int ls, const int ss, const int ms);