% Chapter 4: three experiments with floating-point arithmetic.
% Run ch04_roundoff from the folder containing this file.
% The extra printed digits expose rounding; they do not add accuracy.

clear;

% 1. Representation and precision.
a = 0.1;
b = 0.2;
c = 0.3;
s = a + b;
one_single = single(1.0);
delta_single = single(1.0e-8);
single_sum = one_single + delta_single;
double_sum = 1.0 + 1.0e-8;

fprintf('1. Representation and precision\n');
fprintf('%-30s%24.16e\n', '0.1 + 0.2:', s);
fprintf('%-30s%24.16e\n', '0.3:', c);
fprintf('%-30s%24.16e\n', 'Difference:', s-c);
fprintf('%-30s%24.16e\n', 'Gap above 1 (single):', double(eps(single(1))));
fprintf('%-30s%24.16e\n', 'Gap above 1 (double):', eps(1));
fprintf('%-30s%24.16e\n', '1 + 1e-8 (single):', double(single_sum));
fprintf('%-30s%24.16e\n', '1 + 1e-8 (double):', double_sum);

% 2. Cancellation. The stable formula is a reference, not exact arithmetic.
fprintf('\n2. Cancellation: sqrt(1+x) - 1\n');
fprintf('%12s%24s%24s\n', 'x', 'Direct formula', 'Stable formula');
for k = 2:2:18
    x = 10.0^(-k);
    direct = sqrt(1.0+x) - 1.0;
    stable = x / (sqrt(1.0+x) + 1.0);
    fprintf('%12.2e%24.16e%24.16e\n', x, direct, stable);
end

% 3. Order of summation. Store each intermediate result explicitly.
a = 1.0e16;
b = 1.0;
c = -1.0e16;
temp = a + b;
first = temp + c;
temp = a + c;
second = temp + b;

fprintf('\n3. Order of summation\n');
fprintf('%-30s%24.16e\n', '(a + b) + c:', first);
fprintf('%-30s%24.16e\n', '(a + c) + b:', second);
fprintf('%-30s%24.16e\n', 'Exact sum:', 1.0);