Repository: sympy/sympy

Issue:
Crash when calculating density of sum of normal and normal to the 5:th power
Calculating the density of the sum of the fifth power of a normal distributed random variable and the random variable raises. Lower powers seems to work, but not higher. I think, although I am not sure, that the calculations are correct. It is just that the calculation at some point leads to a `ConditionSet` perhaps via a `Piecewise`.

```python

import sympy.stats as stats

# Crashes
X = stats.Normal('X', 0, sympy.Symbol("sigma")**2)
stats.density(X**5 + X)

# Works
X = stats.Normal('X', 0, sympy.Symbol("sigma")**2)
stats.density(X**3 + X)
```

```
Traceback (most recent call last):

  File "/tmp/ipykernel_17323/2537738759.py", line 1, in <module>
    stats.density(X**5 + X)

  File "/usr/lib/python3/dist-packages/sympy/stats/rv.py", line 922, in density
    return Density(expr, condition).doit(evaluate=evaluate, **kwargs)

  File "/usr/lib/python3/dist-packages/sympy/stats/rv.py", line 869, in doit
    result = pspace(expr).compute_density(expr, **kwargs)

  File "/usr/lib/python3/dist-packages/sympy/stats/crv.py", line 529, in compute_density
    fy = sum(fx(g) * abs(g.diff(y)) for g in gs)

TypeError: 'ConditionSet' object is not iterable
```

Discussion:
sympy-bot: 

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

:x: There was an issue with the release notes. **Please do not close this pull request;** instead edit the description after reading the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

* No release notes were found. Please edit the PR description and write the release notes under `<!-- BEGIN RELEASE NOTES -->`.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->
    Fixes #24915 
    Continuing the apparently stale pull request #25373, but applying the last review comments.
    #### Brief description of what is fixed or changed
    For the density of the sum of normal and normal to the 5th power or higher power, the ValueError exception is raised with informative messages instead of calculating it as per the Abel-Ruffini theorem.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

    * physics.units
      * Corrected a semantical error in the conversion between volt and statvolt which
        reported the volt as being larger than the statvolt.

    or if no release note(s) should be included use:

    NO ENTRY

    See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
    information on how to write release notes. The bot will check your release
    notes automatically to see if they are formatted correctly. -->
    - stats 
      - The `compute_density` function adheres to the Abel-Ruffini theorem, raising a ValueError.


    <!-- BEGIN RELEASE NOTES -->

    <!-- END RELEASE NOTES -->


</details><p>

github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
| Change   | Before [a00718ba]    | After [36bcecfe]    |   Ratio | Benchmark (Parameter)                                                |
|----------|----------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 67.6±1ms             | 44.6±0.2ms          |    0.66 | integrate.TimeIntegrationRisch02.time_doit(10)                       |
| -        | 66.9±0.7ms           | 43.3±0.3ms          |    0.65 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 19.0±0.1μs           | 29.4±0.3μs          |    1.55 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 5.28±0.01ms          | 2.80±0.02ms         |    0.53 | logic.LogicSuite.time_load_file                                      |
| -        | 73.1±0.8ms           | 28.4±0.05ms         |    0.39 | polys.TimeGCD_GaussInt.time_op(1, 'dense')                           |
| -        | 25.6±0.1ms           | 16.8±0.08ms         |    0.66 | polys.TimeGCD_GaussInt.time_op(1, 'expr')                            |
| -        | 71.8±0.6ms           | 28.9±0.2ms          |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'sparse')                          |
| -        | 253±3ms              | 125±0.2ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'dense')                           |
| -        | 254±3ms              | 124±0.6ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'sparse')                          |
| -        | 646±3ms              | 371±2ms             |    0.57 | polys.TimeGCD_GaussInt.time_op(3, 'dense')                           |
| -        | 652±2ms              | 372±2ms             |    0.57 | polys.TimeGCD_GaussInt.time_op(3, 'sparse')                          |
| -        | 492±3μs              | 286±3μs             |    0.58 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense')            |
| -        | 1.78±0.01ms          | 1.05±0ms            |    0.59 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense')            |
| -        | 5.76±0.03ms          | 3.06±0.02ms         |    0.53 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense')            |
| -        | 442±4μs              | 229±1μs             |    0.52 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense')               |
| -        | 1.46±0.01ms          | 677±3μs             |    0.46 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense')               |
| -        | 4.92±0.02ms          | 1.65±0.01ms         |    0.34 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense')               |
| -        | 373±2μs              | 205±0.9μs           |    0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense')                |
| -        | 2.43±0.01ms          | 1.22±0ms            |    0.5  | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense')                |
| -        | 9.94±0.08ms          | 4.33±0.02ms         |    0.44 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense')                |
| -        | 353±1μs              | 168±0.8μs           |    0.48 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense')            |
| -        | 2.50±0.03ms          | 890±7μs             |    0.36 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense')            |
| -        | 9.42±0.08ms          | 2.62±0.02ms         |    0.28 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense')            |
| -        | 1.03±0.01ms          | 428±3μs             |    0.42 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense')           |
| -        | 1.73±0.01ms          | 503±0.8μs           |    0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 5.88±0.05ms          | 1.78±0.01ms         |    0.3  | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense')           |
| -        | 8.36±0.02ms          | 1.48±0.01ms         |    0.18 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 286±2μs              | 65.3±0.3μs          |    0.23 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 3.45±0.05ms          | 398±3μs             |    0.12 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense')              |
| -        | 3.96±0.03ms          | 275±0.9μs           |    0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 7.11±0.05ms          | 1.25±0.01ms         |    0.18 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense')              |
| -        | 8.68±0.04ms          | 826±4μs             |    0.1  | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 4.99±0.01ms          | 2.97±0.01ms         |    0.6  | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 12.2±0.1ms           | 6.53±0.03ms         |    0.53 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense')  |
| -        | 22.0±0.2ms           | 9.01±0.04ms         |    0.41 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 5.19±0.03ms          | 875±1μs             |    0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 12.8±0.09ms          | 7.04±0.03ms         |    0.55 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 102±0.9ms            | 25.6±0.04ms         |    0.25 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense')     |
| -        | 164±0.5ms            | 53.6±0.2ms          |    0.33 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 172±0.7μs            | 112±0.3μs           |    0.65 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense')      |
| -        | 359±2μs              | 221±1μs             |    0.62 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse')     |
| -        | 4.24±0.04ms          | 832±5μs             |    0.2  | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense')      |
| -        | 5.22±0.07ms          | 383±2μs             |    0.07 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 20.0±0.4ms           | 2.77±0.01ms         |    0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense')      |
| -        | 23.2±0.2ms           | 630±8μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 480±1μs              | 140±0.4μs           |    0.29 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 4.67±0.05ms          | 622±7μs             |    0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense')  |
| -        | 5.26±0.04ms          | 142±1μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 12.9±0.06ms          | 1.27±0.01ms         |    0.1  | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense')  |
| -        | 14.0±0.1ms           | 146±1μs             |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 132±1μs              | 74.1±0.5μs          |    0.56 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 249±1μs              | 88.7±1μs            |    0.36 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 24.0±0.1ms           | 10.1±0.02ms         |    0.42 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |
| -        | 28.1±0.2ms           | 15.1±0.06ms         |    0.54 | solve.TimeSparseSystem.time_linsolve_Aaug(20)                        |
| -        | 54.1±0.2ms           | 24.4±0.05ms         |    0.45 | solve.TimeSparseSystem.time_linsolve_Aaug(30)                        |
| -        | 27.9±0.1ms           | 15.1±0.06ms         |    0.54 | solve.TimeSparseSystem.time_linsolve_Ab(20)                          |
| -        | 54.2±0.2ms           | 24.4±0.2ms          |    0.45 | solve.TimeSparseSystem.time_linsolve_Ab(30)                          |

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

Patch:
diff --git a/.mailmap b/.mailmap
index 42af95d451..43b63d0271 100644
--- a/.mailmap
+++ b/.mailmap
@@ -320,6 +320,7 @@ Animesh Sinha <animeshsinha1309@gmail.com>
 Anish Shah <shah.anish07@gmail.com>
 Anjul Kumar Tyagi <anjul.ten@gmail.com>
 Ankit Agrawal <aaaagrawal@iitb.ac.in>
+Ankit Kumar Singh <ankitdiswar10@gmail.com>
 Ankit Raj Pandey <pandeyan@grinnell.edu> Ankit Pandey <pandeyan@grinnell.edu>
 Ansh Mishra <anshmishra471@gmail.com>
 Anthony Scopatz <scopatz@gmail.com>
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 36e1a26e14..0a51840296 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -531,11 +531,11 @@ def compute_density(self, expr, **kwargs):
 
         gs = solveset(expr - y, self.value, S.Reals)
 
-        if isinstance(gs, Intersection) and S.Reals in gs.args:
-            gs = list(gs.args[1])
-
-        if not gs:
-            raise ValueError("Can not solve %s for %s"%(expr, self.value))
+        if isinstance(gs, Intersection):
+            if len(gs.args) == 2 and gs.args[0] is S.Reals:
+                gs = gs.args[1]
+        if not gs.is_FiniteSet:
+            raise ValueError("Can not solve %s for %s" % (expr, self.value))
         fx = self.compute_density(self.value)
         fy = sum(fx(g) * abs(g.diff(y)) for g in gs)
         return Lambda(y, fy)


Test patch:
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index e8484a9419..b2c4206b5c 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -1577,3 +1577,7 @@ def test_issue_16318():
     # test compute_expectation function of the SingleContinuousDomain
     N = SingleContinuousDomain(x, Interval(0, 1))
     raises(ValueError, lambda: SingleContinuousDomain.compute_expectation(N, x+1, {x, y}))
+
+def test_compute_density():
+    X = Normal('X', 0, Symbol("sigma")**2)
+    raises(ValueError, lambda: density(X**5 + X))