Ray vs AABB algorithm that also gives which side was hit?
Given an intersection algorithm like the one you linked which first generates a $t$ value for every plane of the AABB, then computes $t_{min}$ and $t_{max}$ from those, you simply need to figure out which of the six initial values correspond to your final $t_{min}$. There will be two matches if your ray hits an edge, three if it hits a corner. In that case, this may include a back face (see the fix below).

This will not require a tolerance as $t_{min}$ is always one of the six values. In fact, it is always one of three intermediate values and you can use that to make your check simpler, faster and better. Given code, from the link above, such as:

float tmin = max(max(min(t1, t2), min(t3, t4)), min(t5, t6));


you could rewrite as:

float tminx = min(t1, t2);
float tminy = min(t3, t4);
float tmin = max(max(tminx, tminy), min(t5, t6));
int axis = 2;
if (tmin == tminx) axis = 0;
if (tmin == tminy) axis = 1;


As an added bonus, unlike testing all six values, this will not match the back face in the event that a far edge is hit.