[Rd] fft fails for lengths 392, 588, 968, 980 .... (PR#1429)

[brahm@alum.mit.edu]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._