Loss of Significance
Contents
This script demonstrates loss of significant digits and how to fix it for the specific example of evaluation of
for increasingly smaller values of .
Number of trials
N = 15;
Result using double-precision
disp('Result using double precision') for n=0:N x = 10^(-n); f = sqrt(x^2+1)-1; disp(sprintf('x=%f, sqrt(x^2+1)-1 = %e',x,f)) end
Result using double precision x=1.000000, sqrt(x^2+1)-1 = 4.142136e-001 x=0.100000, sqrt(x^2+1)-1 = 4.987562e-003 x=0.010000, sqrt(x^2+1)-1 = 4.999875e-005 x=0.001000, sqrt(x^2+1)-1 = 4.999999e-007 x=0.000100, sqrt(x^2+1)-1 = 5.000000e-009 x=0.000010, sqrt(x^2+1)-1 = 5.000000e-011 x=0.000001, sqrt(x^2+1)-1 = 5.000445e-013 x=0.000000, sqrt(x^2+1)-1 = 4.884981e-015 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000
Using single-precision
Here is where we see the problems
pause disp('Result using single precision') for n=0:N x = 10^(-n); f = single(sqrt(x^2+1))-1; disp(sprintf('x=%f, sqrt(x^2+1)-1 = %e',x,f)) end
Result using single precision x=1.000000, sqrt(x^2+1)-1 = 4.142135e-001 x=0.100000, sqrt(x^2+1)-1 = 4.987597e-003 x=0.010000, sqrt(x^2+1)-1 = 4.994869e-005 x=0.001000, sqrt(x^2+1)-1 = 4.768372e-007 x=0.000100, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000010, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000001, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000 x=0.000000, sqrt(x^2+1)-1 = 0.000000e+000
Using single-precision, but fixed
pause disp('Result using single precision, but fixed') for n=0:N x = 10^(-n); f = single(x^2)/(single(sqrt(x^2+1))+1); disp(sprintf('x=%f, x^2/(sqrt(x^2+1)+1) = %e',x,f)) end
Result using single precision, but fixed x=1.000000, x^2/(sqrt(x^2+1)+1) = 4.142135e-001 x=0.100000, x^2/(sqrt(x^2+1)+1) = 4.987562e-003 x=0.010000, x^2/(sqrt(x^2+1)+1) = 4.999875e-005 x=0.001000, x^2/(sqrt(x^2+1)+1) = 4.999999e-007 x=0.000100, x^2/(sqrt(x^2+1)+1) = 5.000000e-009 x=0.000010, x^2/(sqrt(x^2+1)+1) = 5.000000e-011 x=0.000001, x^2/(sqrt(x^2+1)+1) = 5.000000e-013 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-015 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-017 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-019 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-021 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-023 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-025 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-027 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-029 x=0.000000, x^2/(sqrt(x^2+1)+1) = 5.000000e-031