Let's write our own DFT function for a comparison

I've taken the code for this from the documentation for Fourier. Let's fist compute the values symbolically

(* Calculate the DFT symbolically *)

x = ls; n = Length[x]; u /: Subscript[u, r_] := x[[r]];

FullSimplify[Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]]

{0, 3/2, 0, 0, 0, 0, 0, 0, 3/2}

We can run this same code with different precisions and compare to the "correct" results above.

(*With precision = 8*)

Block[{$MaxPrecision = 8, $MinPrecision = 8},

x = SetPrecision[ls, 8]; n = Length[x]; u /: Subscript[u, r_] := x[[r]];

Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]

]

{-4.4041855*10^-35, 1.5000000, 1.8255231*10^-34, 2.2020928*10^-35,

1.8255231*10^-34, 1.8255231*10^-34, 2.2020928*10^-35, 1.8255231*10^-34, 1.5000000}

(*With boring old machine precision*)

x = N[ls]; n = Length[x]; u /: Subscript[u, r_] := x[[r]];

Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]

{-1.11022*10^-16, 1.5 + 0. I, -9.25186*10^-17 + 0. I,

0. + 0. I, -1.11022*10^-16 + 0. I, -1.11022*10^-16 + 0. I,

0. + 0. I, -9.25186*10^-17 + 0. I, 1.5 + 0. I}

So the built-in compiled version and a version we built in Mathematica do the same thing. If we run with a bunch of different precisions, we can see how the results differ based on the accuracy of the numbers used (X axis is precision used. Y axis is how off the actual results the calculations end up being. Closer to 0 is better ):

Higher precision is generally good, but answering "Why is not always the case that that is true" can be difficult. After all, if I asked "why does 0.1+0.2 give 0.30000000000000004 in floating point computation?" what answer could be given except to walk step-by-step through how the addition plays out in floating point in all the painful detail?

Higher precision usually means better answers, yes, but there's nothing that guarantees this that I know of.