block_replicate promotes float32/float16 input to float64 when conserve_sum=True #20360

Description

[this was mostly written by Claude at my direction]

With the default conserve_sum=True, astropy.nddata.block_replicate returns a float64 array for float32 or float16 input (and complex128 for complex64), while float64 input round-trips and conserve_sum=False leaves every dtype unchanged. Since FITS image data is overwhelmingly float32, block_replicate(image, n) silently doubles the in-memory size of the result for the common case.

The promotion traces to astropy/nddata/blocks.py, block_replicate:

data, block_size = _process_block_inputs(data, block_size)
for i in range(data.ndim):
    data = np.repeat(data, block_size[i], axis=i)

if conserve_sum:
    # in-place division can fail due to dtype casting rule
    data = data / np.prod(block_size)

return data

_process_block_inputs turns block_size into an ndarray via np.atleast_1d(block_size).astype(int), so np.prod(block_size) is an np.int64 scalar (on a platform with 64-bit default integers). Under NEP 50, dividing a float32 array by an int64 scalar promotes the result to float64 (and float16 and complex64 likewise promote to their next-wider type). This only happens in the conserve_sum=True (default) branch, since conserve_sum=False returns the repeated array with no further arithmetic, which is why the dtype is preserved there.

The behaviour is the same on main at 96774df (2026-09-10), where block_replicate is identical to the 8.0.1 release.

Expected behavior

The output dtype matches the input dtype regardless of conserve_sum, as it already does for float64 input and for conserve_sum=False.

How to Reproduce

>>> import numpy as np
>>> from astropy.nddata import block_replicate
>>> block_replicate(np.ones((2, 2), dtype=np.float32), 2).dtype
dtype('float64')
>>> block_replicate(np.ones((2, 2), dtype=np.float32), 2, conserve_sum=False).dtype
dtype('float32')
>>> block_replicate(np.ones((2, 2), dtype=np.float16), 2).dtype
dtype('float64')
>>> block_replicate(np.ones((2, 2), dtype=np.float64), 2).dtype
dtype('float64')
>>> block_replicate(np.ones((2, 2), dtype=np.complex64), 2).dtype
dtype('complex128')

Versions

import astropy
astropy.system_info()

platform
--------
platform.platform() = 'macOS-26.6.2-arm64-arm-64bit'
platform.version() = 'Darwin Kernel Version 25.6.0: Fri Jul 31 19:18:49 PDT 2026; root:xnu-12377.161.14~5/RELEASE_ARM64_T6000'
platform.python_version() = '3.12.13'

packages
--------
astropy              8.0.1
numpy                2.5.2
scipy                1.18.1
matplotlib           --
pandas               --
pyerfa               2.0.1.5