Updated calls to np.prod to use dtype=np.int64. (quantumlib#4363)

Fixes quantumlib#4346 .

Liberally upgraded a lot of the calls to `np.prod` to use ~~`dtype=int`~~ `dtype=np.int64` that didn't already have it. I opted to keep the calls to `np.prod` instead of doing something like `reduce(mul, <stuff>)` because using basic numpy, we can only really comfortably create and index arrays with sizes up to the max value of a C type `long`.

~~Confirmed `dtype=int` will give us that:~~
~~https://numpy.org/doc/stable/reference/arrays.dtypes.html#specifying-and-constructing-data-types~~
~~https://numpy.org/doc/stable/reference/arrays.scalars.html#numpy.int_~~

cirq/circuits/circuit.py

@@ -1000,7 +1000,7 @@ def unitary(
 
         # Force qubits to have dimension at least 2 for backwards compatibility.
         qid_shape = self.qid_shape(qubit_order=qs)
-        side_len = np.product(qid_shape, dtype=int)
+        side_len = np.prod(qid_shape, dtype=np.int64)
 
         state = qis.eye_tensor(qid_shape, dtype=dtype)
 
@@ -1083,7 +1083,7 @@ def final_state_vector(
 
         # Force qubits to have dimension at least 2 for backwards compatibility.
         qid_shape = self.qid_shape(qubit_order=qs)
-        state_len = np.product(qid_shape, dtype=int)
+        state_len = np.prod(qid_shape, dtype=np.int64)
 
         state = qis.to_valid_state_vector(initial_state, qid_shape=qid_shape, dtype=dtype).reshape(
             qid_shape

cirq/experiments/fidelity_estimation.py

@@ -185,7 +185,7 @@ def xeb_fidelity(
         ValueError: Circuit is inconsistent with qubit order or one of the
             bitstrings is inconsistent with the number of qubits.
     """
-    dim = np.product(circuit.qid_shape())
+    dim = np.prod(circuit.qid_shape(), dtype=np.int64)
 
     if isinstance(bitstrings, tuple):
         bitstrings = list(bitstrings)

cirq/experiments/fidelity_estimation_test.py

@@ -25,7 +25,7 @@ def sample_noisy_bitstrings(
     circuit: cirq.Circuit, qubit_order: Sequence[cirq.Qid], depolarization: float, repetitions: int
 ) -> np.ndarray:
     assert 0 <= depolarization <= 1
-    dim = np.product(circuit.qid_shape())
+    dim = np.prod(circuit.qid_shape(), dtype=np.int64)
     n_incoherent = int(depolarization * repetitions)
     n_coherent = repetitions - n_incoherent
     incoherent_samples = np.random.randint(dim, size=n_incoherent)

cirq/interop/quirk/cells/input_rotation_cells.py

@@ -162,7 +162,7 @@ def _apply_unitary_(self, args: 'cirq.ApplyUnitaryArgs'):
 
         target_axes = transposed_args.axes[: len(self.base_operation.qubits)]
         control_axes = transposed_args.axes[len(self.base_operation.qubits) :]
-        control_max = np.product([q.dimension for q in self.register]).item()
+        control_max = np.prod([q.dimension for q in self.register], dtype=np.int64).item()
 
         for i in range(control_max):
             operation = self.base_operation ** (self.exponent_sign * i / control_max)

cirq/ops/arithmetic_operation.py

@@ -192,12 +192,12 @@ def _apply_unitary_(self, args: 'cirq.ApplyUnitaryArgs'):
                 shape.append(1)
                 overflow_sizes.append(register + 1)
             else:
-                size = int(np.product([q.dimension for q in register]).item())
+                size = int(np.prod([q.dimension for q in register], dtype=np.int64).item())
                 shape.append(size)
                 input_ranges.append(range(size))
                 overflow_sizes.append(size)
 
-        leftover = args.target_tensor.size // np.product(shape).item()
+        leftover = args.target_tensor.size // np.prod(shape, dtype=np.int64).item()
         new_shape = (*shape, leftover)
 
         transposed_args = args.with_axes_transposed_to_start()

cirq/ops/controlled_operation.py

@@ -137,7 +137,7 @@ def _extend_matrix(self, sub_matrix: np.ndarray) -> np.ndarray:
         for control_vals in itertools.product(*self.control_values):
             active = (*(v for v in control_vals), *(slice(None),) * sub_n) * 2
             tensor[active] = sub_tensor
-        return tensor.reshape((np.prod(qid_shape, dtype=int).item(),) * 2)
+        return tensor.reshape((np.prod(qid_shape, dtype=np.int64).item(),) * 2)
 
     def _unitary_(self) -> Union[np.ndarray, NotImplementedType]:
         sub_matrix = protocols.unitary(self.sub_operation, None)

cirq/ops/fourier_transform.py

@@ -117,7 +117,7 @@ def _apply_unitary_(self, args: 'cirq.ApplyUnitaryArgs'):
         if isinstance(self.exponent, sympy.Basic):
             return NotImplemented
 
-        n = int(np.product([args.target_tensor.shape[k] for k in args.axes]))
+        n = int(np.prod([args.target_tensor.shape[k] for k in args.axes], dtype=np.int64))
         for i in range(n):
             p = 1j ** (4 * i / n * self.exponent)
             args.target_tensor[args.subspace_index(big_endian_bits_int=i)] *= p

cirq/ops/identity.py

@@ -77,7 +77,7 @@ def _has_unitary_(self) -> bool:
         return True
 
     def _unitary_(self) -> np.ndarray:
-        return np.identity(np.prod(self._qid_shape, dtype=int).item())
+        return np.identity(np.prod(self._qid_shape, dtype=np.int64).item())
 
     def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs') -> Optional[np.ndarray]:
         return args.target_tensor

cirq/ops/matrix_gates.py

@@ -61,7 +61,7 @@ def __init__(
         self._matrix = matrix
         self._qid_shape = tuple(qid_shape)
         self._name = name
-        m = int(np.prod(self._qid_shape))
+        m = int(np.prod(self._qid_shape, dtype=np.int64))
         if self._matrix.shape != (m, m):
             raise ValueError(
                 'Wrong matrix shape for qid_shape.\n'

cirq/ops/measurement_gate.py

@@ -130,7 +130,7 @@ def _measurement_key_(self):
         return self.key
 
     def _kraus_(self):
-        size = np.prod(self._qid_shape, dtype=int)
+        size = np.prod(self._qid_shape, dtype=np.int64)
 
         def delta(i):
             result = np.zeros((size, size))

cirq/ops/projector.py

@@ -68,7 +68,7 @@ def matrix(self, projector_qids: Optional[Iterable[raw_types.Qid]] = None) -> cs
             for qid in projector_qids
         ]
 
-        total_d = np.prod([qid.dimension for qid in projector_qids])
+        total_d = np.prod([qid.dimension for qid in projector_qids], dtype=np.int64)
 
         ones_idx = []
         for idx in itertools.product(*idx_to_keep):

cirq/ops/random_gate_channel.py

@@ -95,7 +95,9 @@ def _mixture_(self):
         if mixture is None:
             return None
 
-        do_nothing = np.eye(np.product(protocols.qid_shape(self.sub_gate)), dtype=np.float64)
+        do_nothing = np.eye(
+            np.prod(protocols.qid_shape(self.sub_gate), dtype=np.int64), dtype=np.float64
+        )
         result = [(p * float(self.probability), m) for p, m in mixture]
         result.append((1 - float(self.probability), do_nothing))
         return result
@@ -108,7 +110,9 @@ def _kraus_(self):
         if channel is None:
             return NotImplemented
 
-        do_nothing = np.eye(np.product(protocols.qid_shape(self.sub_gate)), dtype=np.float64)
+        do_nothing = np.eye(
+            np.prod(protocols.qid_shape(self.sub_gate), dtype=np.int64), dtype=np.float64
+        )
         result = [e * np.sqrt(self.probability) for e in channel]
         result.append(np.sqrt(1 - float(self.probability)) * do_nothing)
         return result

cirq/ops/random_gate_channel_test.py

@@ -160,7 +160,7 @@ def num_qubits(self) -> int:
 def assert_channel_sums_to_identity(val):
     m = cirq.kraus(val)
     s = sum(np.conj(e.T) @ e for e in m)
-    np.testing.assert_allclose(s, np.eye(np.product(cirq.qid_shape(val))), atol=1e-8)
+    np.testing.assert_allclose(s, np.eye(np.prod(cirq.qid_shape(val), dtype=np.int64)), atol=1e-8)
 
 
 def test_channel():

cirq/ops/raw_types.py

@@ -586,7 +586,7 @@ def _commutes_(
 
         # Don't create gigantic matrices.
         shape = protocols.qid_shape_protocol.qid_shape(circuit12)
-        if np.product(shape) > 2 ** 10:
+        if np.prod(shape, dtype=np.int64) > 2 ** 10:
             return NotImplemented  # coverage: ignore
 
         m12 = protocols.unitary_protocol.unitary(circuit12, default=None)

cirq/protocols/unitary_protocol.py

@@ -178,7 +178,7 @@ def _strat_unitary_from_apply_unitary(val: Any) -> Optional[np.ndarray]:
 
     if result is NotImplemented or result is None:
         return result
-    state_len = np.prod(val_qid_shape, dtype=int)
+    state_len = np.prod(val_qid_shape, dtype=np.int64)
     return result.reshape((state_len, state_len))
 
 
@@ -199,5 +199,5 @@ def _strat_unitary_from_decompose(val: Any) -> Optional[np.ndarray]:
     # Package result.
     if result is None:
         return None
-    state_len = np.prod(val_qid_shape, dtype=int)
+    state_len = np.prod(val_qid_shape, dtype=np.int64)
     return result.reshape((state_len, state_len))

cirq/qis/channels.py

@@ -21,7 +21,7 @@
 
 def kraus_to_choi(kraus_operators: Sequence[np.ndarray]) -> np.ndarray:
     """Returns the unique Choi matrix corresponding to a Kraus representation of a channel."""
-    d = np.prod(kraus_operators[0].shape)
+    d = np.prod(kraus_operators[0].shape, dtype=np.int64)
     c = np.zeros((d, d), dtype=np.complex128)
     for k in kraus_operators:
         v = np.reshape(k, d)

cirq/qis/measures.py

@@ -44,7 +44,7 @@ def _validate_int_state(state: int, qid_shape: Optional[Tuple[int, ...]]) -> Non
             f'but {state} was given.'
         )
     if qid_shape is not None:
-        dim = np.prod(qid_shape)
+        dim = np.prod(qid_shape, dtype=np.int64)
         if state >= dim:
             raise ValueError(
                 'Invalid state for given qid shape: '
@@ -155,10 +155,10 @@ def _numpy_arrays_to_state_vectors_or_density_matrices(
 ) -> Tuple[np.ndarray, np.ndarray]:
     if state1.ndim > 2 or (state1.ndim == 2 and state1.shape[0] != state1.shape[1]):
         # State tensor, convert to state vector
-        state1 = np.reshape(state1, (np.prod(state1.shape).item(),))
+        state1 = np.reshape(state1, (np.prod(state1.shape, dtype=np.int64).item(),))
     if state2.ndim > 2 or (state2.ndim == 2 and state2.shape[0] != state2.shape[1]):
         # State tensor, convert to state vector
-        state2 = np.reshape(state2, (np.prod(state2.shape).item(),))
+        state2 = np.reshape(state2, (np.prod(state2.shape, dtype=np.int64).item(),))
     if state1.ndim == 2 and state2.ndim == 2:
         # Must be square matrices
         if state1.shape == state2.shape:
@@ -171,32 +171,44 @@ def _numpy_arrays_to_state_vectors_or_density_matrices(
                 )
             if state1.shape == qid_shape:
                 # State tensors, convert to state vectors
-                state1 = np.reshape(state1, (np.prod(qid_shape).item(),))
-                state2 = np.reshape(state2, (np.prod(qid_shape).item(),))
+                state1 = np.reshape(state1, (np.prod(qid_shape, dtype=np.int64).item(),))
+                state2 = np.reshape(state2, (np.prod(qid_shape, dtype=np.int64).item(),))
         elif state1.shape[0] < state2.shape[0]:
             # state1 is state tensor and state2 is density matrix.
             # Convert state1 to state vector
-            state1 = np.reshape(state1, (np.prod(state1.shape).item(),))
+            state1 = np.reshape(state1, (np.prod(state1.shape, dtype=np.int64).item(),))
         else:  # state1.shape[0] > state2.shape[0]
             # state2 is state tensor and state1 is density matrix.
             # Convert state2 to state vector
-            state2 = np.reshape(state2, (np.prod(state2.shape).item(),))
-    elif state1.ndim == 2 and state2.ndim < 2 and np.prod(state1.shape) == np.prod(state2.shape):
+            state2 = np.reshape(state2, (np.prod(state2.shape, dtype=np.int64).item(),))
+    elif (
+        state1.ndim == 2
+        and state2.ndim < 2
+        and np.prod(state1.shape, dtype=np.int64) == np.prod(state2.shape, dtype=np.int64)
+    ):
         # state1 is state tensor, convert to state vector
-        state1 = np.reshape(state1, (np.prod(state1.shape).item(),))
-    elif state1.ndim < 2 and state2.ndim == 2 and np.prod(state1.shape) == np.prod(state2.shape):
+        state1 = np.reshape(state1, (np.prod(state1.shape, dtype=np.int64).item(),))
+    elif (
+        state1.ndim < 2
+        and state2.ndim == 2
+        and np.prod(state1.shape, dtype=np.int64) == np.prod(state2.shape, dtype=np.int64)
+    ):
         # state2 is state tensor, convert to state vector
-        state2 = np.reshape(state2, (np.prod(state2.shape).item(),))
+        state2 = np.reshape(state2, (np.prod(state2.shape, dtype=np.int64).item(),))
 
     if validate:
-        dim1: int = state1.shape[0] if state1.ndim == 2 else np.prod(state1.shape).item()
-        dim2: int = state2.shape[0] if state2.ndim == 2 else np.prod(state2.shape).item()
+        dim1: int = (
+            state1.shape[0] if state1.ndim == 2 else np.prod(state1.shape, dtype=np.int64).item()
+        )
+        dim2: int = (
+            state2.shape[0] if state2.ndim == 2 else np.prod(state2.shape, dtype=np.int64).item()
+        )
         if dim1 != dim2:
             raise ValueError('Mismatched dimensions in given states: ' f'{dim1} and {dim2}.')
         if qid_shape is None:
             qid_shape = (dim1,)
         else:
-            expected_dim = np.prod(qid_shape)
+            expected_dim = np.prod(qid_shape, dtype=np.int64)
             if dim1 != expected_dim:
                 raise ValueError(
                     'Invalid state dimension for given qid shape: '

cirq/qis/states.py

@@ -96,7 +96,7 @@ def __init__(
             qid_shape = infer_qid_shape(data)
         self._data = data
         self._qid_shape = qid_shape
-        self._dim = np.prod(self.qid_shape, dtype=int).item()
+        self._dim = np.prod(self.qid_shape, dtype=np.int64).item()
         if validate:
             self.validate(dtype=dtype, atol=atol)
 
@@ -263,7 +263,7 @@ def quantum_state(
                 'Please specify the qid shape explicitly using '
                 'the qid_shape argument.'
             )
-        dim = np.prod(qid_shape, dtype=int).item()
+        dim = np.prod(qid_shape, dtype=np.int64).item()
         if not 0 <= state < dim:
             raise ValueError(
                 f'Computational basis state is out of range.\n'
@@ -281,7 +281,7 @@ def quantum_state(
         if qid_shape is None:
             qid_shape = infer_qid_shape(state)
         if data.ndim == 1 and data.dtype.kind != 'c':
-            if len(qid_shape) == np.prod(qid_shape, dtype=int):
+            if len(qid_shape) == np.prod(qid_shape, dtype=np.int64):
                 raise ValueError(
                     'Because len(qid_shape) == product(qid_shape), it is '
                     'ambiguous whether the given state contains '
@@ -412,7 +412,7 @@ def infer_qid_shape(*states: 'cirq.QUANTUM_STATE_LIKE') -> Tuple[int, ...]:
         )
 
     # check if the shape is compatible with the states specified as integers
-    if integer_states and np.prod(qid_shape, dtype=int) <= max(integer_states):
+    if integer_states and np.prod(qid_shape, dtype=np.int64) <= max(integer_states):
         raise ValueError(
             'Failed to infer the qid shape of the given states. '
             f'The given integer state {max(integer_states)} is too high for the '
@@ -566,7 +566,7 @@ def _intersection_explicit_with_unfactorized_qid_shapes(
     return {
         qid_shape
         for qid_shape in explicit_qid_shapes
-        if np.prod(qid_shape, dtype=int) == unfactorized_total_dimension
+        if np.prod(qid_shape, dtype=np.int64) == unfactorized_total_dimension
     }
 
 
@@ -698,7 +698,7 @@ def density_matrix_from_state_vector(
         sum_inds.tolist(),
         indices + sum_inds[indices].tolist(),
     )
-    new_shape = np.prod([shape[i] for i in indices], dtype=int)
+    new_shape = np.prod([shape[i] for i in indices], dtype=np.int64)
 
     return rho.reshape((new_shape, new_shape))
 
@@ -887,10 +887,10 @@ def validate_normalized_state_vector(
                 dtype, state_vector.dtype
             )
         )
-    if state_vector.size != np.prod(qid_shape, dtype=int):
+    if state_vector.size != np.prod(qid_shape, dtype=np.int64):
         raise ValueError(
             'state_vector has incorrect size. Expected {} but was {}.'.format(
-                np.prod(qid_shape, dtype=int), state_vector.size
+                np.prod(qid_shape, dtype=np.int64), state_vector.size
             )
         )
     norm = np.sum(np.abs(state_vector) ** 2)
@@ -914,7 +914,7 @@ def validate_qid_shape(
     size = state_vector.size
     if qid_shape is None:
         qid_shape = (2,) * (size.bit_length() - 1)
-    if size != np.prod(qid_shape, dtype=int):
+    if size != np.prod(qid_shape, dtype=np.int64):
         raise ValueError(
             'state_vector.size ({}) is not a power of two or is not a product '
             'of the qid shape {!r}.'.format(size, qid_shape)
@@ -1006,7 +1006,7 @@ def validate_density_matrix(
             f'Incorrect dtype for density matrix: Expected {dtype} '
             f'but has dtype {density_matrix.dtype}.'
         )
-    expected_shape = (np.prod(qid_shape, dtype=int),) * 2
+    expected_shape = (np.prod(qid_shape, dtype=np.int64),) * 2
     if density_matrix.shape != expected_shape:
         raise ValueError(
             f'Incorrect shape for density matrix: Expected {expected_shape} '
@@ -1083,6 +1083,6 @@ def eye_tensor(half_shape: Tuple[int, ...], *, dtype: 'DTypeLike') -> np.ndarray
     Returns:
         The created numpy array with shape `half_shape + half_shape`.
     """
-    identity = np.eye(np.prod(half_shape, dtype=int).item(), dtype=dtype)
+    identity = np.eye(np.prod(half_shape, dtype=np.int64).item(), dtype=dtype)
     identity.shape = half_shape * 2
     return identity

cirq/sim/density_matrix_simulator.py

@@ -352,7 +352,7 @@ def density_matrix(self, copy=True):
             state = self._merged_sim_state
             if state is not None:
                 matrix = state.target_tensor
-                size = int(np.sqrt(np.prod(matrix.shape, dtype=int)))
+                size = int(np.sqrt(np.prod(matrix.shape, dtype=np.int64)))
                 self._density_matrix = np.reshape(matrix, (size, size))
         return self._density_matrix.copy() if copy else self._density_matrix
 
@@ -435,7 +435,7 @@ def __init__(
         super().__init__(
             params=params, measurements=measurements, final_simulator_state=final_simulator_state
         )
-        size = np.prod(protocols.qid_shape(self), dtype=int)
+        size = np.prod(protocols.qid_shape(self), dtype=np.int64)
         self.final_density_matrix = np.reshape(
             final_simulator_state.density_matrix.copy(), (size, size)
         )