|
Packit |
67cb25 |
/* fft/factorize.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_fft_complex.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "factorize.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fft_complex_factorize (const size_t n,
|
|
Packit |
67cb25 |
size_t *nf,
|
|
Packit |
67cb25 |
size_t factors[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t complex_subtransforms[] =
|
|
Packit |
67cb25 |
{7, 6, 5, 4, 3, 2, 0};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* other factors can be added here if their transform modules are
|
|
Packit |
67cb25 |
implemented. The end of the list is marked by 0. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int status = fft_factorize (n, complex_subtransforms, nf, factors);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fft_halfcomplex_factorize (const size_t n,
|
|
Packit |
67cb25 |
size_t *nf,
|
|
Packit |
67cb25 |
size_t factors[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t halfcomplex_subtransforms[] =
|
|
Packit |
67cb25 |
{5, 4, 3, 2, 0};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int status = fft_factorize (n, halfcomplex_subtransforms, nf, factors);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fft_real_factorize (const size_t n,
|
|
Packit |
67cb25 |
size_t *nf,
|
|
Packit |
67cb25 |
size_t factors[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t real_subtransforms[] =
|
|
Packit |
67cb25 |
{5, 4, 3, 2, 0};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int status = fft_factorize (n, real_subtransforms, nf, factors);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fft_factorize (const size_t n,
|
|
Packit |
67cb25 |
const size_t implemented_subtransforms[],
|
|
Packit |
67cb25 |
size_t *n_factors,
|
|
Packit |
67cb25 |
size_t factors[])
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t nf = 0;
|
|
Packit |
67cb25 |
size_t ntest = n;
|
|
Packit |
67cb25 |
size_t factor;
|
|
Packit |
67cb25 |
size_t i = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("length n must be positive integer", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
factors[0] = 1;
|
|
Packit |
67cb25 |
*n_factors = 1;
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* deal with the implemented factors first */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (implemented_subtransforms[i] && ntest != 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
factor = implemented_subtransforms[i];
|
|
Packit |
67cb25 |
while ((ntest % factor) == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ntest = ntest / factor;
|
|
Packit |
67cb25 |
factors[nf] = factor;
|
|
Packit |
67cb25 |
nf++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
i++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* deal with any other even prime factors (there is only one) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
factor = 2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while ((ntest % factor) == 0 && (ntest != 1))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ntest = ntest / factor;
|
|
Packit |
67cb25 |
factors[nf] = factor;
|
|
Packit |
67cb25 |
nf++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* deal with any other odd prime factors */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
factor = 3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (ntest != 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
while ((ntest % factor) != 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
factor += 2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ntest = ntest / factor;
|
|
Packit |
67cb25 |
factors[nf] = factor;
|
|
Packit |
67cb25 |
nf++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that the factorization is correct */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t product = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < nf; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
product *= factors[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (product != n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("factorization failed", GSL_ESANITY);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*n_factors = nf;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
fft_binary_logn (const size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t ntest ;
|
|
Packit |
67cb25 |
size_t binary_logn = 0 ;
|
|
Packit |
67cb25 |
size_t k = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (k < n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
k *= 2;
|
|
Packit |
67cb25 |
binary_logn++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ntest = (1 << binary_logn) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n != ntest )
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return -1 ; /* n is not a power of 2 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return binary_logn;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|