erf(1b-20) -> 1.387778780781446b-17
The correct answer is closer to 1.128379167095513b-20.
This is caused by bfloat-erf computing erf using 1-gamma_incomplete(1/2,x^2)/sqrt(%pi). For small x, gamma_incomplete(1/2,x^2) is very close to sqrt(%pi), so we lose lots of precision using this formula.
For small x, we should just use the taylor series.