task: metadata: install: pip install -e .
patch: diff --git a/unyt/_array_functions.py b/unyt/_array_functions.py
index 62aca69..3bd811b 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))
 
 

problem_statement: 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,
)
```

repo: yt-project/unyt
test_cmd: pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
test_patch: diff --git a/unyt/tests/test_array_functions.py b/unyt/tests/test_array_functions.py
index a039939..291db7d 100644
--- a/unyt/tests/test_array_functions.py
+++ b/unyt/tests/test_array_functions.py
@@ -9,7 +9,7 @@ import pytest
 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 @@ class TestHistograms:
         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