Trying to build a function to compute divided difference for arbitrary list of points

Well, instead of evaluating an individual entry as $ f_{ij} $ once at a time, here I provide an implementation generating the whole upper-triangular $ f $ matrix. So here we go:

Clear[fmat]
fmat[pts : {{_, _} ..}] := Module[{n = Length[pts], kernels, xs, ys, denomenators, entries},
                                  {xs, ys} = pts\[Transpose]; 
                                  kernels = NestList[Insert[#, 0, 2] &, {-1, 1}, n - 2];
                                  denomenators = ListCorrelate[#, xs] & /@ kernels;
                                  entries = FoldList[Differences[#]/#2 &, ys, denomenators];
                                  SparseArray[MapIndexed[Band[Prepend[#2, 1]] -> # &, entries], n]
                                 ]

To check its validity, one can see a symbolic result and compare it with that on the Wiki page:

Then below calculation is trivial

Also, I have tested on lists of random reals with a larger shape, say $ 1000\times 2 $. The efficiency of generating the $ 1000\times 1000 $ $ f $ matrix is acceptable, at least.

---

P.S. I learnt that the calculation of divided difference is related to Newton interpolation. So if your purpose is to do so interpolations, you can directly use built-in functions like Interpolation, InterpolatingPolynomial, etc.

ClearAll[dividedDif]
dividedDif = NestList[#/(#[[1, 1]] + #[[-1, -1]] &[Denominator /@ (List @@ #)]) & /@ 
     Differences[#] &, Divide @@ (Differences /@ #), Length[#[[1]]] - 2][[All, 1 ]] &;

Prepend[dividedDif @ #, #[[1, 1]]] & @ {{r, s, t, u, v}, {a, b, c, d, e}}