Wrong sign in numerical integration, possible precision issue

I need to integrate the following function:

where z > 0. The problem is that the integrand is very small for large z and high precision is required in the integration. So far, I have written the integrand as

double integrand__W(double x, double z){
    double arg = z*z/(4.0*x);
    
    double num = exp(arg+x)+1;
    double den1 = expm1(arg);
    double den2 = exp(x);

    num = isinf(num) ? arg+x : log(num);
    den1 = isinf(den1) ? arg : log(den1);
    den2 = x; //log(exp(x))=x
    double t1 = num-den1-den2;
    
    num = exp(x);
    double den = exp(x)+1;
    double t2 = isinf(den) ? exp(-x) : num/(den*den);
    
    return t1*t2;
}

For numerical integration, I'm using Cubature, a simple C-package for adaptive multidimensional integration:

//integrator
struct fparams {
    double z;
};

int inf_W(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval){
    struct fparams * fp = (struct fparams *)fdata;
    double z = fp->z;
    double t = x[0]; 
    double aux = integrand__W(a_int+t*pow(1.0-t, -1.0), z)*pow(1.0-t, -2.0);
    if (!isnan(aux) && !isinf(aux))
    {
        fval[0] = aux;
    }
    else
    {
        fval[0] = 0.0;
    }
    return 0;
}

//range integration 1D
    size_t maxEval = 1e7;
    double xl[1] = { 0 };
    double xu[1] = { 1 };

    double W, W_ERR;
    struct fparams params = {z};
    hcubature(1, inf_W, &params, 1, xl, xu, maxEval, 0, 1e-5, ERROR_INDIVIDUAL, &W, &W_ERR);
    cout << "z: " << z << " | " << W << " , " << W_ERR << endl;

where the integration over the semi-infinite interval is possible by a change of variables:

Analytically, I know that the integrated is non-negative, so the integral itself should be non-negative. However, I'm getting some incorrect results due to a lack of accuracy:

z: 100 | -3.97632e-17 , 1.24182e-16

In Mathematica, working with high precision, I can get the desired result:

w[x_, z_] := E^x/(E^x + 1)^2 Log[(E^(z^2/(4 x)) + E^-x)/(E^(z^2/(4 x)) - 1)]

W[z_?NumericQ] := NIntegrate[w[x, z], {x, 0, ∞},
  WorkingPrecision -> 40,
  Method -> "LocalAdaptive"]

W[100]

(* 4.679853458969239635780655689865016458810*10^-43 *)

My question: Is there any way to write my integrand such that I can reach the required precision? Thanks.

After asking the same question to a different community, I got two suggestions that seem to work:

Avoiding subtractive cancellation

Manipulate the integral a little bit first:

and then rewrite the integrand as

double integrand__W(double x, double z){
    double arg = z*z/(4.0*x);
    
    double t1 = log1p((exp(-x)+1)/expm1(arg));
    
    double num = exp(x);
    double den = exp(x)+1;
    double t2 = isinf(den) ? exp(-x) : num/(den*den);
    
    return t1*t2;
}

Use of Exp-Sinh quadrature

This integration scheme is provided by the Boost library:

#include <iostream>
#include <cmath>
#include <boost/math/quadrature/exp_sinh.hpp>

using boost::math::quadrature::exp_sinh;
using std::exp;
using std::expm1;
using std::log;

int main() {
    exp_sinh<double> integrator;
    double z = 100.0;
    auto f = [z](double x) {
        double k1 = 1.0/(2 + exp(-x) +exp(x));
        double t = z*z/(4*x);
        double log_arg;
        if (t > 1) {
            log_arg = (1 + exp(-x)*exp(-t))/(1 - exp(-t));
        } else {
            log_arg = (exp(t) + exp(-x))/expm1(t);
        }
        return k1*log(log_arg);
    };
    double termination = sqrt(std::numeric_limits<double>::epsilon());
    double error;
    double L1;
    double Q = integrator.integrate(f, termination, &error, &L1);
    std::cout << "Q = " << Q << ", error estimate: " << error << "\n";
}