If anyone ever needs to compute the integral of softplus (2-fold integral of sigmoid), it's -Li_2(-exp(x)) where Li_2 is the dilogarithm, and you can find a very good and fast numerical approximation to the dilogarithm of real argument in https://arxiv.org/abs/2201.01678.
Comparison of methods for the calculation of the real dilogarithm regarding instruction-level parallelism

We compare different methods for the computation of the real dilogarithm regarding their ability for using instruction-level parallelism when executed on appropriate CPUs. As a result we present an instruction-level-aware method and compare it to existing implementations.

arXiv.org