Numerical Error source when dealing with integer series

Question

I am currently trying to compute the value of the first Fibonacci number recursively. the idea is as follow:

• Compute $f_{n}$ and $f_{n-1}$ for $n = 2,...,100$,
• Compute $f_k$ for $k = n−2, n−3, \dots, 0$,
• Call $\hat{f}_0(n)$ the value of $f_0$ computed as a function of $n$.

Here is the Matlab code I used:

N = 101;
n = linspace(0, 100, N);

% Fbonacci sequence generation
fn = fibonacci(n+1);

f0 = zeros(1, 99);

% backward computation of \hat{f_0} computation
for i = 2:100
    fi = fn(i+1);
    fi_1 = fn(i);
    fk = cat(2, zeros(1, i - 1), [fi_1, fi]);
    
    for k = i-2:-1:0
        a = fk(k + 3);
        b = fk(k + 2);
        fk(k+1) = a - b;
    end
    f0(i-1) = fk(1);
end

figure(2)
f1 = semilogy(linspace(2, 100, 99), abs(f0));
hold on
f2 = semilogy([0, 100], [1, 1]);
hold on
f3 = semilogy([0, 100], [2/eps, 2/eps]);
grid on

xlabel('$n$', 'Interpreter','latex')
legend([f1, f2, f3], {'$\hat{f}_0$', '$1$: analytical value', '$2/eps$'}, 'Interpreter', 'latex');

I can't find an explication to the errors found on the value of $\hat{f}_0$ in the second plot. Fibonacci numbers are integers, and integers are supposed to be represented exactly in machines, or am I wrong? Because in this case, when computing the value of $\hat{f}_0$ recursively, there is no rounding or truncation. The value of $f(100)$ is of order 20, which is far away from $realmax$. I don't see any source of errors.

Any ideas?

Answer 1

Integers can be represented "exactly" in Matlab, however a lot of times Matlab will choose to work with floating point double precision, which may not necessarily represent all integers exactly (IIRC it is exact up to about $2^{53}$).

Indeed, even when you try typing an "integer" into matlab, it will tell you that it defaults to doubles:

>> class(1)

ans =

    'double'

You can try to force Matlab to use ints by casting it to a fixed-width integer, for example using int64(your_value), or the symbolic library wide integer sym(your_value)

When I tried forcing the fibonacci function to use fixed-width integer math, it complains that the function expects a floating point type or a symbol. Here are the results to calling fibonacci with various type arguments:

>> fibonacci(int64(101))
Error using fibonacci (line 16)
Invalid data type. Argument must be double, single, or sym.

>> fibonacci(101)

ans =

   5.7315e+20

>> fibonacci(sym(101))
 
ans =
 
573147844013817084101

>> fibonacci(101)-fibonacci(sym(101))
 
ans =
 
-14533

The fact that Matlab can't compute fibonacci(101) exactly with doubles is not unexpected; the answer requires a 69-bit mantissa, while doubles only have a 53-bit mantissa.