I believe this code does what you want:

#include <iostream>
#include <cmath>

float shortest_distance_mod(double A, double B, double M) {
  return std::abs(std::fmod(std::abs(A - B) + M * .5, M) - M * .5);
}

int main() {
  std::cout << shortest_distance_mod(3.0, 4.0, 4.0) << '\n';
}

I should mention that there are two other of U(1) that deserve consideration: 2x2 special orthogonal matrices and unit-length complex numbers. The latter is particularly useful. You can think of it as representing rotations using the vector (cos(alpha), sin(alpha)). Making them complex numbers makes sense because complex multiplication then corresponds to composition of rotations (i.e., sum of angles).

Through the years I have posted several snippets of code using this representation instead of angles. I'll be happy to produce code for whatever operations you currently do with angles, if you are interested. The code is typically more concise, it is easier to get right and it requires fewer trigonometric operations than using angles.

Thanks for the code snippet. Though I must say I don't understand it very well. I see that A is current location, B is the desired, M is probably the length of the circle and abs(A - B) is the distance, but I don't see why we have to add + M * .5 and that two times (second time with subtraction), use mod() and than again abs()?

How I originally thought to solve the problem is to first transform it to radian space. I thought to multiply my domain x[0, 4] like x * 2pi. This would get me to [0, 25.12] which is I think not what I wanted. I think I should somehow convert [0, 4] to [0, 6.28] if I'm correct and than I'll be able to use sin(), cos() to get the distances and other things?

I feel like I'm so close to understand this but somehow I'm not quite there. Could you please explain this a bit further, why do we use modulo here (@alvaro solution) for the expected value? How does that returns us to the desired range?

Edit: Oh, I think I see, we just scratch the extra value, the value that got us to 8pi (4 times around the circle).
Now it looks like to me that it's mandatory to use modulo after a multiplication so we could get back to transformed space, is this true?

Give A-B a name, like x. Then plot these functions:

x

abs(x)

abs(x) + 2

fmod(abs(x) + 2, 4)

fmod(abs(x) + 2, 4) - 2

abs(fmod(abs(x) + 2, 4) - 2

You'll see that the last one has the shape you need: https://www.wolframalpha.com/input/?i=plot+of+(abs(fmod(abs(x)+%2B+2,+4)+-+2)+for+x+between+-8+and+8

Perhaps that sequence of transformations will make sense then.

Ok, I plotted every function. I can see where you were going with this through I don't understand it yet completely. Maybe I'll get to it later.
@swiftcoder mentioned that you are using mod() to wrap the function back to desired range but you were not converting it to radians in the first place so I'm left a bit confused. I see that you had to bring it back because of +2 (+M*.5) but it's still not entirely clear to me.

It came to my mind that when I tried to convert [0, 4] to [0, 6.28] radians I can do it with percentages too. If we say that the first part is in domain K and second in domain L, then if we have some angle A_k in domain K, we can divide it by K Lenght, that is 4. This will give us the percentage p. If we multiply L Lenght * p it will give us the same angle in domain L, that is a angle A_l in radians.
Now, I know that the next calculation is wrong, but I guess that we could do something like we did before to get the same answer, try to multiply A_k * 6.28 and mod() this with K/L Lenght. This doesn't give us A_l but I'm sure that we could do something similar to arrive to it?

I don't understand where your confusion is coming from. Why are you so bent on using radians? In the question you posed about computing the smallest distance between 3 and 4 when wrapping modulo 4, radians are nowhere to be found; not in the question, and not in the answer. You just need a saw-tooth pattern, and the formula I gave you produces it.