irk_geometric_vec now handles p==1.0 without MKL error

Closes issue #15.

Note that MKL samples from Geometric(p) are supported on
the set of non-negative 32-bit integers, while samples from
numpy.random are spported on 64-bit positive integers.

```
In [1]: import mkl_random as rnd

In [2]: rnd.geometric(1, 10)
Out[2]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int32)

In [3]: import numpy as np

In [4]: rnd2 = np.random.default_rng()

In [5]: rnd2
Out[5]: Generator(PCG64) at 0x2AF349B17A50

In [6]: rnd2.geometric(1, 10)
Out[6]: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

In [7]: rnd2.geometric(1, 10)-1
Out[7]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [8]: np.bincount(rnd2.geometric(0.78, 10**5)-1)
Out[8]:
array([77870, 17233,  3786,   879,   177,    44,     7,     3,     0,
           0,     1])

In [9]: np.bincount(rnd.geometric(0.78, 10**5))
Out[9]: array([78121, 17088,  3743,   821,   172,    42,    10,     2,     1])

In [10]: np.bincount(rnd.geometric(0.78, 10**5))
Out[10]: array([78029, 17225,  3658,   848,   179,    48,    10,     3])
```

Author: oleksandr-pavlyk

mkl_random/mklrand.pyx

@@ -4894,8 +4894,8 @@ cdef class RandomState:
 
         fp = PyFloat_AsDouble(p)
         if not PyErr_Occurred():
-            if fp < 0.0:
-                raise ValueError("p < 0.0")
+            if fp <= 0.0:
+                raise ValueError("p <= 0.0")
             if fp > 1.0:
                 raise ValueError("p > 1.0")
             return vec_discd_array_sc(self.internal_state, irk_geometric_vec, size, fp,
@@ -4905,7 +4905,7 @@ cdef class RandomState:
 
 
         op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY)
-        if np.any(np.less(op, 0.0)):
+        if np.any(np.less_equal(op, 0.0)):
             raise ValueError("p < 0.0")
         if np.any(np.greater(op, 1.0)):
             raise ValueError("p > 1.0")

mkl_random/src/mkl_distributions.cpp

@@ -1169,17 +1169,26 @@ irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p)
     if(len < 1)
         return;
 
-    while(len > MKL_INT_MAX) {
-        err = viRngGeometric(VSL_RNG_METHOD_GEOMETRIC_ICDF, state->stream, MKL_INT_MAX, res, p);
-        assert(err == VSL_STATUS_OK);
-
-        res += MKL_INT_MAX;
-        len -= MKL_INT_MAX;
-    }
+    if ((0.0 < p) && (p < 1.0)) {
+        while(len > MKL_INT_MAX) {
+            err = viRngGeometric(VSL_RNG_METHOD_GEOMETRIC_ICDF, state->stream, MKL_INT_MAX, res, p);
+            assert(err == VSL_STATUS_OK);
 
-    err = viRngGeometric(VSL_RNG_METHOD_GEOMETRIC_ICDF, state->stream, len, res, p);
-    assert(err == VSL_STATUS_OK);
+            res += MKL_INT_MAX;
+            len -= MKL_INT_MAX;
+        }
 
+        err = viRngGeometric(VSL_RNG_METHOD_GEOMETRIC_ICDF, state->stream, len, res, p);
+        assert(err == VSL_STATUS_OK);
+    } else {
+	if (p==1.0) {
+	    npy_intp i;
+	    for(i=0; i < len; ++i) res[i] = 0;
+	} else {
+	    assert(p >= 0.0);
+	    assert(p <= 1.0);
+	}
+    }
 }
 
 void