Repository: yt-project/unyt

Issue:
Segmentation fault with `u.Msun`
* unyt version: https://github.com/kyleaoman/unyt/tree/hist_bugfix commit 13b6da058e5387443adc2b3d8a52baaf4cff4eaf
* Python version: 3.12
* Operating System: Linux (Fedora 40)

### Description

Trying to histogram a unyt array initialized with `unyt_array(arr, units)` in density mode when units are `u.Msun` crashes with a segmentation fault... but not when units are `u.solMass`, and not when the array is initialized as `np.array(...) * u.Msun`.

### What I Did

This is rearing its head in the updated handling of the `histogramdd` argument `density` but is pretty clearly more fundamental than anything I've directly introduced. Here's as much as I've managed to isolate it so far:

```
import numpy as np
import unyt as u

sample = [
    u.unyt_array(np.arange(3), u.Msun),
]
np.histogramdd(
    sample,
    density=True,
)
```
That exits with a segmentation fault.

This however exits cleanly:
```
import numpy as np
import unyt as u

sample = [
    u.unyt_array(np.arange(3), u.solMass),
]
np.histogramdd(
    sample,
    density=True,
)
```

So does this:
```
import numpy as np
import unyt as u

sample = [
    np.arange(3) * u.Msun,
]
np.histogramdd(
    sample,
    density=True,
)
```

Discussion:
@neutrinoceros: Can you share the exact version of Python you're using (down to the patch level), your provider (Python.org, conda, something else ?) and the version of numpy you are running with ? unyt itself is pure Python so I highly doubt that the memory access error is on our side.@neutrinoceros: FWIW, running your first example on my system produces a `RecursionError`, not a segfault. This is very likely showing a problem in unyt itself.@neutrinoceros: Note that the key difference betwee `unyt.msun` and `unyt.solMass` is that the former is a quantity (whose value is stretching the limits of floating point arithmetic) while the latter is a unit.@neutrinoceros: I'm able to avoid this `RecursionError` with the following patch
```patch
diff --git a/unyt/_array_functions.py b/unyt/_array_functions.py
index 62aca69..a752fd8 100644
--- a/unyt/_array_functions.py
+++ b/unyt/_array_functions.py
@@ -288,7 +288,8 @@ def _histogramdd(
     if density:
         for s in sample:
             counts /= getattr(s, "units", 1)
-    counts *= getattr(weights, "units", 1)
+    if weights is not None and hasattr(weights, "units"):
+        counts *= weights.units
     return counts, tuple(_bin * getattr(s, "units", 1) for _bin, s in zip(bins, sample))
```
(it also seems to work when I pass `weights`, regardless if it has units or not). I'm going to propose patching this case and similar ones from your previous PR.@kyleaoman: That seems like it should be equivalent, wonder what's triggering the recursion... also wonder if it's a deeper issue that might crop up elsewhere, or if this fix is enough.

(Also just to mention - when I first encountered this in a pytest run it came with a partial stack trace before the segfault, that indeed looked like something recursing a lot.)@neutrinoceros: > That seems like it should be equivalent, wonder what's triggering the recursion... 

multipliying an array by a scalar `1` calls in some ufunc/array func logic, but at this point we're already deep in `unyt_array.__array_function__`. I think that's where the issue is coming from, though I must admit I didn't look much closer after I found a working solution. PR incoming.

Solution patch:
diff --git a/unyt/_array_functions.py b/unyt/_array_functions.py
index 62aca69b..3bd811b7 100644
--- a/unyt/_array_functions.py
+++ b/unyt/_array_functions.py
@@ -286,9 +286,16 @@ def _histogramdd(
     # don't getattr(..., "units", NULL_UNIT) because e.g. we don't want
     # a unyt_array if weights are not a unyt_array and not density
     if density:
+        divider_units = 1 * NULL_UNIT
         for s in sample:
-            counts /= getattr(s, "units", 1)
-    counts *= getattr(weights, "units", 1)
+            if not hasattr(s, "units"):
+                continue
+            divider_units *= s.units
+        if divider_units != NULL_UNIT:
+            counts /= divider_units
+
+    if weights is not None and hasattr(weights, "units"):
+        counts *= weights.units
     return counts, tuple(_bin * getattr(s, "units", 1) for _bin, s in zip(bins, sample))
 
 


Test patch:
diff --git a/unyt/tests/test_array_functions.py b/unyt/tests/test_array_functions.py
index a039939e..291db7d8 100644
--- a/unyt/tests/test_array_functions.py
+++ b/unyt/tests/test_array_functions.py
@@ -9,7 +9,7 @@
 from numpy.testing import assert_allclose
 from packaging.version import Version
 
-from unyt import A, K, cm, degC, delta_degC, g, km, rad, s
+from unyt import A, K, Msun, cm, degC, delta_degC, g, km, rad, s
 from unyt._array_functions import (
     _HANDLED_FUNCTIONS as HANDLED_FUNCTIONS,
     _UNSUPPORTED_FUNCTIONS as UNSUPPORTED_FUNCTIONS,
@@ -831,6 +831,12 @@ def test_histogramdd_with_weights_and_dimless_arr(self):
         assert not hasattr(ywbins2, "units")
         assert not hasattr(zwbins2, "units")
 
+    @pytest.mark.parametrize("weights", [None, [0, 1, 2], [0, 1, 2] * cm])
+    def test_histogramdd_recursion(self, weights):
+        # regression test for https://github.com/yt-project/unyt/issues/540
+        sample = [unyt_array(np.arange(3), Msun)]
+        np.histogramdd(sample, density=True, weights=weights)
+
     def test_histogram_bin_edges(self):
         rng = np.random.default_rng()
         arr = rng.normal(size=1000) * cm