It took me a good few hours but I've finally figured it out. The idea is that you figure out which point on the bounding box is closest to the ray, and then compute distance from that point to the ray. It can also compute barycentric coordinate on the ray (t), but if it intersects the bounding box, it will be in some arbitrary point inside of it.
float rayAABBdist ( vec2 p, vec2 v, float tt, float ll, float bb, float rr )
{
float tx1 = ( ll - p.x ) / v.x;
float tx2 = ( rr - p.x ) / v.x;
float ty1 = ( tt - p.y ) / v.y;
float ty2 = ( bb - p.y ) / v.y;
float p1 = max ( 0.0, max ( min ( tx1, tx2 ), min ( ty1, ty2 ) ) );
float p2 = max ( 0.0, min ( max ( tx1, tx2 ), max ( ty1, ty2 ) ) );
float x = max ( min ( ( p.x + v.x * p1 + p.x + v.x * p2 ) / 2, rr ), ll );
float y = max ( min ( ( p.y + v.y * p1 + p.y + v.y * p2 ) / 2, bb ), tt );
float t = max ( 0.0, ( ( x - p.x ) * v.x + ( y - p.y ) * v.y ) / ( v.x * v.x + v.y * v.y ) );
x = p.x + v.x * t - x;
y = p.y + v.y * t - y;
return x * x + y * y;
}
A variant that handles a line segment is the same except it also clamps coefficients to 1.0 on the top. It's not the most optimal solution, but it works. It is pretty expensive though, so depending on your application you might just use plain raycast instead. I've tried some different approaches but they all fail in edge cases.
Two years looks like pretty late but if anyone's gonna web search for this problem, they will find this solution.