% matlab script to illustrate the use of compensated summation on
% vectors of random numbers

clear all, close all

% initialize array of single-precision random numbers in (0,1)
m    = 20;
rng('default');
x    = rand(m,1,'single');

% (a) sum variables in the order in which they are listed
%     sum is accumulated in double precision

sumD = sum(x,'double');

% (b) sum variables in the order in which they are listed
%     sum is accumulated in single precision

sumS = sum(x, 'native');

% (c) sum variables in the order in which they are listed using
%     compensated summation
%     sum is accumulated in double precision

sumC = compSum(x);

% (d) sum variables in order of ascending magnitude
%     sum is accumulated in single precision

[xAsc,I]   = sort(abs(x), 'ascend');
xAsc       = x(I);
sumAsc     = sum(xAsc, 'native');

% (e) sum variables in order of descending magnitude
%     sum is accumulated in single precision

xDes   = flipud(xAsc);
sumDes = sum(xDes, 'native');

% let's look at the results:
% we can assume sumD is the exact answer

disp(['Error in single-precision sum: ', num2str(abs(sumD-sumS))])
disp(['Error in compensated sum: ', num2str(abs(sumD-sumC))])
disp(['Error in ascending-order sum: ', num2str(abs(sumD-sumAsc))])
disp(['Error in descending-order sum: ', num2str(abs(sumD-sumDes))])
