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.