BTW, on the principle of "where else did we make the same mistake",
I looked through the other aggregates using float8_regr_accum.
Most seem okay, but float8_regr_intercept does this:
PG_RETURN_FLOAT8((Sy - Sx * Sxy / Sxx) / N);
Seems to me that expression is also prone to internal
overflow/underflow. Underflow probably isn't a huge issue,
since the result will reduce to Sy/N which is likely to be good
enough. But can we do anything about overflow?
One simple change that might make things better is to compute
PG_RETURN_FLOAT8((Sy - Sx * (Sxy / Sxx)) / N);
on the theory that the sums of products are likely to both be large.
regards, tom lane