function [Gp, Rp] = estgr (gamma, delta, G, d, alpha, p)

% ------------------------------------------------------------------------
% Estimate lower/upper bounds for the expression d'G (G'G + alpha I)^p G'd
% using Gauss/Gauss-Radau quadrature rules.
%
% References:
% [1] G. H. GOLUB and U. von MATT,
%     "Generalized Cross-Validation for Large-Scale Problems", 1996,
%     TR-96-28
%
% [2] G. H. GOLUB and U. von MATT,
%     "Tikhonov Regularization for Large Scale Problems", 1997, SCCM-97-03
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

n         = length       (gamma);
normd     = norm         (d,2);
[Q, u, v] = lbidiagqr    (gamma, delta, sqrt(alpha));
Bk        = sparse       (1:n,   1:n, gamma, n+1, n) + ...
            sparse       (2:n+1, 1:n, delta, n+1, n);

b         = lbidiagqtx   (Q, eye(2*n+1, 1)*normd);
xi        = ubidiagsolve (u, v, b(1:n));

if     p == -1, Gp = Bk *xi; Gp = Gp(1)*normd;
elseif p == -2, Gp = xi'*xi;
elseif p == -3
  b   = lbidiagqtx   (Q, [zeros(n+1,1); xi/sqrt(alpha)]);
  eta = ubidiagsolve (u, v, b(1:n));
  Gp  = xi'*eta;
end

% compute \tilde{U}_k
[Q, u, v] = lbidiagqr (gamma, delta, 0) ;
u(n) = 0;

[Q, u, v] = ubidiagqr    (u, v, sqrt(alpha));
b         = ubidiagqtx   (Q, [zeros(n,1); eye(n,1)]);
eta       = ubidiagsolve (u, v, b(1:n));
zeta      = lbidiagsolve (u, v, eta);
Rp        = norm(G'*d)^2*norm (zeta, 2)^2/alpha;
