[SciPy-user] fft numerical precision

Ryan Krauss ryanlists at gmail.com
Wed Feb 15 16:27:14 CST 2006


Sorry to bring this back up again (this thread has been dormant since
12/21/05), but I think this is the key road block between me and
really understanding this problem (and I think it went over my head
when it was originally posted).  If double precision numbers use 52
bits for the mantissa, what does it mean to say that 2**-52 is the
approximate fractional precision of double precision floats?  How does
this affect optimization routines like fzero and fmin so that their
accuracy seems limited to 2**-52 when the smallest possible floating
point number should be more like 2**-1074.  Does this somehow relate
to how floating point subtraction is done (as relates to my post from
earlier today).

Thanks again,

Ryan


=========================
Hi Ryan,

Note that 1e-16 is the approximate _fractional_ precision of double
precision floats.  You can get absolute numbers much smaller than that.

There are a bunch of good references to floating point accuracy etc.
Some of them are linked here:
http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html

Note that the FFTW website lists FFT accuracies in both single and
double precision here:

http://www.fftw.org/accuracy/

Scott
- Show quoted text -

On Wednesday 21 December 2005 02:53 pm, Ryan Krauss wrote:
> I have run into this as well.  1e-16 seems large to me for double
> precision.  Robert can you explain this a bit more please.  Is there
> an IEEE spec or something that specifies how much of the 64bits is
> for exponent and how much is for the mantissa (I think that is the
> right word)?   I was playing with some FORTRAN code and it seems like
> there was a big difference with complex vs. double complex.  It seems
> like 1e-16 was the magnitude floor of complex and double complex was
> better.
>
> Ryan
>
> On 12/21/05, Robert Kern <robert.kern at gmail.com> wrote:
> > Darren Dale wrote:
> > > I am running into problems related, I think, to numerical
> > > precision of fast Fourier transforms. If I Fourier transform a
> > > gaussian distribution:
> > >
> > > absolute(fft(stats.norm.pdf(linspace(-10,10,512), loc=0,
> > > scale=1)))
> > >
> > > I find a floor of about 1e-16. Does anyone know of a way to
> > > improve the precision?
> >
> > 1e-16 is the best you can do with double precision.
> >
> > --
> > Robert Kern
> > robert.kern at gmail.com
> >
> > "In the fields of hell where the grass grows high
> >  Are the graves of dreams allowed to die."
> >   -- Richard Harter
> >
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user at scipy.net
> > http://www.scipy.net/mailman/listinfo/scipy-user
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user

--
Scott M. Ransom            Address:  NRAO
Phone:  (434) 296-0320               520 Edgemont Rd.
email:  sransom at nrao.edu             Charlottesville, VA 22903 USA
GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989

On 12/21/05, Robert Kern <robert.kern at gmail.com> wrote:
> Ryan Krauss wrote:
> > I have run into this as well.  1e-16 seems large to me for double
> > precision.  Robert can you explain this a bit more please.  Is there
> > an IEEE spec or something that specifies how much of the 64bits is for
> > exponent and how much is for the mantissa (I think that is the right
> > word)?   I was playing with some FORTRAN code and it seems like there
> > was a big difference with complex vs. double complex.  It seems like
> > 1e-16 was the magnitude floor of complex and double complex was
> > better.
>
> IEEE-754 double precision has 11 bits for exponent, 1 bit for sign, and 52 for
> the rest of the mantissa. 2.**-52 ~= 2.22e-16
>
> Goldberg, David. What Every Computer Scientist Should Know About Floating Point
> Arithmetic. 1991.
> http://www.physics.ohio-state.edu/~dws/grouplinks/floating_point_math.pdf
>
> --
> Robert Kern
> robert.kern at gmail.com
>
> "In the fields of hell where the grass grows high
>  Are the graves of dreams allowed to die."
>   -- Richard Harter
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
>



More information about the SciPy-user mailing list