Hi team, I recently became aware of gh-18754 <https://github.com/scipy/scipy/pull/18754>, a PR to make the factorial function in special return a 0d array when given 0d input. This brings it in line with the behavior of factorial2, and factorialk, which was settled down in gh-15600 <https://github.com/scipy/scipy/issues/15600>. The longstanding behavior for ufuncs in NumPy and SciPy has been to return a scalar when the output should otherwise be expected to be a 0d-array. All ufuncs generated in scipy.special through the standard infrastructure will have the expected behavior. The factorial functions are however not implemented with the standard infrastructure in order to allow exact calculation with arbitrary-precision integers. I would prefer to have consistency between the factorial functions and the other ufuncs in special. There however isn't a consensus on what the correct behavior should be. I would appreciate any input anyone has to offer, especially from those familiar with the reasons behind the current default convention to convert 0d-arrays to scalar before returning an array from the NumPy C API into Python. Feel free to join the discussion at gh-18768 <https://github.com/scipy/scipy/issues/18768>. Thanks everyone! Albert