factorial for large numbers (#628)

* impl [apps/similarity] - orbit_cardinality

factorials of numbers larger than 170 do not fit into a int, hence return float using the qarg `exact=True` when needed

* format with black

* refactor [apps] - similarity if/else

remove unneeded `else` statement

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update CHANGELOG.md

* test [apps] - similarity

test case for large numbers

* format [tests] - black

Co-authored-by: Josh Izaac <josh146@gmail.com>

.github/CHANGELOG.md

@@ -146,7 +146,7 @@
 
 
 * Added native support in the Fock backend for the MZgate.
-  [#610](https://github.com/XanaduAI/strawberryfields/issues/610)
+  [(#610)](https://github.com/XanaduAI/strawberryfields/issues/610)
 
 * `measure_threshold` is now supported in the `bosonic` backend.
   [(#618)](https://github.com/XanaduAI/strawberryfields/pull/618)
@@ -175,6 +175,9 @@
 * Measurement arguments are now stored as non-keyword arguments, instead of keyword arguments, in
   the resulting Blackbird program when using the `io.to_blackbird()` converter function.
   [(#622)](https://github.com/XanaduAI/strawberryfields/pull/622)
+
+* Now factorials of numbers larger than 170 are calculated approximately and stored on a `float`.
+  [(#628)](https://github.com/XanaduAI/strawberryfields/pull/628)
 
 <h3>Documentation</h3>
 

strawberryfields/apps/similarity.py

@@ -284,7 +284,7 @@ def orbits(photon_number: int) -> Generator[list, None, None]:
         yield sorted(a[: k + 1], reverse=True)
 
 
-def orbit_cardinality(orbit: list, modes: int) -> int:
+def orbit_cardinality(orbit: list, modes: int) -> Union[int, float]:
     """Gives the number of samples belonging to the input orbit.
 
     For example, there are three possible samples in the orbit [2, 1, 1] with three modes: [1, 1,
@@ -304,7 +304,13 @@ def orbit_cardinality(orbit: list, modes: int) -> int:
     """
     sample = orbit + [0] * (modes - len(orbit))
     counts = list(Counter(sample).values())
-    return int(factorial(modes) / np.prod(factorial(counts)))
+
+    # factorials of numbers larger than 170 do not fit into a int,
+    # hence return float using the qarg `exact=True`
+    if modes > 170:
+        return factorial(modes, exact=True) / np.prod(factorial(counts, exact=True))
+
+    return int(factorial(modes, exact=False) / np.prod(factorial(counts, exact=False)))
 
 
 def event_cardinality(photon_number: int, max_count_per_mode: int, modes: int) -> int:

tests/apps/test_similarity.py

@@ -194,6 +194,7 @@ def test_sample_max_count(self, photon_num, modes_dim, count):
     [(1, 1, 1, 1), 5, 5],
     [(1, 1, 2), 5, 30],
     [(1, 2, 3), 5, 60],
+    [(2, 1), 171, 29070.0],
 ]
 
 
@@ -202,7 +203,7 @@ def test_orbit_cardinality(orbit, max_photon, expected):
     """Test if function ``strawberryfields.apps.similarity.orbit_cardinality`` returns the
     correct number of samples for some hard-coded examples."""
 
-    assert similarity.orbit_cardinality(list(orbit), max_photon) == expected
+    assert np.allclose(similarity.orbit_cardinality(list(orbit), max_photon), expected)
 
 
 events = [