θ = 0.1
A = [ cos(θ) sin(θ); -sin(θ) cos(θ)]
N = 100
p = [1,2,3,4,5]
b = [1,0]

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

recurrence(A, b, p)

# compute G(x) = ∑ x[2]^2, and dG/dp, by adjoint method
function recurrence_adjoint(A, b, p)
    X = Array(Float64, 2, length(p)+1)
    X[:,1] = b
    G = X[2,1]^2
    for i = 1:length(p)
        X[:,i+1] = A * X[:,i] + [0,p[i]]
        G += X[2,i+1]^2
    end
    λ = [0, 2*X[2,end]]
    dGdp = Array(Float64, length(p))
    for i = length(p):-1:1
        dGdp[i] = λ[2]
        λ = A'λ + [0, 2*X[2,i]]
    end
    return G, dGdp
end

recurrence_adjoint(A, b, p)

function recurrence_FD(A, b, p, δ=1e-5)
    G = recurrence(A, b, p)
    dGdp = Array(Float64, length(p))
    p′ = copy!(similar(p, Float64), p) # modify a copy of p rather than p
    for i = 1:length(p)
        p′[i] = p[i] + δ
        G₊ = recurrence(A, b, p′)
        p′[i] = p[i] - δ
        G₋ = recurrence(A, b, p′)
        dGdp[i] = (G₊ - G₋) / (2δ)
        p′[i] = p[i]
    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")