Computes the extended GCD and Bezout coefficients for two numbers.
Uses the Extended Euclidean algorithm to find (g, s, t, sub_direction) where g = gcd(a, b).
The relationship between inputs and outputs is:
If sub_direction is true: g = s * a - t * b
If sub_direction is false: g = t * b - s * a
Returns a tuple (g, s, t, sub_direction) where g is the GCD and (s, -t) or (-s, t) are the
Bezout coefficients (according to sub_direction).