Quadruple precision floating-point accumulation of a double precision vector p is emulated numerically in Matlab:
function t = orosumvec(p, recurs) x = cumsum(p); z = x - p; u = (z - x) + p; v = [0 x(1:end-1)] - z; if ~recurs t = sum([u v x(end)]); else t = orosumvec([u v x(end)], 0); end end
This is verified by Advanpix Multiprecision Toolbox.
Matlab's built-in variable precision arithmetic, vpa() from Mathworks Symbolic Math Toolbox,
produces erroneous results in 2018. One recursive call is necessary and sufficient.
Customer support service by UserEcho