APL(NARS), 86 chars

{(g x y)←{(u v)←⍵⋄v=0:u,1,0⋄(g x y)←∇v,v∣u⋄g,y,x-y×⌊u÷v}⍵⋄x≥0:x y⋄(a b)←⍵÷g⋄(b+x),y-a}

Calculate egcd (extended gcd) with a recursive function, than aggiust the result in the way the first output number is >= 0. I don't know if it is right.

test:

  f←{(g x y)←{(u v)←⍵⋄v=0:u,1,0⋄(g x y)←∇v,v∣u⋄g,y,x-y×⌊u÷v}⍵⋄x≥0:x y⋄(a b)←⍵÷g⋄(b+x),y-a}
  f 42 12
1 ¯3 
  f 4096 84
4 ¯195 
  f 5 3
2 ¯3 
  f 1155 405
20 ¯57 
  f 23780 177
20 ¯2687 
  f 11235813  112358
8643 ¯864301 

Janet, 126 bytes

|((fn e[b a x u y v](case b 0[(mod u(/ $1 a))(mod v(/(- $)a))](e(% a b)b(- u(*(div a b)x))x(- v(*(div a b)y))y)))$1 $ 0 1 1 0)

Uses the extended Euclidean algorithm. The trick to make the first coefficient positive is taken from @Sherlock9’s Python answer.

Python 3, 101 106 bytes

Edit: Added some improvements and corrections suggested by Bruce_Forte.

An answer that uses the extended Euclidean algorithm. It's a bit clunky in places though, and I hope to golf it some more. I could convert to Python 2 to save a byte on integer division (//) but I'm not sure how Python 2's % modulus operator works with a negative second argument, as that is crucial for getting the output right.

def e(a,b):
 r=b;x=a;s=z=0;t=y=1
 while r:q=x/r;x,r=r,x%r;y,s=s,y-q*s;z,t=t,z-q*t
 return y%(b/x),z%(-a/x)

Ungolfed:

def e(a, b):
    r = b
    x = a    # becomes gcd(a, b)
    s = 0
    y = 1    # the coefficient of a
    t = 1
    z = 0    # the coefficient of b
    while r:
        q = x / r
        x, r = r, x % r
        y, s = s, y - q * s
        z, t = t, z - q * t
    return y % (b / x), z % (-a / x) # modulus in this way so that y is positive and z is negative

Ruby, 83 bytes

I think there are a few ways to fine-tune and golf this solution, but I like it so far. Maybe I'll try an extended Euclidean algorithm solution next.

->x,y{a=b=0;y.downto(0).map{|u|(-x..0).map{|v|a,b=u,v if u*x+v*y==x.gcd(y)}};p a,b}

How it works

This code starts with a loop of u from y down to 0, with an inner loop of v from -x to 0, inside which we check every u and v if u*x+v*y == gcd(x, y). Since there could be multiple matches along the way (this uses a very exhaustive search), we start far from 0 so that when we get the last of the multiple matches, it's the one where |u| and |v| are closest to 0.

def bezout(x,y)
  a=b=0
  y.downto(0).each do |u|
    (-x..0).each do |v|
      if u*x + v*y == x.gcd(y)
        a,b=u,v
      end
    end
  end
  p a,b
end

Mathematica, 80 bytes

f@l_:=Mod@@NestWhile[{Last@#,{1,-Quotient@@(#.l)}.#}&,{{1,0},{0,1}},Last@#.l>0&]

Explanation:

Extended Euclidean algorithm is used here, in a Nest style. The method that the coefficients are stored in arrays makes it possible to use Dot.

Another possible representation is simply using symbolic expression, like u a - v b with {a->19, b->17}. Such representation makes use of the feature of Mathematica and is interesting, but it is much longer in bytes.

Test cases:

f[{5, 3}]              (* {2, -3} *)
f[{42, 12}]            (* {1, -3} *)
f[{11235813, 112358}]  (* {8643, -864301} *)

MATL, 37 40 bytes

ZdXK2Gw/:1G*1GK/t_w2$:XI2G*!+K=2#fIb)

Uses release (9.3.1), which is earlier than this challenge.

This is a brute force approach, so it may not work for large inputs.

Try it online! The online compiler is based on a newer release, but produces the same results.

Explanation

Zd            % implicitly input A and B. Compute their GCD. Call that C
XK            % copy C to clipboard K
2Gw/:         % vector [1, 2, ..., B/C]
1G*           % multiply that vector by A
1GK/t_w2$:    % vector [-A/C, -A/C+1 ..., A/C]
XI            % copy to clipboard I
2G*           % multiply that vector by B
!+            % all pairwise sums of elements from those vectors
K=2#f         % find row and column indices of sum that equals C
Ib)           % index second vector with column (indexing first vector with
              % row is not needed, because that vector is of the form [1, 2, ...])

Haskell, 51 bytes

a#b=[(u,-v)|v<-[1..],u<-[1..v],gcd a b==u*a-v*b]!!0

Usage example: 27998 # 6461 -> (3,-13).

This is a brute force approach which finds all combinations of u and v that are valid solutions ordered by u and picks the first one. This takes some time to run for large |v|.