Extended Euclidean Algorithm: backward vs. forward

I'm confused about how to do the extended algorithm. For example, here's the gcd part gcd(8000,7001)

$$\begin{align}8000 &= 7001\cdot1 + 999\\ 7001&=999\cdot 7+8\\ 999&=8\cdot 124+7\\ 8&=7\cdot 1+1\\ 7&=1\cdot 7+0\\ \gcd(8000,7001)&=1\end{align}$$

Now, the extended algorithm

$$\begin{align}1 &= 8 - 1\cdot7\\ &= 8 - 1\cdot(999 - 8\cdot124)\\ &= -1\cdot999 + 8\cdot125\end{align}$$

How do you get this 8*125 from the previous line?

$\begin{eqnarray}\!\text{By the distributive law}\ \ && \,8\:\!\ \overbrace{ -\, 1\cdot(999\,-\,8\cdot 124)}^{\textstyle -1\,(a\!-\!b) = -a\! +\! b\ }\\ &=&\ 8\cdot\color{#c00}1\, -\, 999\, +\, 8\cdot\color{#c00}{124}\\ &=&\ 8\cdot\color{#0a0}{ 125} - 999\ \ \,{\rm by}\ \ \color{#c00}{124 + 1} = \color{#0a0}{125}\end{eqnarray}$

This common "$\rm\color{#90f}{back}$-substitution" extended Euclidean gcd algorithm is notoriously error-prone. Better, this $\rm\color{#90f}{forward}$ method is simpler to compute and easier to remember. It keeps track of each remainder's expression as a linear combination of the gcd arguments. Executing it here we obtain that below, where each line $\,\ a\ \ b\ \ c\ \,$ means $\ a = 8000\, b + 7001\, c.\ $

$$\begin{array}{rrr} 8000 & 1 & 0\\ 7001 & 0 & 1\\ 999 & 1 & -1\\ 8 & -7& 8\\ -1 & \!\!\color{#c00}{876} & \!\!\!\color{#0a0}{-1001}\\ \end{array}$$

Here is another example computing $\rm\ gcd(141,19),\,$ with the equations written explicitly

$$\rm\begin{eqnarray} [\![1]\!]\ \ \ \, \color{#C00}{141}\!\ &=&\,\ \ 1&\cdot& 141\, +\ 0&\cdot& 19 \\ [\![2]\!]\quad\ \color{#C00}{19}\ &=&\,\ \ 0&\cdot& 141\, +\ 1&\cdot& 19 \\ \color{#940}{[\![1]\!]-7\,[\![2]\!]}\, \rightarrow\, [\![3]\!]\quad\ \ \ \color{#C00}{ 8}\ &=&\,\ \ 1&\cdot& 141\, -\ 7&\cdot& 19 \\ \color{#940}{[\![2]\!]-2\,[\![3]\!]}\,\rightarrow\,[\![4]\!]\quad\ \ \ \color{#C00}{3}\ &=& {-}2&\cdot& 141\, + 15&\cdot& 19 \\ \color{#940}{[\![3]\!]-3\,[\![4]\!]}\,\rightarrow\,[\![5]\!]\quad \color{#C00}{{-}1}\ &=&\,\ \ 7&\cdot& 141\, -\color{#0A0}{ 52}&\cdot& \color{#0A0}{19} \\ {\rm negating}\ \Rightarrow\ \ \ \ \ \ {1}\ &=& {-}7&\cdot& 141\, +\color{#0A0}{ 52}&\cdot& \color{#0A0}{19} \end{eqnarray}\qquad\qquad\qquad\quad$$

The prior Bezout equation $\Rightarrow 141^{-1}\equiv \color{c00}{-7}\pmod{\!19},\,$ & $\,\color{#0a0}{19^{-1}\!\equiv 52}\pmod{\!141}\,$ by reducing the Bezout equation $\bmod19\,$ and $\bmod 141\,$ resp., as explained here. Thus we see that using the extended Euclidean algorithm to compute the gcd Bezout equation yields one method of computing modular inverses (and fractions). See here and here for many more examples of this and various other methods.

$$ \begin{align} 8000&=7001\cdot1+999\\ 7001&=999\cdot7+8\\ 999&=8\cdot124+7\\ 8&=7\cdot1+1\\[12pt] 1 &=8-7\cdot1\\ &=8-1(999-8\cdot124)\\ &=8\cdot125-1\cdot999\\ &=125(7001-7\cdot999)-1(8000-7001\cdot1)\\ &=126\cdot7001-1\cdot8000-875\cdot999\\ &=126\cdot7001-1\cdot8000-875\cdot(8000-7001)\\ &=1001\cdot7001-876\cdot8000 \end{align} $$ As Bill Dubuque mentions, the back-substitution method is hard to follow, and therefore, error-prone. There is also Euclid-Wallis version of the Extended Euclidean Algorithm: $$ \begin{array}{r} &&1&7&124&1&7\\\hline \color{#C00000}{1}&0&1&-7&869&\color{#C00000}{-876}&7001\\ 0&\color{#00A000}{1}&-1&8&-993&\color{#00A000}{1001}&-8000\\ \color{#C00000}{8000}&\color{#00A000}{7001}&999&8&7&\color{#0000FF}{1}&0\\ \end{array} $$ Which gives $1001\cdot7001-876\cdot8000=1$