Repository: sympy/sympy
Issue category: Bug Report

Issue:
wigner_6j and wigner_9j violating numpy behavioral subtyping
Dear Developers,

I get an error when I run the following code snippet using wigner_6j.

```
import numpy as np
from sympy.physics.wigner import wigner_6j

y = np.float64(5.0)
x = wigner_6j(5, 5, y, 5, 5, 5)
```


This produces the following error message:

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 5
      2 from sympy.physics.wigner import wigner_6j
      4 y = np.float64(5.0)
----> 5 x = wigner_6j(5, 5, y, 5, 5, 5)

File ~\code\.venv\jlab\lib\site-packages\sympy\physics\wigner.py:572, in wigner_6j(j_1, j_2, j_3, j_4, j_5, j_6, prec)
    479 def wigner_6j(j_1, j_2, j_3, j_4, j_5, j_6, prec=None):
    480     r"""
    481     Calculate the Wigner 6j symbol `\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)`.
    482 
   (...)
    569 
    570     """
    571     res = (-1) ** int(j_1 + j_2 + j_4 + j_5) * \
--> 572         racah(j_1, j_2, j_5, j_4, j_3, j_6, prec)
    573     return res

File ~\code\.venv\jlab\lib\site-packages\sympy\physics\wigner.py:451, in racah(aa, bb, cc, dd, ee, ff, prec)
    398 def racah(aa, bb, cc, dd, ee, ff, prec=None):
    399     r"""
    400     Calculate the Racah symbol `W(a,b,c,d;e,f)`.
    401 
   (...)
    449     - Jens Rasch (2009-03-24): initial version
    450     """
--> 451     prefac = _big_delta_coeff(aa, bb, ee, prec) * \
    452         _big_delta_coeff(cc, dd, ee, prec) * \
    453         _big_delta_coeff(aa, cc, ff, prec) * \
    454         _big_delta_coeff(bb, dd, ff, prec)
    455     if prefac == 0:
    456         return S.Zero

File ~\code\.venv\jlab\lib\site-packages\sympy\physics\wigner.py:372, in _big_delta_coeff(aa, bb, cc, prec)
    369 # the triangle test will only pass if a) all 3 values are ints or
    370 # b) 1 is an int and the other two are half-ints
    371 if not int_valued(aa + bb - cc):
--> 372     raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
    373 if not int_valued(aa + cc - bb):
    374     raise ValueError("j values must be integer or half integer and fulfill the triangle relation")

ValueError: j values must be integer or half integer and fulfill the triangle relation
```

A similar issue is found for wigner_9j. A similar issue was found for the wigner_3j method in #27219, and was fixed in [#27228](https://github.com/sympy/sympy/pull/27228).

As an aside maybe there could be some underlying shared function to fix wigner_3j, wigner_6j, wigner_9j, clebsch_gordan, etc., Instead of repeatedly patching.


Hints:

Patch:
diff --git a/sympy/physics/wigner.py b/sympy/physics/wigner.py
index 3dcd2f8f4690..5874acbb61c8 100644
--- a/sympy/physics/wigner.py
+++ b/sympy/physics/wigner.py
@@ -571,6 +571,8 @@ def wigner_6j(j_1, j_2, j_3, j_4, j_5, j_6, prec=None):
     algebra system [Rasch03]_.
 
     """
+    j_1, j_2, j_3, j_4, j_5, j_6 = map(sympify, \
+                [j_1, j_2, j_3, j_4, j_5, j_6])
     res = (-1) ** int(j_1 + j_2 + j_4 + j_5) * \
         racah(j_1, j_2, j_5, j_4, j_3, j_6, prec)
     return res
@@ -627,6 +629,8 @@ def wigner_9j(j_1, j_2, j_3, j_4, j_5, j_6, j_7, j_8, j_9, prec=None):
     for finite precision arithmetic and only useful for a computer
     algebra system [Rasch03]_.
     """
+    j_1, j_2, j_3, j_4, j_5, j_6, j_7, j_8, j_9 = map(sympify, \
+                [j_1, j_2, j_3, j_4, j_5, j_6, j_7, j_8, j_9])
     imax = int(min(j_1 + j_9, j_2 + j_6, j_4 + j_8) * 2)
     imin = imax % 2
     sumres = 0


Test patch:
diff --git a/sympy/physics/tests/test_clebsch_gordan.py b/sympy/physics/tests/test_clebsch_gordan.py
index 2085d4a13479..e4313e3e412d 100644
--- a/sympy/physics/tests/test_clebsch_gordan.py
+++ b/sympy/physics/tests/test_clebsch_gordan.py
@@ -67,6 +67,10 @@ def test_clebsch_gordan_numpy():
 
 
 def test_wigner():
+    try:
+        import numpy as np
+    except ImportError:
+        skip("numpy not installed")
     def tn(a, b):
         return (a - b).n(64) < S('1e-64')
     assert tn(wigner_9j(1, 1, 1, 1, 1, 1, 1, 1, 0, prec=64), Rational(1, 18))
@@ -74,6 +78,8 @@ def tn(a, b):
         70)/(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105))
     assert wigner_6j(5, 5, 5, 5, 5, 5) == Rational(1, 52)
     assert tn(wigner_6j(8, 8, 8, 8, 8, 8, prec=64), Rational(-12219, 965770))
+    assert wigner_6j(1, 1, 1, 1.0, np.float64(1.0), 1) == Rational(1, 6)
+    assert wigner_6j(3.0, np.float32(3), 3.0, 3, 3, 3) == Rational(-1, 14)
     # regression test for #8747
     half = S.Half
     assert wigner_9j(0, 0, 0, 0, half, half, 0, half, half) == half
@@ -89,6 +95,14 @@ def tn(a, b):
                       2, 4, 4,
                       5, 7 * half, 7 * half)
             == -sqrt(Rational(3481, 5042614500)))
+    assert (wigner_9j(5, 5, 5.0,
+                      np.float64(5.0), 5, 5,
+                      5, 5, 5)
+            == 0)
+    assert (wigner_9j(1.0, 2.0, 3.0,
+                      3, 2, 1,
+                      2, 1, 3)
+            == -4*sqrt(70)/11025)
 
 
 def test_gaunt():