Repository: sympy/sympy

Issue:
Incorrect simplification of expression with Piecewise
## Description
When having a sum of `Piecewise` functions in the form below, the application of `simplify` or `piecewise_fold` give incorrect results. It seems like the second term of `expr` is dropped during the simplification. Hence, evaluating the function with and without simplifications give different outputs. When subdividing `expr` into the individual summands and simplifying them separately, the result is correct.

The chosen expression might look unneccessary complicated, but I found no simpler example as removing one of the denominators results in a successful simplification.

## Minimal example

```python
import sympy
from sympy.parsing.sympy_parser import parse_expr

a = sympy.symbols('a', positive=True)

# Expression that shall be simplified
expr = parse_expr("Piecewise((1, a<=1), (0, True))+ Piecewise((1, a<=1), (0, True))/(2*a) 
                   + (Piecewise((1, a>1), (0, True)))/a",local_dict={"a":a})
# A correct simplification of `expr`
expr_correct = parse_expr("Piecewise((1 + 1/(2*a), a<=1), (1/a, True))",local_dict={"a":a})
# Evaluation of the expression at a = 1, with and without simplifications
print("correct: {}, not simplified: {}, simplified: {}, piecewisefold: {}".format(
       expr_correct.xreplace({a:1}),expr.xreplace({a:1}),expr.simplify().xreplace({a:1}),
       sympy.piecewise_fold(expr).xreplace({a:1})))
# OUTPUT:
# correct: 3/2, not simplified: 3/2, simplified: 1, piecewisefold: 1
print(expr.simplify())
# OUTPUT:
# Piecewise((1/a, a > 1), (1, True))
```

## Used version
```
>>> import sys, sympy; print(sympy.__version__, sys.version)
1.10.1 3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:38:29) [Clang 13.0.1 ] 
```

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 6b04fe9fda..6f0e24a24c 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1017,7 +1017,6 @@ Peter Enenkel <peter.enenkel+git@gmail.com> <peter.enenkel+github@gmail.com>
 Peter Schmidt <peter@peterjs.com>
 Peter Stangl <peter.stangl@ph.tum.de>
 Petr Kungurtsev <corwinat@gmail.com> Corwinpro <corwinat@gmail.com>
-Phil LeMaitre <phil_lemaitre@live.ca> <112662371+AntiPhoton47@users.noreply.github.com>
 Phil Ruffwind <rf@rufflewind.com>
 Philippe Bouafia <philippe.bouafia@ensea.fr>
 Phillip Berndt <phillip.berndt@googlemail.com>
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 93cd5540e2..9671ec42db 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -1300,9 +1300,7 @@ def _piecewise_collapse_arguments(_args):
                 newargs[-1] = ExprCondPair(expr, orcond)
                 continue
             elif newargs[-1].cond == cond:
-                newargs[-1] = ExprCondPair(expr, cond)
                 continue
-
         newargs.append(ExprCondPair(expr, cond))
     return newargs
 
diff --git a/sympy/physics/wigner.py b/sympy/physics/wigner.py
index ef404b98ac..edf71a530c 100644
--- a/sympy/physics/wigner.py
+++ b/sympy/physics/wigner.py
@@ -26,10 +26,6 @@
 .. [Liberatodebrito82] 'FORTRAN program for the integral of three
   spherical harmonics', A. Liberato de Brito,
   Comput. Phys. Commun., Volume 25, pp. 81-85 (1982)
-.. [Homeier96] 'Some Properties of the Coupling Coefficients of Real
-  Spherical Harmonics and Their Relation to Gaunt Coefficients',
-  H. H. H. Homeier and E. O. Steinborn J. Mol. Struct., Volume 368,
-  pp. 31-37 (1996)
 
 Credits and Copyright
 =====================
@@ -47,8 +43,6 @@
 
 - Oscar Gerardo Lazo Arjona (2017-06-18): added Wigner D matrices
 
-- Phil Adam LeMaitre (2022-09-19): added real Gaunt coefficient
-
 Copyright (C) 2008 Jens Rasch <jyr2000@gmail.com>
 
 """
@@ -60,14 +54,12 @@
 from sympy.core.symbol import Dummy
 from sympy.core.sympify import sympify
 from sympy.functions.combinatorial.factorials import (binomial, factorial)
-from sympy.functions.elementary.complexes import re
 from sympy.functions.elementary.exponential import exp
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (cos, sin)
 from sympy.functions.special.spherical_harmonics import Ynm
 from sympy.matrices.dense import zeros
 from sympy.matrices.immutable import ImmutableMatrix
-from sympy.utilities.misc import as_int
 
 # This list of precomputed factorials is needed to massively
 # accelerate future calculations of the various coefficients
@@ -707,32 +699,40 @@ def gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None):
 
     Jens Rasch (2009-03-24): initial version for Sage.
     """
-    l_1, l_2, l_3, m_1, m_2, m_3 = [
-        as_int(i) for i in (l_1, l_2, l_3, m_1, m_2, m_3)]
-
-    if l_1 + l_2 - l_3 < 0:
+    if int(l_1) != l_1 or int(l_2) != l_2 or int(l_3) != l_3:
+        raise ValueError("l values must be integer")
+    if int(m_1) != m_1 or int(m_2) != m_2 or int(m_3) != m_3:
+        raise ValueError("m values must be integer")
+
+    sumL = l_1 + l_2 + l_3
+    bigL = sumL // 2
+    a1 = l_1 + l_2 - l_3
+    if a1 < 0:
+        return S.Zero
+    a2 = l_1 - l_2 + l_3
+    if a2 < 0:
         return S.Zero
-    if l_1 - l_2 + l_3 < 0:
+    a3 = -l_1 + l_2 + l_3
+    if a3 < 0:
         return S.Zero
-    if -l_1 + l_2 + l_3 < 0:
+    if sumL % 2:
         return S.Zero
     if (m_1 + m_2 + m_3) != 0:
         return S.Zero
     if (abs(m_1) > l_1) or (abs(m_2) > l_2) or (abs(m_3) > l_3):
         return S.Zero
-    bigL, remL = divmod(l_1 + l_2 + l_3, 2)
-    if remL % 2:
-        return S.Zero
 
     imin = max(-l_3 + l_1 + m_2, -l_3 + l_2 - m_1, 0)
     imax = min(l_2 + m_2, l_1 - m_1, l_1 + l_2 - l_3)
 
-    _calc_factlist(max(l_1 + l_2 + l_3 + 1, imax + 1))
+    maxfact = max(l_1 + l_2 + l_3 + 1, imax + 1)
+    _calc_factlist(maxfact)
 
-    ressqrt = sqrt((2 * l_1 + 1) * (2 * l_2 + 1) * (2 * l_3 + 1) * \
+    argsqrt = (2 * l_1 + 1) * (2 * l_2 + 1) * (2 * l_3 + 1) * \
         _Factlist[l_1 - m_1] * _Factlist[l_1 + m_1] * _Factlist[l_2 - m_2] * \
         _Factlist[l_2 + m_2] * _Factlist[l_3 - m_3] * _Factlist[l_3 + m_3] / \
-        (4*pi))
+        (4*pi)
+    ressqrt = sqrt(argsqrt)
 
     prefac = Integer(_Factlist[bigL] * _Factlist[l_2 - l_1 + l_3] *
                      _Factlist[l_1 - l_2 + l_3] * _Factlist[l_1 + l_2 - l_3])/ \
@@ -753,160 +753,6 @@ def gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None):
     return res
 
 
-def real_gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None):
-    r"""
-    Calculate the real Gaunt coefficient.
-
-    Explanation
-    ===========
-
-    The real Gaunt coefficient is defined as the integral over three
-    real spherical harmonics:
-
-    .. math::
-        \begin{aligned}
-        \operatorname{RealGaunt}(l_1,l_2,l_3,m_1,m_2,m_3)
-        &=\int Z^{m_1}_{l_1}(\Omega)
-         Z^{m_2}_{l_2}(\Omega) Z^{m_3}_{l_3}(\Omega) \,d\Omega \\
-        \end{aligned}
-
-    Alternatively, it can be defined in terms of the standard Gaunt
-    coefficient by relating the real spherical harmonics to the standard
-    spherical harmonics via a unitary transformation `U`, i.e.
-    `Z^{m}_{l}(\Omega)=\sum_{m'}U^{m}_{m'}Y^{m'}_{l}(\Omega)` [Homeier96]_.
-    The real Gaunt coefficient is then defined as
-
-    .. math::
-        \begin{aligned}
-        \operatorname{RealGaunt}(l_1,l_2,l_3,m_1,m_2,m_3)
-        &=\int Z^{m_1}_{l_1}(\Omega)
-         Z^{m_2}_{l_2}(\Omega) Z^{m_3}_{l_3}(\Omega) \,d\Omega \\
-        &=\sum_{m'_1 m'_2 m'_3} U^{m_1}_{m'_1}U^{m_2}_{m'_2}U^{m_3}_{m'_3}
-         \operatorname{Gaunt}(l_1,l_2,l_3,m'_1,m'_2,m'_3)
-        \end{aligned}
-
-    The unitary matrix `U` has components
-
-    .. math::
-        \begin{aligned}
-        U^m_{m'} = \delta_{|m||m'|}*(\delta_{m'0}\delta_{m0} + \frac{1}{\sqrt{2}}\big[\Theta(m)
-        \big(\delta_{m'm}+(-1)^{m'}\delta_{m'-m}\big)+i\Theta(-m)\big((-1)^{-m}
-        \delta_{m'-m}-\delta_{m'm}*(-1)^{m'-m}\big)\big])
-        \end{aligned}
-
-    where `\delta_{ij}` is the Kronecker delta symbol and `\Theta` is a step
-    function defined as
-
-    .. math::
-        \begin{aligned}
-        \Theta(x) = \begin{cases} 1 \,\text{for}\, x > 0 \\ 0 \,\text{for}\, x \leq 0 \end{cases}
-        \end{aligned}
-
-    Parameters
-    ==========
-
-    l_1, l_2, l_3, m_1, m_2, m_3 :
-        Integer.
-
-    prec - precision, default: ``None``.
-        Providing a precision can
-        drastically speed up the calculation.
-
-    Returns
-    =======
-
-    Rational number times the square root of a rational number.
-
-    Examples
-    ========
-
-    >>> from sympy.physics.wigner import real_gaunt
-    >>> real_gaunt(2,2,4,-1,-1,0)
-    -2/(7*sqrt(pi))
-    >>> real_gaunt(10,10,20,-9,-9,0).n(64)
-    -0.00002480019791932209313156167...
-
-    It is an error to use non-integer values for `l` and `m`::
-        real_gaunt(2.8,0.5,1.3,0,0,0)
-        Traceback (most recent call last):
-        ...
-        ValueError: l values must be integer
-        real_gaunt(2,2,4,0.7,1,-3.4)
-        Traceback (most recent call last):
-        ...
-        ValueError: m values must be integer
-
-    Notes
-    =====
-
-    The real Gaunt coefficient inherits from the standard Gaunt coefficient,
-    the invariance under any permutation of the pairs `(l_i, m_i)` and the
-    requirement that the sum of the `l_i` be even to yield a non-zero value.
-    It also obeys the following symmetry rules:
-
-    - zero for `l_1`, `l_2`, `l_3` not fulfiling the condition
-      `l_1 \in \{l_{\text{max}}, l_{\text{max}}-2, \ldots, l_{\text{min}}\}`,
-      where `l_{\text{max}} = l_2+l_3`,
-
-      .. math::
-          \begin{aligned}
-          l_{\text{min}} = \begin{cases} \kappa(l_2, l_3, m_2, m_3) & \text{if}\,
-          \kappa(l_2, l_3, m_2, m_3) + l_{\text{max}}\, \text{is even} \\
-          \kappa(l_2, l_3, m_2, m_3)+1 & \text{if}\, \kappa(l_2, l_3, m_2, m_3) +
-          l_{\text{max}}\, \text{is odd}\end{cases}
-          \end{aligned}
-
-      and `\kappa(l_2, l_3, m_2, m_3) = \max{\big(|l_2-l_3|, \min{\big(|m_2+m_3|,
-      |m_2-m_3|\big)}\big)}`
-
-    - zero for an odd number of negative `m_i`
-
-    Algorithms
-    ==========
-
-    This function uses the algorithms of [Homeier96]_ and [Rasch03]_ to
-    calculate the value of the real Gaunt coefficient exactly. Note that
-    the formula used in [Rasch03]_ contains alternating sums over large
-    factorials and is therefore unsuitable for finite precision arithmetic
-    and only useful for a computer algebra system [Rasch03]_. However, this
-    function can in principle use any algorithm that computes the Gaunt
-    coefficient, so it is suitable for finite precision arithmetic in so far
-    as the algorithm which computes the Gaunt coefficient is.
-    """
-    l_1, l_2, l_3, m_1, m_2, m_3 = [
-        as_int(i) for i in (l_1, l_2, l_3, m_1, m_2, m_3)]
-
-    # check for quick exits
-    if sum(1 for i in (m_1, m_2, m_3) if i < 0) % 2:
-        return S.Zero  # odd number of negative m
-    if (l_1 + l_2 + l_3) % 2:
-        return S.Zero  # sum of l is odd
-    lmax = l_2 + l_3
-    lmin = max(abs(l_2 - l_3), min(abs(m_2 + m_3), abs(m_2 - m_3)))
-    if (lmin + lmax) % 2:
-        lmin += 1
-    if lmin not in range(lmax, lmin - 2, -2):
-        return S.Zero
-
-    kron_del = lambda i, j: 1 if i == j else 0
-    s = lambda e: -1 if e % 2 else 1  #  (-1)**e to give +/-1, avoiding float when e<0
-    A = lambda a, b: (-kron_del(a, b)*s(a-b) + kron_del(a, -b)*
-                      s(b)) if b < 0 else 0
-    B = lambda a, b: (kron_del(a, b) + kron_del(a, -b)*s(a)) if b > 0 else 0
-    C = lambda a, b: kron_del(abs(a), abs(b))*(kron_del(a, 0)*kron_del(b, 0) +
-                                          (B(a, b) + I*A(a, b))/sqrt(2))
-    ugnt = 0
-    for i in range(-l_1, l_1+1):
-        U1 = C(i, m_1)
-        for j in range(-l_2, l_2+1):
-            U2 = C(j, m_2)
-            U3 = C(-i-j, m_3)
-            ugnt = ugnt + re(U1*U2*U3)*gaunt(l_1, l_2, l_3, i, j, -i-j)
-
-    if prec is not None:
-        ugnt = ugnt.n(prec)
-    return ugnt
-
 
 class Wigner3j(Function):
 


Test patch:
diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py
index d7b202d03b..ee6a044b5b 100644
--- a/sympy/functions/elementary/tests/test_piecewise.py
+++ b/sympy/functions/elementary/tests/test_piecewise.py
@@ -1548,3 +1548,8 @@ def test_issue_22533():
     x = Symbol('x', real=True)
     f = Piecewise((-1 / x, x <= 0), (1 / x, True))
     assert integrate(f, x) == Piecewise((-log(x), x <= 0), (log(x), True))
+
+
+def test_issue_24072():
+    assert Piecewise((1, x > 1), (2, x <= 1), (3, x <= 1)
+        ) == Piecewise((1, x > 1), (2, True))
diff --git a/sympy/physics/tests/test_clebsch_gordan.py b/sympy/physics/tests/test_clebsch_gordan.py
index d6a41f1b6b..2e77bf2f7a 100644
--- a/sympy/physics/tests/test_clebsch_gordan.py
+++ b/sympy/physics/tests/test_clebsch_gordan.py
@@ -1,4 +1,4 @@
-from sympy.core.numbers import (I, pi, Rational)
+from sympy.core.numbers import (I, pi)
 from sympy.core.singleton import S
 from sympy.core.symbol import symbols
 from sympy.functions.elementary.exponential import exp
@@ -7,8 +7,8 @@
 from sympy.functions.special.spherical_harmonics import Ynm
 from sympy.matrices.dense import Matrix
 from sympy.physics.wigner import (clebsch_gordan, wigner_9j, wigner_6j, gaunt,
-        real_gaunt, racah, dot_rot_grad_Ynm, wigner_3j, wigner_d_small, wigner_d)
-from sympy.testing.pytest import raises
+        racah, dot_rot_grad_Ynm, wigner_3j, wigner_d_small, wigner_d)
+from sympy.core.numbers import Rational
 
 # for test cases, refer : https://en.wikipedia.org/wiki/Table_of_Clebsch%E2%80%93Gordan_coefficients
 
@@ -18,41 +18,223 @@ def test_clebsch_gordan_docs():
     assert clebsch_gordan(Rational(3, 2), S.Half, 1, Rational(-1, 2), S.Half, 0) == -sqrt(2)/2
 
 
-def test_clebsch_gordan():
-    # Argument order: (j_1, j_2, j, m_1, m_2, m)
-
-    h = S.One
-    k = S.Half
-    l = Rational(3, 2)
-    i = Rational(-1, 2)
-    n = Rational(7, 2)
-    p = Rational(5, 2)
-    assert clebsch_gordan(k, k, 1, k, k, 1) == 1
-    assert clebsch_gordan(k, k, 1, k, k, 0) == 0
-    assert clebsch_gordan(k, k, 1, i, i, -1) == 1
-    assert clebsch_gordan(k, k, 1, k, i, 0) == sqrt(2)/2
-    assert clebsch_gordan(k, k, 0, k, i, 0) == sqrt(2)/2
-    assert clebsch_gordan(k, k, 1, i, k, 0) == sqrt(2)/2
-    assert clebsch_gordan(k, k, 0, i, k, 0) == -sqrt(2)/2
-    assert clebsch_gordan(h, k, l, 1, k, l) == 1
-    assert clebsch_gordan(h, k, l, 1, i, k) == 1/sqrt(3)
-    assert clebsch_gordan(h, k, k, 1, i, k) == sqrt(2)/sqrt(3)
-    assert clebsch_gordan(h, k, k, 0, k, k) == -1/sqrt(3)
-    assert clebsch_gordan(h, k, l, 0, k, k) == sqrt(2)/sqrt(3)
-    assert clebsch_gordan(h, h, S(2), 1, 1, S(2)) == 1
-    assert clebsch_gordan(h, h, S(2), 1, 0, 1) == 1/sqrt(2)
-    assert clebsch_gordan(h, h, S(2), 0, 1, 1) == 1/sqrt(2)
-    assert clebsch_gordan(h, h, 1, 1, 0, 1) == 1/sqrt(2)
-    assert clebsch_gordan(h, h, 1, 0, 1, 1) == -1/sqrt(2)
-    assert clebsch_gordan(l, l, S(3), l, l, S(3)) == 1
-    assert clebsch_gordan(l, l, S(2), l, k, S(2)) == 1/sqrt(2)
-    assert clebsch_gordan(l, l, S(3), l, k, S(2)) == 1/sqrt(2)
-    assert clebsch_gordan(S(2), S(2), S(4), S(2), S(2), S(4)) == 1
-    assert clebsch_gordan(S(2), S(2), S(3), S(2), 1, S(3)) == 1/sqrt(2)
-    assert clebsch_gordan(S(2), S(2), S(3), 1, 1, S(2)) == 0
-    assert clebsch_gordan(p, h, n, p, 1, n) == 1
-    assert clebsch_gordan(p, h, p, p, 0, p) == sqrt(5)/sqrt(7)
-    assert clebsch_gordan(p, h, l, k, 1, l) == 1/sqrt(15)
+def test_clebsch_gordan1():
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 1
+    j = 1
+    m_1 = S.Half
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = -1
+    j = 1
+    m_1 = Rational(-1, 2)
+    m_2 = Rational(-1, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 0
+    j = 1
+    m_1 = S.Half
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 0
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 0
+    j = 1
+    m_1 = S.Half
+    m_2 = Rational(-1, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(2)/2
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 0
+    j = 0
+    m_1 = S.Half
+    m_2 = Rational(-1, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(2)/2
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 0
+    j = 1
+    m_1 = Rational(-1, 2)
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(2)/2
+
+    j_1 = S.Half
+    j_2 = S.Half
+    m = 0
+    j = 0
+    m_1 = Rational(-1, 2)
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == -sqrt(2)/2
+
+def test_clebsch_gordan2():
+    j_1 = S.One
+    j_2 = S.Half
+    m = Rational(3, 2)
+    j = Rational(3, 2)
+    m_1 = 1
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+    j_1 = S.One
+    j_2 = S.Half
+    m = S.Half
+    j = Rational(3, 2)
+    m_1 = 1
+    m_2 = Rational(-1, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(3)
+
+    j_1 = S.One
+    j_2 = S.Half
+    m = S.Half
+    j = S.Half
+    m_1 = 1
+    m_2 = Rational(-1, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(2)/sqrt(3)
+
+    j_1 = S.One
+    j_2 = S.Half
+    m = S.Half
+    j = S.Half
+    m_1 = 0
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == -1/sqrt(3)
+
+    j_1 = S.One
+    j_2 = S.Half
+    m = S.Half
+    j = Rational(3, 2)
+    m_1 = 0
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(2)/sqrt(3)
+
+    j_1 = S.One
+    j_2 = S.One
+    m = S(2)
+    j = S(2)
+    m_1 = 1
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+
+    j_1 = S.One
+    j_2 = S.One
+    m = 1
+    j = S(2)
+    m_1 = 1
+    m_2 = 0
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+
+    j_1 = S.One
+    j_2 = S.One
+    m = 1
+    j = S(2)
+    m_1 = 0
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+    j_1 = S.One
+    j_2 = S.One
+    m = 1
+    j = 1
+    m_1 = 1
+    m_2 = 0
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+    j_1 = S.One
+    j_2 = S.One
+    m = 1
+    j = 1
+    m_1 = 0
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == -1/sqrt(2)
+
+def test_clebsch_gordan3():
+    j_1 = Rational(3, 2)
+    j_2 = Rational(3, 2)
+    m = S(3)
+    j = S(3)
+    m_1 = Rational(3, 2)
+    m_2 = Rational(3, 2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+
+    j_1 = Rational(3, 2)
+    j_2 = Rational(3, 2)
+    m = S(2)
+    j = S(2)
+    m_1 = Rational(3, 2)
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+    j_1 = Rational(3, 2)
+    j_2 = Rational(3, 2)
+    m = S(2)
+    j = S(3)
+    m_1 = Rational(3, 2)
+    m_2 = S.Half
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+def test_clebsch_gordan4():
+    j_1 = S(2)
+    j_2 = S(2)
+    m = S(4)
+    j = S(4)
+    m_1 = S(2)
+    m_2 = S(2)
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+
+    j_1 = S(2)
+    j_2 = S(2)
+    m = S(3)
+    j = S(3)
+    m_1 = S(2)
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
+
+    j_1 = S(2)
+    j_2 = S(2)
+    m = S(2)
+    j = S(3)
+    m_1 = 1
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 0
+
+def test_clebsch_gordan5():
+    j_1 = Rational(5, 2)
+    j_2 = S.One
+    m = Rational(7, 2)
+    j = Rational(7, 2)
+    m_1 = Rational(5, 2)
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
+
+
+    j_1 = Rational(5, 2)
+    j_2 = S.One
+    m = Rational(5, 2)
+    j = Rational(5, 2)
+    m_1 = Rational(5, 2)
+    m_2 = 0
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(5)/sqrt(7)
+
+    j_1 = Rational(5, 2)
+    j_2 = S.One
+    m = Rational(3, 2)
+    j = Rational(3, 2)
+    m_1 = S.Half
+    m_2 = 1
+    assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(15)
 
 
 def test_wigner():
@@ -112,40 +294,6 @@ def gaunt_ref(l1, l2, l3, m1, m2, m3):
                                 assert abs(g) < threshold
                             if (l1 + l2 + l3) % 2:
                                 assert abs(g) < threshold
-    assert gaunt(1, 1, 0, 0, 2, -2) is S.Zero
-
-
-def test_realgaunt():
-    # All non-zero values corresponding to l values from 0 to 2
-    for l in range(3):
-        for m in range(-l, l+1):
-            assert real_gaunt(0, l, l, 0, m, m) == 1/(2*sqrt(pi))
-    assert real_gaunt(1, 1, 2, 0, 0, 0) == sqrt(5)/(5*sqrt(pi))
-    assert real_gaunt(1, 1, 2, 1, 1, 0) == -sqrt(5)/(10*sqrt(pi))
-    assert real_gaunt(2, 2, 2, 0, 0, 0) == sqrt(5)/(7*sqrt(pi))
-    assert real_gaunt(2, 2, 2, 0, 2, 2) == -sqrt(5)/(7*sqrt(pi))
-    assert real_gaunt(2, 2, 2, -2, -2, 0) == -sqrt(5)/(7*sqrt(pi))
-    assert real_gaunt(1, 1, 2, -1, 0, -1) == sqrt(15)/(10*sqrt(pi))
-    assert real_gaunt(1, 1, 2, 0, 1, 1) == sqrt(15)/(10*sqrt(pi))
-    assert real_gaunt(1, 1, 2, 1, 1, 2) == sqrt(15)/(10*sqrt(pi))
-    assert real_gaunt(1, 1, 2, -1, 1, -2) == -sqrt(15)/(10*sqrt(pi))
-    assert real_gaunt(1, 1, 2, -1, -1, 2) == -sqrt(15)/(10*sqrt(pi))
-    assert real_gaunt(2, 2, 2, 0, 1, 1) == sqrt(5)/(14*sqrt(pi))
-    assert real_gaunt(2, 2, 2, 1, 1, 2) == sqrt(15)/(14*sqrt(pi))
-    assert real_gaunt(2, 2, 2, -1, -1, 2) == -sqrt(15)/(14*sqrt(pi))
-
-    assert real_gaunt(-2, -2, -2, -2, -2, 0) is S.Zero  # m test
-    assert real_gaunt(-2, 1, 0, 1, 1, 1) is S.Zero  # l test
-    assert real_gaunt(-2, -1, -2, -1, -1, 0) is S.Zero  # m and l test
-    assert real_gaunt(-2, -2, -2, -2, -2, -2) is S.Zero  # m and k test
-    assert real_gaunt(-2, -1, -2, -1, -1, -1) is S.Zero  # m, l and k test
-
-    x = symbols('x', integer=True)
-    v = [0]*6
-    for i in range(len(v)):
-        v[i] = x  # non literal ints fail
-        raises(ValueError, lambda: real_gaunt(*v))
-        v[i] = 0
 
 
 def test_racah():