+
+ % The recursive definition of r_next is prone to accumulate
+ % roundoff error. When sqrt(n) divides k, we recompute the
+ % residual to minimize this error. This modification was suggested
+ % by the second reference.
+ if (mod(k, sqrt_n) == 0)
+ r_next = Q*x_next - b;
+ else
+ r_next = rk + (alpha_k * Qdk);
+ end
+