Autocorrelation via real FFT and real IFFT

I am working on a fast autocorrelation function using the FFT (similar to scipy.signal.fftconvolve). The algorithm is pretty straight forward but I feel like I'm leaving an optimization on the table and I can't figure out how to make it work. Here is an outline of the (usual?) algorithm:

1. Run a real signal $$x$$ through the discrete Fourier transform: $$X = \mathfrak{F}\{x\}$$.
2. Replace the entries in $$X$$ with their norm: $$Y = X \cdot X^*$$. (Don't care about DC for now.)
3. Run $$Y$$ through the inverse discrete Fourier transform the get the autocorrelated real signal: $$y = \mathfrak{F}^{-1}\{Y\}$$

With the forward FFT I can make use of the fact that I have a real signal, so the computation is faster than for complex signals of the same length. So far so good. Here is where I'm stuck: $$Y$$ is also a real signal, since the product of a complex number with its conjugate is always real. Right now I am using a backward real FFT function that expects a complex signal and returns a real signal. Is there a trick similar to the forward real FFT? Something like a Decimation In Time (DIT) algorithm for the input of the IFFT.

Here is an example in Python:

import numpy as np
x = [2, 1, 3, 1, 3, 1, 2, 1]
X = np.fft.rfft(x)
Y = X * np.conj(X)
y = np.fft.irfft(Y)  # <- feel like this could be improved, Y is a real signal.
print(x, X, Y, y, sep='\n')

# x: [2 1 3 1 3 1 2 1]
# X: [14.+0.j -1.-1.j  0.-0.j -1.+1.j  6.+0.j]
# Y: [196.+0.j   2.+0.j   0.+0.j   2.+0.j  36.+0.j]
# y: [30. 20. 29. 20. 28. 20. 29. 20.]


There is such a function, though numpy might not expose it: it's still a real-valued-input FFT! So, it's, aside from an exponent's sign, the same as your rfft. Matter of fact, FFT libraries like FFTw do expose a backwards r2c FFT variant.
And fun fact: your input to the rfft is real, so step 1. is a NOP (no-operation). Step 2 is your rfft, and step 3 is still relatively trivial to compute.