# Precision while calculating Fourier transform

I'm trying to understand how precision works in Mathematica. Particularly I'm calculating discrete Fourier transform using the Fourier function and calculating it "manually". Here is my code:

range = Range[0, 2 Pi - 0.0000001, 2 Pi/8192] // N[#, 20] &;

wf = (Sin[256 #] + Cos[512 #]) &;

fft = Fourier[wf /@ range, FourierParameters -> {-1, 1}];
FFT[f_, wf_, range_] := Mean[Exp[I*(f - 1)*range]*(wf /@ range)];

fft[[{257, 513}]]
FFT[#, wf, range] & /@ {257, 513}


The result is:

{0.*10^-16 + 0.500000000000000 I, 0.500000000000000 + 0.*10^-16 I}
{0.*10^-17 + 0.5000000000000000 I, 0.5000000000000000 + 0.*10^-17 I}


So I expect to get both answers to be accurate up to 20 digits. But I'm puzzled by 0.*10^-16 and 0.*10^-17. What do they mean? how the exponent is chosen? why are they different for two methods?

What is the correct way to ask for 20 digits precision of the final result?

Edit: I'm closing this question as it seems that there is no definite answer and the main issue has nothing to do with Fourier transform, but just general problem of setting specific precision for numeric calculation.

• When carrying out a numerical calculation, precision will typically change (most often it will be reduced, as in your example). Even though your two functions are theoretically identical, they perform different numerical calculations and you cannot expect them to provide the same precision. Expressions of the form 0*10^n indicate that the number they correspond to is somewhere in the interval $[-10^n,10^n]$ – Stelios Jun 1 '15 at 18:18
• @Stelios Since I specifically requested 20 digits of precision why would I get two different answers? – BlacKow Jun 1 '15 at 18:25
• You requested 20 digits precision only for the elements of range. This does not guarantee that subsequent numerical manipulations of these elements will also give results with the same precision. Apparently, the two functions you consider do not perform the exact same computations in the exact same order. (Fourier employs a computational efficient algorithm with much better scaling properties than your FFT.) However, for all practical purposes the results are the same. – Stelios Jun 1 '15 at 18:56
• What is the correct way to ask for 20 digits precision of the final result? I understand that they are the same for practical reasons. I want to be sure I understand how accurate the answer is and how I can increase the accuracy. – BlacKow Jun 1 '15 at 19:03
• I guess the only thing you can do is set an appropriately large precision for the initial elements. I don't think there are any rules on how to do that, you just set the precision by trial and error. For your example, setting the initial precision to 26 gives output with precisions close to 21 and 23 for the two functions (as given by Precision function). Note that there are some numerical functions, e.g., FindRoot, for which you can specify the precision/accuracy goal of the output and issue an automatic warning when this is not achieved. – Stelios Jun 1 '15 at 19:21

You can disable error tracking as below (effectively working in fixed precision).

range = N[Range[0, 2 Pi - 0.0000001, 2 Pi/8192],20];

FFT[f_, wf_, range_] := Mean[Exp[I*(f - 1)*range]*(wf /@ range)];

Block[{$MinPrecision = 20,$MaxPrecision = 20},
FFT[#, wf, range] & /@ {257, 513}
]

(* Out[360]= {-5.7615176179584663593*10^-40 + 0.50000000000000000000 I,
0.50000000000000000000 + 8.8117221835711259773*10^-40 I} *)


Same story if you use built-in Fourier.

Block[{$MinPrecision = 20,$MaxPrecision = 20},
fft = Fourier[wf /@ range, FourierParameters -> {-1, 1}];
fft[[{257, 513}]]]

(* Out[367]= {-3.9424691514700133400*10^-40 + 0.50000000000000000000 I,
0.50000000000000000000 + 3.3631163143795609702*10^-40 I} *)


Whether this ise of fixed precision is a good idea or not depends on what you know about the error of the specific computation, and in particular whether it is something automated tracking will grossly exaggerate and that fixed precision will not mess up by much.

--- edit ---

I'll respond to the question in the comments here since this is a bit long for a comment itself.

There are ways, e.g. as above, to force a given precision on the output. As for obtaining it in a way that actually comes with a guarantee of correctness, in general you have to begin with sufficiently high input precision and work from there. You could, for example, do the computation in exact arithmetic and then apply N[#,20]& to the result. This might be very slow however.

Or do some estimate of the precision loss (not too hard in this case) and begin with enough precision to emerge with the desired precision. A comment by @Stelios describes this idea in a general way. For your DFT example it is fairly straightforward to get an overestimate. You have 2^13 elements of magnitudes between 0 and 2, the DFT matrix entries are all of unit size, so if the precision is epsilon then you lose, at most, 2^14*epsilon. That comes to a little over 4 digits. So start with 25 digits and you are guaranteed to have 20+ correct digits when finished (whether you use fixed precision as I did or let the precision tracking do its thing).

If you use built in Fourier the situation is slightly more complicated. A naive expectation, which I once entertained, is that the precision loss should be less than when doing a straight matrix-times-vector DFT. This is because fewer operations overall are performed. As it happens, precision loss is typically slightly greater and I believe that is due to the way the butterfly operations proceed. That is to say, I think each element of the result has been computed with more operations than is done in the DFT, and the total is far less because so much gets reused in the FFT butterfly process. I confess I never remember FFT innards well enough to be sure if I have the correct explanation though.

Another thing to consider is whether you are interested in the maximal error, as is (closely) approximated by precision tracking, or in the expected error e.g. if the input error comes from iid random error (say normally distributed). If the latter then you can do the computation in fixed precision and adjust the resulting precision downward by an a priori estimate of the expected loss. For this example I think it becomes something like 2^7*epsilon, that is to say, the expected loss factor in front of epsilon is around the square root of the worst-case loss.

--- end edit ---

• Ok. That is interesting. Thank you for pointing out the way to use fixed precision. Obviously the final precision depends on algorithm in this approach. Are you saying that there is no way to request the precision of final result? Even when the calculation is deterministic. – BlacKow Jun 1 '15 at 21:23
• Re obtaining a desired precision, see edit. – Daniel Lichtblau Jun 2 '15 at 15:12
• I totally understand what you are saying. I wonder why Mathematica can't do this for me? I ask for 20 digits, it checks its calculation and gives me the answer with requested precision. Instead I have to figure out what is the precision loss myself and start with higher precision to accommodate my needs. As I mentioned in my edit, I believe it has nothing to do with Fourier transform. – BlacKow Jun 2 '15 at 15:22
• Isn't Mathematica now using FFTW under the hood these days? – J. M.'s discontentment Jun 2 '15 at 15:22
• As for automating the precision output, one of the methods I mentioned does exactly this. If you give an exact expression and request 20 digits (via N[#,20]& say) then Mathematica will do the estimation work to make sure the result is as desired. But this might be slow. What it cannot do is give 20 validated correct digits if you do not begin with sufficiently many e.g. if the input is not exact. In the example, starting with 20 (or even 22) digits means the accumulation of error will be too much to provide a result with 20 validated digits. – Daniel Lichtblau Jun 2 '15 at 15:59