# Why does Fourier give a large error in MachinePrecision?

Consider this code:

f[x_] = Cos[2 π x];

N0 = 10;

ls = Most @ Array[f, N0, {0, 1}]
(*
{1, Cos[(2 π)/9], Sin[π/18], -(1/2), -Cos[π/9], -Cos[π/9],
-(1/2), Sin[π/18], Cos[(2 π)/9]}
*)


Since the function is a periodic function, I would expect the Fourier transform give only two non-zero coefficients. Using Fourier I get

Fourier[ls]
(*
{-6.4763*10^-17 + 0. I, 1.5 + 5.55112*10^-17 I, 1.60247*10^-17 + 2.77556*10^-17 I,
3.23815*10^-17 - 8.01234*10^-18 I, -1.60247*10^-17 - 2.77556*10^-17 I,
-1.60247*10^-17 + 2.77556*10^-17 I, 3.23815*10^-17 + 8.01234*10^-18 I,
1.60247*10^-17 - 2.77556*10^-17 I, 1.5 - 5.55112*10^-17 I}
*)


there is an error about 10^-17.

However, if I change the precision, I get

Fourier[SetPrecision[ls, 8]]
(*
{0.*10^-8, 1.5000000 + 0.*10^-8 I, 0.*10^-8 + 0.*10^-8 I, 0.*10^-8 + 0.*10^-8 I,
0.*10^-8 + 0.*10^-8 I, 0.*10^-8 + 0.*10^-8 I, 0.*10^-8 + 0.*10^-8 I,
0.*10^-8 + 0.*10^-8 I, 1.5000000 + 0.*10^-8 I}
*)

Fourier[SetPrecision[ls, 16]]
(*
{0.*10^-16, 1.500000000000000 + 0.*10^-16 I, 0.*10^-16 + 0.*10^-16 I,
0.*10^-16 + 0.*10^-16 I, 0.*10^-16 + 0.*10^-16 I, 0.*10^-16 + 0.*10^-16 I,
0.*10^-16 + 0.*10^-16 I, 0.*10^-16 + 0.*10^-16 I, 1.500000000000000 + 0.*10^-16 I}
*)

Fourier[SetPrecision[ls, 32]]
(*
{0.*10^-32, 1.5000000000000000000000000000000 + 0.*10^-32 I,0.*10^-32 + 0.*10^-32 I,
0.*10^-32 + 0.*10^-32 I, 0.*10^-32 + 0.*10^-32 I, 0.*10^-32 + 0.*10^-32 I,
0.*10^-32 + 0.*10^-32 I, 0.*10^-32 + 0.*10^-32 I,
1.5000000000000000000000000000000 + 0.*10^-32 I}
*)


So why I get the error only in the machine precision case? Are these errors in Fourier[ls] coming from Mathematica or related to how numbers are represented in computers?

Edit:

If we fix the precision in the calculation, we can see that the errors are smaller compared to Machine precision calculation.

Block[{$MaxPrecision = 8,$MinPrecision = 8}, Fourier[SetPrecision[ls, 8]]]
(* {-4.3981121*10^-35, 1.5000000, 1.8230938*10^-34, 2.1989581*10^-35, 1.8230938*10^-34, 1.8230938*10^-34, 2.1989581*10^-35, 1.8230938*10^-34, 1.5000000} *)

Block[{$MaxPrecision = 16,$MinPrecision = 16}, Fourier[SetPrecision[ls, 16]]]
(* {-4.406340574057345*10^-35, 1.500000000000000, 1.816158363592948*10^-34, 2.203268244891241*10^-35, 1.816158363592948*10^-34, 1.816158363592948*10^-34, 2.203268244891241*10^-35, 1.816158363592948*10^-34, 1.500000000000000} *)

-
In absolute terms, it isn't a large error; it represents about 1 unit in the last place, which is pretty good. It's just large in a relative sense because the values are supposed to be zero, but in this case I think there isn't any alternative to working in arbitrary precision so that you can (1) make the absolute value of 1 ULP arbitrarily small and (2) know whether the values returned have any significant digits or not. –  Oleksandr R. May 10 '14 at 19:07
@OleksandrR. How can we do (2)? And do you know why Fourier[SetPrecision[ls, 16]] seems give smaller errors? –  xslittlegrass May 10 '14 at 19:12
The errors are not necessarily smaller; it is just that Mma realises that the result has no significant digits and is zero to within the available precision. The actual error or uncertainty in the results is more or less the same whether you use MachinePrecision or $MachinePrecision (though actually the total accrued error seems a bit worse for the arbitrary-precision calculation). In Mma, using arbitrary precision automatically turns on precision tracking, so you accomplish both (1) and (2) at the same time. – Oleksandr R. May 10 '14 at 19:19 @OleksandrR. How can I see the error in the results are more or less the same? And how precision tracking could help to decrease the errors if the calculation are running at a fixed precision? – xslittlegrass May 10 '14 at 20:22 Compare the precision of the output with that of the input, and you'll see that about half a decimal digit is lost. Precision tracking does not decrease the error; it just causes e.g.$0.5 \pm 1.0\$ to be reported as 0 rather than 0.5. This is a distinction without a difference. –  Oleksandr R. May 10 '14 at 20:33

When you use machine precision, it is often a good idea to apply Chop to your result. In this case,

 {-6.4763*10^-17 + 0. I, 1.5 + 5.55112*10^-17 I, 1.60247*10^-17 + 2.77556*10^-17 I,
3.23815*10^-17 - 8.01234*10^-18 I, -1.60247*10^-17 - 2.77556*10^-17 I,
-1.60247*10^-17 + 2.77556*10^-17 I, 3.23815*10^-17 + 8.01234*10^-18 I,
1.60247*10^-17 - 2.77556*10^-17 I, 1.5 - 5.55112*10^-17 I} // Chop


gives

{0, 1.5, 0, 0, 0, 0, 0, 0, 1.5}


which is probably what you expected.

-
Yes I think this works for most of case, but I'm dealing with a problem (called high harmonic generation) which cares a lot on the small numbers. Sorry for the confusing part. –  xslittlegrass May 10 '14 at 18:32
@xslittlegrass. Then you probably can't use machine arithmetic. –  m_goldberg May 10 '14 at 18:46
So are you saying that the errors are from the they way that the numbers are presented in computers, and not from Mathematica? –  xslittlegrass May 10 '14 at 19:05
@xslittlegrass. Yes –  m_goldberg May 10 '14 at 19:39