% 4-22-2008
% Round-off Error example
clear all
close all

% Determine the roots of x^2 + 2x + 1 = 0 (Wilkinson's Polynomial)
a =[2432902008176640000
    -8752948036761600000
    13803759753640704000     
    -12870931245150988800
    8037811822645051776
    -3599979517947607200
    1206647803780373360
    -311333643161390640
    63030812099294896
    -10142299865511450
    1307535010540395
    -135585182899530
    11310276995381
    -756111184500
    40171771630
    -1672280820
    53327946
    -1256850
    20615
    -210
    1]
a=flipud(a);

plot(roots(a),zeros(20),'k.','MarkerSize',30)
hold on

% Plot 100 random variations on the wilkinson polynomial
for j=1:100
    a =[2432902008176640000
    -8752948036761600000
    13803759753640704000     
    -12870931245150988800
    8037811822645051776
    -3599979517947607200
    1206647803780373360
    -311333643161390640
    63030812099294896
    -10142299865511450
    1307535010540395
    -135585182899530
    11310276995381
    -756111184500
    40171771630
    -1672280820
    53327946
    -1256850
    20615
    -210
    1];
    a=flipud(a);
    % Perturb all 20 coefficients randomaly by a factor no greater than
    % 10^(-10)
    for i=1:20
        a(i) = a(i)*(1 + 10^(-10)*unifrnd(-1,1));
    end
    plot(roots(a),'b.','MarkerSize',10)
end
hold off