θ = 0.1
A = [ cos(θ) sin(θ); -sin(θ) cos(θ)]
N = 100
p = rand(N)
b = [1,0]

# compute xᴺ
function recurrence(A, b, p)
    x = b
    for i = 1:length(p)
        x = A * x + [0,p[i]]
    end
    return x
end

recurrence(A, b, p)

# compute g(x) = x[2]^2 from last iteration, and dg/dp, by adjoint method
function recurrence_adjoint(A, b, p)
    X = Array(Float64, 2, length(p)+1)
    X[:,1] = b
    for i = 1:length(p)
        X[:,i+1] = A * X[:,i] + [0,p[i]]
    end
    λ = [0, 2*X[2,end]]
    dgdp = Array(Float64, length(p))
    for i = length(p):-1:1
        dgdp[i] = λ[2]
        λ = A'λ
    end
    return X[2,end]^2, dgdp
end

recurrence_adjoint(A, b, p)

function recurrence_FD(A, b, p, δ=1e-5)
    g = recurrence(A, b, p)[2]^2
    dgdp = Array(Float64, length(p))
    for i = 1:length(p)
        pold = p[i]
        p[i] = pold + δ
        g₊ = recurrence(A, b, p)[2]^2
        p[i] = pold - δ
        g₋ = recurrence(A, b, p)[2]^2
        dgdp[i] = (g₊ - g₋) / (2δ)
        p[i] = pold
    end
    return g, dgdp
end

norm(recurrence_FD(A, b, p)[2] - recurrence_adjoint(A, b, p)[2])

using PyPlot
δ = logspace(-15, -1, 100)
loglog(δ, [norm(recurrence_adjoint(A, b, p)[2] - recurrence_FD(A, b, p, d)[2]) for d in δ], ".-")
xlabel("finite-difference step δ")
ylabel("norm of error in dg/dp")