Prepare a MATLAB/Octave script that sums n random, single-precision floating-poi
ID: 2259218 • Letter: P
Question
Prepare a MATLAB/Octave script that sums n random, single-precision floating-point numbers xi, uniformly distributed in the interval [0, 1]. Sum the numbers in each of the following ways (using single-precision variables unless specifically indicated otherwise):
(a) Sum the numbers in the order in which they were generated, using a double-precision variable to accumulate the sum.
(b) Sum the numbers in the order in which they were generated, this time using a single-precision accumulator.
(c) Use the following compensated summation algorithm, again using only single precision, to sum the numbers in the order in which they were generated:
s = x1
c=0
for i = 1 to n do
y = xi c t=s+y
c = (t s) y s=t
end for
1
(d) Sum the numbers in order of increasing magnitude (using the builtin sort function).
(e) Sum the numbers in decreasing magnitude (reversing the order obtained in the previous part).
Once you check for correctness, extend the previous program to answer the following:
(i) Experiment with different values of n and chart your results. Note that you may need to use
fairly large values of n to see substantial differences.
(ii) How do the methods rank in terms of accuracy, and why?
(iii) How do the methods compare in cost? Time the routines for a large n and provide a com- parison.
(iv) Briefly explain how/why the compensated summation algorithm works.
Explanation / Answer
MATLAB code
close all
clear
clc
n = 10;
x = single(rand(1,n));
% Part (a)
sum_double = double(0);
for i=1:n
sum_double = sum_double + x(i);
end
fprintf('(a) Sum [double] = %f ', sum_double);
% Part (b)
sum_single = single(0);
for i=1:n
sum_single = sum_single + x(i);
end
fprintf('(b) Sum [single] = %f ', sum_single);
% Part (c)
s = x(1);
c = single(0);
for i=2:n
y = x(i) - c;
t = s + y;
c = (t - s) - y;
s = t;
end
fprintf('(c) Sum = %f ', s);
% Part (d)
x = sort(x);
sum_sort = single(0);
for i=1:n
sum_sort = sum_sort + x(i);
end
fprintf('(d) Sum = %f ', sum_sort);
output
(a) Sum [double] = 5.322150
(b) Sum [single] = 5.322150
(c) Sum = 5.322150
(d) Sum = 5.322150