One thought: would the deci_sqrt approach help with value ranges where the values are well within float limits, but the squares of the values are not? E.g., on my machine, I currently get errors for both of the following:

>>> xs = [random.normalvariate(0.0, 1e200) for _ in range(10**6)]
>>> statistics.stdev(xs)

>>> xs = [random.normalvariate(0.0, 1e-200) for _ in range(10**6)]
>>> statistics.stdev(xs)

It's hard to imagine that there are too many use-cases for values of this size, but it still feels a bit odd to be constrained to only half of the dynamic range of float.