use scipy double factorial

Author: lazyoracle

src/isotropic/e2.py

@@ -6,9 +6,28 @@
 import jax.random as random
 from jax import Array
 from jax.typing import ArrayLike
+from scipy.special import factorial2
 
 from isotropic.utils.bisection import get_theta
-from isotropic.utils.distribution import double_factorial
+
+
+def double_factorial_ratio(n: int, m: int) -> float:
+    """
+    Compute the ratio of double factorials n!! / m!!.
+
+    Parameters
+    ----------
+    n : int
+        The numerator double factorial.
+    m : int
+        The denominator double factorial.
+
+    Returns
+    -------
+    float
+        The ratio of the double factorials.
+    """
+    return factorial2(n) / factorial2(m)
 
 
 def F_j(theta_j: float, j: int, d: int) -> Array:
@@ -30,12 +49,11 @@ def F_j(theta_j: float, j: int, d: int) -> Array:
         The value of the function F_j evaluated at theta_j.
     """
     dj = d - j
-    num = double_factorial(dj - 2)
-    den = double_factorial(dj - 1)
+    numoverden = double_factorial_ratio(dj - 2, dj - 1)
 
     def F_odd(_):
-        C_j = double_factorial(dj - 1) / (2 * double_factorial(dj - 2))
-        prefactor = C_j * num / den
+        C_j = (1 / 2) * double_factorial_ratio(dj - 1, dj - 2)
+        prefactor = C_j * numoverden
         k_max = (dj - 2) // 2  # upper bound for k in range
         k_vals = jnp.arange(0, k_max + 1)
 
@@ -54,8 +72,8 @@ def product_term(k):
         return prefactor - C_j * jnp.cos(theta_j) * sum_terms
 
     def F_even(_):
-        C_j = double_factorial(dj - 1) / (jnp.pi * double_factorial(dj - 2))
-        prefactor = C_j * (num / den) * theta_j
+        C_j = (1 / jnp.pi) * double_factorial_ratio(dj - 1, dj - 2)
+        prefactor = C_j * numoverden * theta_j
         k_max = (dj - 1) // 2
         k_vals = jnp.arange(1, k_max + 1)
 

src/isotropic/utils/distribution.py

@@ -46,6 +46,8 @@ def double_factorial_ratio(num: int, den: int) -> Array:
     -----
     For very large numbers, this is numerically stable only when |num - den| is ~5.
     """
+    if abs(num - den) > 4:
+        raise ValueError("num and den should be close to each other")
     num_elems = jnp.arange(num, 0, -2, dtype=jnp.uint64)
     den_elems = jnp.arange(den, 0, -2, dtype=jnp.uint64)
 

tests/test_e2.py

@@ -1,8 +1,8 @@
 import jax.numpy as jnp
 from jax import Array, random
+from scipy.special import factorial2
 
 from isotropic.e2 import F_j, get_e2_coeffs
-from isotropic.utils.distribution import double_factorial
 
 
 def test_get_e2():
@@ -41,9 +41,9 @@ def test_F_j_even():
     result = F_j(theta_j, j, d)
 
     # manually calculate expected result
-    C_j = double_factorial(d - j - 1) / (jnp.pi * double_factorial(d - j - 2))
-    num = double_factorial(d - j - 2)
-    den = double_factorial(d - j - 1)
+    C_j = factorial2(d - j - 1) / (jnp.pi * factorial2(d - j - 2))
+    num = factorial2(d - j - 2)
+    den = factorial2(d - j - 1)
     prefactor = C_j * (num / den) * theta_j
 
     # k goes from 1 to (d - j - 1) // 2, i.e., 1 to 2
@@ -65,9 +65,9 @@ def test_F_j_odd():
     result = F_j(theta_j, j, d)
 
     # manually calculate expected result
-    C_j = double_factorial(d - j - 1) / (2 * double_factorial(d - j - 2))
-    num = double_factorial(d - j - 2)
-    den = double_factorial(d - j - 1)
+    C_j = factorial2(d - j - 1) / (2 * factorial2(d - j - 2))
+    num = factorial2(d - j - 2)
+    den = factorial2(d - j - 1)
     prefactor = C_j * num / den
 
     # k goes from 0 to (d - j - 2) // 2, i.e., 0 to 2

tests/test_error_generation.py

@@ -13,10 +13,9 @@
 
 
 def test_add_isotropic_error():
-    Phi_original = jnp.asarray([1 + 0j, 1 + 0j], dtype=complex) / jnp.sqrt(
-        2
-    )  # n = 1, d = 3
-    Phi_spherical = statevector_to_hypersphere(Phi_original)  # d+1 = 4
+    Phi = jnp.ones(4, dtype=complex)
+    Phi = Phi / jnp.linalg.norm(Phi)
+    Phi_spherical = statevector_to_hypersphere(Phi)  # d+1 = 4
     basis = get_orthonormal_basis(
         Phi_spherical
     )  # gives d vectors with d+1 elements each
@@ -27,8 +26,13 @@ def test_add_isotropic_error():
     )
     e2 = jnp.expand_dims(coeffs, axis=-1) * basis
 
+    # orthogonality check
+    assert jnp.allclose(jnp.abs(jnp.dot(e2, Phi_spherical)), 0.0), (
+        f"Expected 0.0, got {jnp.abs(jnp.dot(basis, Phi_spherical))}"
+    )
+
     def g(theta):
-        return normal_integrand(theta, d=Phi_spherical.shape[0], sigma=0.9)
+        return normal_integrand(theta, d=Phi_spherical.shape[0], sigma=0.6)
 
     theta_zero = get_theta_zero(x=0.5, g=g)
     Psi_spherical = add_isotropic_error(Phi_spherical, e2=e2, theta_zero=theta_zero)
@@ -38,7 +42,3 @@ def g(theta):
     assert jnp.isclose(jnp.linalg.norm(Psi), 1.0), (
         f"Expected 1.0, got {jnp.linalg.norm(Psi)}"
     )
-
-    # fidelity check
-    fidelity = jnp.abs(jnp.vdot(Phi_original, Psi)) ** 2
-    assert 0.0 <= fidelity <= 1.0, f"Expected fidelity in [0, 1], got {fidelity}"

tests/test_utils.py

@@ -1,5 +1,6 @@
 import jax.numpy as jnp
 from scipy.linalg import null_space
+from scipy.special import factorial2
 
 from isotropic.utils.bisection import get_theta
 from isotropic.utils.distribution import (
@@ -61,39 +62,39 @@ def test_double_factorial():
     # Test even double factorial
     n_even = 6
     result_even = double_factorial(n_even)
-    expected_even = 48.0  # 6!! = 6 * 4 * 2 = 48
+    expected_even = factorial2(n_even)
     assert jnp.isclose(result_even, expected_even), (
         f"Expected {expected_even}, got {result_even}"
     )
 
     # Test odd double factorial
     n_odd = 5
     result_odd = double_factorial(n_odd)
-    expected_odd = 15.0  # 5!! = 5 * 3 * 1 = 15
+    expected_odd = factorial2(n_odd)
     assert jnp.isclose(result_odd, expected_odd), (
         f"Expected {expected_odd}, got {result_odd}"
     )
 
     # Test zero double factorial
     n_zero = 0
     result_zero = double_factorial(n_zero)
-    expected_zero = 1.0
+    expected_zero = factorial2(n_zero)
     assert jnp.isclose(result_zero, expected_zero), (
         f"Expected {expected_zero}, got {result_zero}"
     )
 
 
 def test_double_factorial_ratio():
-    num, den = (2**5) - 1, (2**5) - 2
+    num, den = (2**8) - 1, (2**8) - 2
     ratio_received = double_factorial_ratio(num, den)
-    ratio_expected = double_factorial(num) / double_factorial(den)
+    ratio_expected = factorial2(num) / factorial2(den)
     assert jnp.isclose(ratio_received, ratio_expected), (
         f"Expected {ratio_expected}, got {ratio_received}"
     )
 
-    num, den = (2**5) - 3, (2**5) - 1
+    num, den = (2**8) - 3, (2**8) - 1
     ratio_received = double_factorial_ratio(num, den)
-    ratio_expected = double_factorial(num) / double_factorial(den)
+    ratio_expected = factorial2(num) / factorial2(den)
     assert jnp.isclose(ratio_received, ratio_expected), (
         f"Expected {ratio_expected}, got {ratio_received}"
     )