It’s already implemented here: GitHub - JuliaStats/LogExpFunctions.jl: Julia package for various special functions based on `log` and `exp`.
I assume you know this, but your expression is analytically equivalent to \log\left(\sum_i e^{x_i}\right), i.e., just factoring out the largest term. You could also write it as y = x_\mathrm{max} + \log\left(\sum_i e^{x_i-x_\mathrm{max}}\right) and not worry about any index shenanigans.