Ray vs AABB algorithm that also gives which side was hit? There are plenty of well-known algorithms for determining if a ray hit an Axis-Aligned Bounding Box (AABB), like Andrew Woo's covered in a response here. 

Using exit = origin + (dir * tMax) will give you the exact position of exit using distance through the ray. (and same for tMin and the entry)

However, I need a way to not only do that, but also get if we hit the left, right, top, bottom, front, or rear face of the AABB.

The naive approach would be to do some kind of floating-point comparison given the minpoint/maxpoint and try to see what side (if any) the point exactly equals in each dimension, or exactly equals within a tolerance. That seems error-prone, I was wondering if there was a better way.

Is there an algorithm for this?
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.