R devel - Apr 2002 - fft fails for lengths 392, 588, 968, 980 .... (PR#1429)

R 1.4.1, Linux and Windows

for(i in 1:1000) {
    X <- rnorm(i)
    XX <- fft(fft(X), inverse=T)/i
    if(max(Mod(XX-X)) > 1e-10) print(i)
}

[1] 392
[1] 588
[1] 968
[1] 980

and I then get a segfault during gc().  The answers are way off, with
imaginary parts 1e10 or more.

These numbers are all multiples of 7^2 or 11^2.

(Based on a report to R-help
Date: Thu, 28 Mar 2002 09:37:34 +0100
From: Gabriel Fricout <Gabriel.FRICOUT@cmm.ensmp.fr>
)


-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Gabriel Fricout <Gabriel.FRICOUT@cmm.ensmp.fr> reported a bug in fft(),
and
Brian D. Ripley <ripley@stats.ox.ac.uk> narrowed it down to specific
values of
N (the data vector length).  For N <= 1000, BDR showed the failures occur at:
  N = 392 =  7^2 * 2^2 * 2
  N = 588 =  7^2 * 2^2 * 3
  N = 968 = 11^2 * 2^2 * 2
  N = 980 =  7^2 * 2^2 * 5

BDR's test is as follows:
  R> for (i in 1:1000) {
  R>   X <- rnorm(i)
  R>   XX <- fft(fft(X), inverse=T)/i
  R>   if(max(Mod(XX-X)) > 1e-10) print(i)
  R> }

In src/appl/fft.c, lines 820-822 set "maxf" in subroutine fft_factor. 
"maxf"
is supposed to be the "maximum factor size", and the algorithm
basically grabs
the last non-squared factor and the last squared factor (there are
"kt" squared
factors):

  maxf = nfac[m_fac-kt-1];
  if (kt > 0)
    maxf = imax2(nfac[kt-1], maxf);

But the last squared factor (nfac[kt-1]) is not necessarily the largest, since
the squared factors are added to nfac in this order: 4, 3,5,7,..., 2 (see lines
769-791).  Thus I believe we need to add these two lines after line 822:

  if (kt > 1) maxf = imax2(nfac[kt-2], maxf);
  if (kt > 2) maxf = imax2(nfac[kt-3], maxf);

I don't know for sure that the second line is necessary (I saw no error for
N = 1152 = 4^2 * 3^2 * 2^2 * 2).  But I did rebuild R-devel (2002-04-03) with
these two lines added, and BDR's test ran without errors.  I submit this as
the
bug fix for PR#1429.
-- 
                              -- David Brahm (brahm@alum.mit.edu)

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._