Repository: sympy/sympy

Issue:
Benini distribution allows non valid parameters
Currently, the parameters used in Benini distribution can take non-positive values as well.
According to https://en.wikipedia.org/wiki/Benini_distribution the shape parameters (alpha and beta) and scale parameter(sigma) must be positive.
To reproduce the issue
```
>>> from sympy import Symbol
>>> from sympy.stats import Benini, density
>>> alpha = Symbol('alpha', positive = False)
>>> beta = Symbol('beta', positive = False)
>>> sigma = Symbol('sigma', positive = False)
>>> z = Symbol(z)
>>> X = Benini('X', alpha, beta, sigma) #this works, but it should not
```
Any permutation of signs of `alpha`, `beta` and `sigma` (with at least one of them being non positive) works fine

Discussion:

Patch:
diff --git a/sympy/ntheory/__init__.py b/sympy/ntheory/__init__.py
index bf00652217..31977979af 100644
--- a/sympy/ntheory/__init__.py
+++ b/sympy/ntheory/__init__.py
@@ -7,8 +7,7 @@
 from .primetest import isprime
 from .factor_ import divisors, factorint, multiplicity, perfect_power, \
     pollard_pm1, pollard_rho, primefactors, totient, trailing, divisor_count, \
-    divisor_sigma, factorrat, reduced_totient, primenu, primeomega, \
-         mersenne_prime_exponent, is_perfect, is_mersenne_prime
+    divisor_sigma, factorrat, reduced_totient, primenu, primeomega
 from .partitions_ import npartitions
 from .residue_ntheory import is_primitive_root, is_quad_residue, \
     legendre_symbol, jacobi_symbol, n_order, sqrt_mod, quadratic_residues, \
diff --git a/sympy/ntheory/factor_.py b/sympy/ntheory/factor_.py
index ac3c629836..2d3362822a 100644
--- a/sympy/ntheory/factor_.py
+++ b/sympy/ntheory/factor_.py
@@ -19,14 +19,6 @@
 from .primetest import isprime
 from .generate import sieve, primerange, nextprime
 
-
-# Note: This list should be updated whenever new Mersenne primes are found.
-# Refer: https://www.mersenne.org/
-MERSENNE_PRIME_EXPONENTS = (2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203,
- 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049,
- 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583,
- 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917, 82589933)
-
 small_trailing = [0] * 256
 for j in range(1,8):
     small_trailing[1<<j::1<<(j+1)] = [j] * (1<<(7-j))
@@ -2048,84 +2040,3 @@ def eval(cls, n):
                 raise ValueError("n must be a positive integer")
             else:
                 return sum(factorint(n).values())
-
-
-def mersenne_prime_exponent(nth):
-    """Returns the exponent ``i`` for the nth Mersenne prime (which
-    has the form `2^i - 1`).
-
-    Examples
-    ========
-
-    >>> from sympy.ntheory.factor_ import mersenne_prime_exponent
-    >>> mersenne_prime_exponent(1)
-    2
-    >>> mersenne_prime_exponent(20)
-    4423
-    """
-    n = as_int(nth)
-    if n < 1:
-        raise ValueError("nth must be a positive integer; mersenne_prime_exponent(1) == 2")
-    if n > 51:
-        raise ValueError("There are only 51 perfect numbers; nth must be less than or equal to 51")
-    return MERSENNE_PRIME_EXPONENTS[n - 1]
-
-
-def is_perfect(n):
-    """Returns True if ``n`` is a perfect number, else False.
-
-    A perfect number is equal to the sum of its positive, proper divisors.
-
-    Examples
-    ========
-
-    >>> from sympy.ntheory.factor_ import is_perfect, divisors
-    >>> is_perfect(20)
-    False
-    >>> is_perfect(6)
-    True
-    >>> sum(divisors(6)[:-1])
-    6
-
-    References
-    ==========
-
-    .. [1] http://mathworld.wolfram.com/PerfectNumber.html
-
-    """
-    from sympy.core.power import integer_log
-
-    r, b = integer_nthroot(1 + 8*n, 2)
-    if not b:
-        return False
-    n, x = divmod(1 + r, 4)
-    if x:
-        return False
-    e, b = integer_log(n, 2)
-    return b and (e + 1) in MERSENNE_PRIME_EXPONENTS
-
-
-def is_mersenne_prime(n):
-    """Returns True if  ``n`` is a Mersenne prime, else False.
-
-    A Mersenne prime is a prime number having the form `2^i - 1`.
-
-    Examples
-    ========
-
-    >>> from sympy.ntheory.factor_ import is_mersenne_prime
-    >>> is_mersenne_prime(6)
-    False
-    >>> is_mersenne_prime(127)
-    True
-
-    References
-    ==========
-
-    .. [1] http://mathworld.wolfram.com/MersennePrime.html
-
-    """
-    from sympy.core.power import integer_log
-
-    r, b = integer_log(n + 1, 2)
-    return b and r in MERSENNE_PRIME_EXPONENTS
diff --git a/sympy/printing/mathml.py b/sympy/printing/mathml.py
index d01846e92e..edf666395e 100644
--- a/sympy/printing/mathml.py
+++ b/sympy/printing/mathml.py
@@ -979,66 +979,6 @@ def _print_Basic(self, e):
         mrow.appendChild(brac)
         return mrow
 
-    def _print_Tuple(self, e):
-        mrow = self.dom.createElement('mrow')
-        x = self.dom.createElement('mfenced')
-        for arg in e.args:
-            x.appendChild(self._print(arg))
-        mrow.appendChild(x)
-        return mrow
-
-    def _print_Interval(self, i):
-        mrow = self.dom.createElement('mrow')
-        brac = self.dom.createElement('mfenced')
-        if i.start == i.end:
-            # Most often, this type of Interval is converted to a FiniteSet
-            brac.setAttribute('open', '{')
-            brac.setAttribute('close', '}')
-            brac.appendChild(self._print(i.start))
-        else:
-            if i.left_open:
-                brac.setAttribute('open', '(')
-            else:
-                brac.setAttribute('open', '[')
-
-            if i.right_open:
-                brac.setAttribute('close', ')')
-            else:
-                brac.setAttribute('close', ']')
-            brac.appendChild( self._print(i.start))
-            brac.appendChild( self._print(i.end))
-
-        mrow.appendChild(brac)
-        return mrow
-
-    def _print_Abs(self, expr, exp=None):
-        mrow = self.dom.createElement('mrow')
-        x = self.dom.createElement('mfenced')
-        x.setAttribute('open', '|')
-        x.setAttribute('close', '|')
-        x.appendChild(self._print(expr.args[0]))
-        mrow.appendChild(x)
-        return mrow
-
-    _print_Determinant = _print_Abs
-
-    def _print_re_im(self, c, expr):
-        mrow = self.dom.createElement('mrow')
-        mi = self.dom.createElement('mi')
-        mi.setAttribute('mathvariant', 'fraktur')
-        mi.appendChild(self.dom.createTextNode(c))
-        mrow.appendChild(mi)
-        brac = self.dom.createElement('mfenced')
-        brac.appendChild(self._print(expr))
-        mrow.appendChild(brac)
-        return mrow
-
-    def _print_re(self, expr, exp=None):
-        return self._print_re_im('R', expr.args[0])
-
-    def _print_im(self, expr, exp=None):
-        return self._print_re_im('I', expr.args[0])
-
     def _print_AssocOp(self, e):
         mrow = self.dom.createElement('mrow')
         mi = self.dom.createElement('mi')
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index 2de09f6614..122a1632f8 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -232,6 +232,12 @@ def Arcsin(name, a=0, b=1):
 class BeniniDistribution(SingleContinuousDistribution):
     _argnames = ('alpha', 'beta', 'sigma')
 
+    @staticmethod
+    def check(alpha, beta, sigma):
+        _value_check(alpha > 0, "Shape parameter Alpha must be positive.")
+        _value_check(beta > 0, "Shape parameter Beta must be positive.")
+        _value_check(sigma > 0, "Scale parameter Sigma must be positive.")
+
     @property
     def set(self):
         return Interval(self.sigma, oo)


Test patch:
diff --git a/sympy/ntheory/tests/test_factor_.py b/sympy/ntheory/tests/test_factor_.py
index 6226323e26..2bcb8b5e62 100644
--- a/sympy/ntheory/tests/test_factor_.py
+++ b/sympy/ntheory/tests/test_factor_.py
@@ -13,7 +13,7 @@
     factorrat, reduced_totient)
 from sympy.ntheory.factor_ import (smoothness, smoothness_p,
     antidivisors, antidivisor_count, core, digits, udivisors, udivisor_sigma,
-    udivisor_count, primenu, primeomega, small_trailing, mersenne_prime_exponent, is_perfect, is_mersenne_prime)
+    udivisor_count, primenu, primeomega, small_trailing)
 from sympy.ntheory.generate import cycle_length
 from sympy.ntheory.multinomial import (
     multinomial_coefficients, multinomial_coefficients_iterator)
@@ -562,30 +562,3 @@ def test_primeomega():
     assert primeomega(n)
     assert primeomega(n).subs(n, 2 ** 31 - 1) == 1
     assert summation(primeomega(n), (n, 2, 30)) == 59
-
-
-def test_mersenne_prime_exponent():
-    assert mersenne_prime_exponent(1) == 2
-    assert mersenne_prime_exponent(4) == 7
-    assert mersenne_prime_exponent(10) == 89
-    assert mersenne_prime_exponent(25) == 21701
-    raises(ValueError, lambda: mersenne_prime_exponent(52))
-    raises(ValueError, lambda: mersenne_prime_exponent(0))
-
-
-def test_is_perfect():
-    assert is_perfect(6) is True
-    assert is_perfect(15) is False
-    assert is_perfect(28) is True
-    assert is_perfect(400) is False
-    assert is_perfect(496) is True
-    assert is_perfect(8128) is True
-    assert is_perfect(10000) is False
-
-
-def test_is_mersenne_prime():
-    assert is_mersenne_prime(10) is False
-    assert is_mersenne_prime(127) is True
-    assert is_mersenne_prime(511) is False
-    assert is_mersenne_prime(131071) is True
-    assert is_mersenne_prime(2147483647) is True
diff --git a/sympy/printing/tests/test_mathml.py b/sympy/printing/tests/test_mathml.py
index c6fc19f4f1..1c6b33572c 100644
--- a/sympy/printing/tests/test_mathml.py
+++ b/sympy/printing/tests/test_mathml.py
@@ -2,19 +2,14 @@
     tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, E, I, oo, \
     pi, GoldenRatio, EulerGamma, Sum, Eq, Ne, Ge, Lt, Float, Matrix, Basic, S, \
     MatrixSymbol, Function, Derivative, log
-from sympy.core.containers import Tuple
-from sympy.functions.elementary.complexes import re, im, Abs, conjugate
 from sympy.functions.combinatorial.factorials import factorial, factorial2, binomial
 from sympy.functions.elementary.complexes import conjugate
 from sympy.functions.special.zeta_functions import polylog, lerchphi
 from sympy.logic.boolalg import And, Or, Implies, Equivalent, Xor, Not
-from sympy.matrices.expressions.determinant import Determinant
-from sympy.printing.mathml import mathml, MathMLContentPrinter, MathMLPresentationPrinter, \
-    MathMLPrinter
 from sympy.sets.sets import FiniteSet, Union, Intersection, Complement, SymmetricDifference
 from sympy.stats.rv import RandomSymbol
-from sympy.sets.sets import Interval
-from sympy.stats.rv import RandomSymbol
+from sympy.printing.mathml import mathml, MathMLContentPrinter, MathMLPresentationPrinter, \
+    MathMLPrinter
 from sympy.utilities.pytest import raises
 
 x = Symbol('x')
@@ -915,42 +910,6 @@ def test_presentation_mathml_order():
     assert mml.childNodes[6].childNodes[1].childNodes[0].nodeValue == '3'
 
 
-def test_print_intervals():
-    a = Symbol('a', real=True)
-    assert mpp.doprint(Interval(0, a)) == '<mrow><mfenced close="]" open="["><mn>0</mn><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Interval(0, a, False, False)) == '<mrow><mfenced close="]" open="["><mn>0</mn><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Interval(0, a, True, False)) == '<mrow><mfenced close="]" open="("><mn>0</mn><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Interval(0, a, False, True)) == '<mrow><mfenced close=")" open="["><mn>0</mn><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Interval(0, a, True, True)) == '<mrow><mfenced close=")" open="("><mn>0</mn><mi>a</mi></mfenced></mrow>'
-
-
-def test_print_tuples():
-    a = Symbol('a')
-    assert mpp.doprint(Tuple(0,)) == '<mrow><mfenced><mn>0</mn></mfenced></mrow>'
-    assert mpp.doprint(Tuple(0, a)) == '<mrow><mfenced><mn>0</mn><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Tuple(0, a, a)) == '<mrow><mfenced><mn>0</mn><mi>a</mi><mi>a</mi></mfenced></mrow>'
-    assert mpp.doprint(Tuple(0, 1, 2, 3, 4)) == '<mrow><mfenced><mn>0</mn><mn>1</mn><mn>2</mn><mn>3</mn><mn>4</mn></mfenced></mrow>'
-    assert mpp.doprint(Tuple(0, 1, Tuple(2, 3, 4))) == '<mrow><mfenced><mn>0</mn><mn>1</mn><mrow><mfenced><mn>2</mn><mn>3</mn><mn>4</mn></mfenced></mrow></mfenced></mrow>'
-
-
-def test_print_re_im():
-    x = Symbol('x')
-    assert mpp.doprint(re(x)) == '<mrow><mi mathvariant="fraktur">R</mi><mfenced><mi>x</mi></mfenced></mrow>'
-    assert mpp.doprint(im(x)) == '<mrow><mi mathvariant="fraktur">I</mi><mfenced><mi>x</mi></mfenced></mrow>'
-    assert mpp.doprint(re(x + 1)) == '<mrow><mrow><mi mathvariant="fraktur">R</mi><mfenced><mi>x</mi></mfenced></mrow><mo>+</mo><mn>1</mn></mrow>'
-    assert mpp.doprint(im(x + 1)) == '<mrow><mi mathvariant="fraktur">I</mi><mfenced><mi>x</mi></mfenced></mrow>'
-
-
-def test_print_Abs():
-    x = Symbol('x')
-    assert mpp.doprint(Abs(x)) == '<mrow><mfenced close="|" open="|"><mi>x</mi></mfenced></mrow>'
-    assert mpp.doprint(Abs(x + 1)) == '<mrow><mfenced close="|" open="|"><mrow><mi>x</mi><mo>+</mo><mn>1</mn></mrow></mfenced></mrow>'
-
-
-def test_print_Determinant():
-    assert mpp.doprint(Determinant(Matrix([[1, 2], [3, 4]]))) == '<mrow><mfenced close="|" open="|"><mfenced close="]" open="["><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>3</mn></mtd><mtd><mn>4</mn></mtd></mtr></mtable></mfenced></mfenced></mrow>'
-
-
 def test_presentation_settings():
     raises(TypeError, lambda: mathml(Symbol("x"), printer='presentation',method="garbage"))
 
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 5e86511ae7..1d0523c1cd 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -188,12 +188,25 @@ def test_arcsin():
 
 def test_benini():
     alpha = Symbol("alpha", positive=True)
-    b = Symbol("beta", positive=True)
+    beta = Symbol("beta", positive=True)
     sigma = Symbol("sigma", positive=True)
 
-    X = Benini('x', alpha, b, sigma)
-    assert density(X)(x) == ((alpha/x + 2*b*log(x/sigma)/x)
-                          *exp(-alpha*log(x/sigma) - b*log(x/sigma)**2))
+    X = Benini('x', alpha, beta, sigma)
+    assert density(X)(x) == ((alpha/x + 2*beta*log(x/sigma)/x)
+                          *exp(-alpha*log(x/sigma) - beta*log(x/sigma)**2))
+
+    alpha = Symbol("alpha", positive=False)
+    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
+
+    beta = Symbol("beta", positive=False)
+    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
+
+    alpha = Symbol("alpha", positive=True)
+    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
+
+    beta = Symbol("beta", positive=True)
+    sigma = Symbol("sigma", positive=False)
+    raises(ValueError, lambda: Benini('x', alpha, beta, sigma))
 
 
 def test_beta():