A = {{a₀, a₁}, {a₁, a₂}}, b = {b₀, b₁}, x₀ = {0, 0}

Step 1: r₀ = b - Ax₀; p₀ = r₀ (No expansion)

r₀[0] = b₀

r₀[1] = b₁

Step 2: alpha₀ = <r₀, r₀> / <p₀, Ap₀> (No expansion)

<p₀, Ap₀> = b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁)

alpha₀ = (b₀²+b₁²) / (b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁))

Step 3: r₁ = r₀ - alpha₀*Ap₀ (No expansion)

r₁[0] = (b₀(b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁)) - (b₀²+b₁²)(a₀b₀+a₁b₁)) / (b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁))

r₁[1] = (b₁(b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁)) - (b₀²+b₁²)(a₁b₀+a₂b₁)) / (b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁))

Step 4: x₁ = x₀ + alpha₀*p₀ (No expansion)

x₁[0] = (b₀(b₀² + b₁²)) / (b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁))

x₁[1] = (b₁(b₀² + b₁²)) / (b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁))

Step 5: beta = rsnew / rsold = <r₁, r₁> / <r₀, r₀> (No expansion)

To make it simple to read, let's devote m = b₀(a₀b₀ + a₁b₁) + b₁(a₁b₀ + a₂b₁)

rsnew =

[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²
------------------------------------------------------------------------
m²

beta =

[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²
-------------------------------------------------------------------------
(b₀²+b₁²)m²

Step 6: p₁ = r₁ +beta*p₀ (No expansion)

p₁[0] =

m(b₀²+b₁²)[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)] + b₀{[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²}
------------------------------------------------------------------------------------------------------------------------------------
(b₀²+b₁²)m²

p₁[1] =

m(b₀²+b₁²)[b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)] + b₁{[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²}
------------------------------------------------------------------------------------------------------------------------------------
(b₀²+b₁²)m²

Step 7: alpha₁ = <r₁, r₁> / <p₁, Ap₁> (No expansion)

<p₁, Ap₁> = p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])

alpha₁ =

[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²
------------------------------------------------------------------------------------------------------------
m²[p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])]

Step 8: r₂ = r₁ - alpha₁*Ap₁ (Need to expand for cancellation)

r₂[0] = r₁[0] - alpha₁*(A[0][0]*p1[0]+A[0][1]*p1[1])

r₂[0] =

m[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)][p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])] - (a₀p₁[0]+a₁p₁[1]){[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²}
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m²[p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])]

r₂[1] =

Step 9: x₂ = x₁ + alpha₁*p₁

x₂[0] = x₁[0] + alpha₁*p₁[0]

x₂[0] =

mb₀(b₀²+b₁²)[p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])] + p₁[0]{[b₀m - (b₀²+b₁²)(a₀b₀+a₁b₁)]² + [b₁m - (b₀²+b₁²)(a₁b₀+a₂b₁)]²}
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m²[p₁[0](a₀p₁[0]+a₁p₁[1]) + p₁[1](a₁p₁[0]+a₂p₁[1])]

x₂[1] =

assertion: (Need to expand for cancellation)

bncg[0] = A[0][0]*x[0] + A[0][1]*x[1] = a₀(a₂b₀ - a₁b₁) / (a₀a₂ - a₁²) + a₁(-a₁b₀ + a₀b₁) / (a₀a₂ - a₁²) = b₀(a₀a₂ - a₁²) / (a₀a₂ - a₁²) = b₀

bncg[1] = A[1][0]*x[0] + A[1][1]*x[1] = a₁(a₂b₀ - a₁b₁) / (a₀a₂ - a₁²) + a₂(-a₁b₀ + a₀b₁) / (a₀a₂ - a₁²) = b₁(a₀a₂ - a₁²) / (a₀a₂ - a₁²) = b₁

b[0] = b₀

b[1] = b₁

assert(bncg[i] == b[i])

END