In such a world there are infinite number of paths from S to T. Let's denote the coordinates of T by (Tx, Ty), the coordinates of S by (Sx, Sy), and the size of the world by (Wx, Wy). The wrapped coordinates of T are (Tx + j * Wx, Ty + i * Wy), where i and j are integers, that is, elements of the set {..., -2, -1, 0, 1, 2, ...}. The vectors connecting S to T are (Dx, Dy) := (Tx + j * Wx - Sx, Ty + i * Wy - Sy). For a given (i, j) pair, the distance is the length of the vector, sqrt(Dx * Dx + Dy * Dy), and the direction in radians is atan(Dy / Dx). The shortest path is one of the 9 paths, where i and j are in {-1, 0, 1}:

I would calculate all of these 9 distances, and find the shortest path using simple linear search. (Calculating only the square of distances is enough, so there is no need to take the square root.)

In such a world there are infinite number of paths from S to T. Let's denote the coordinates of T by (Tx, Ty), the coordinates of S by (Sx, Sy), and the size of the world by (Wx, Wy). The wrapped coordinates of T are (Tx + i * Wx, Ty + j * Wy), where i and j are integers, that is, elements of the set {..., -2, -1, 0, 1, 2, ...}. The vectors connecting S to T are (Dx, Dy) := (Tx + i * Wx - Sx, Ty + j * Wy - Sy). For a given (i, j) pair, the distance is the length of the vector, sqrt(Dx * Dx + Dy * Dy), and the direction in radians is atan(Dy / Dx). The shortest path is one of the 9 paths, where i and j are in {-1, 0, 1}:

The i and j values for the shortest path can be determined directly:

int i = Sx - Tx > Wx / 2 ? 1 : Sx - Tx < -Wx / 2 ? -1 : 0;
int j = Sy - Ty > Wy / 2 ? 1 : Sy - Ty < -Wy / 2 ? -1 : 0;

Thank you, @IlmariKaronen, @SamHocevar and @romkyns for your help!