htswanalysis/MACS/lib/gsl/gsl-1.11/fft/factorize.c

   1 /* fft/factorize.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_errno.h>
  22 #include <gsl/gsl_fft_complex.h>
  23
  24 #include "factorize.h"
  25
  26 static int
  27 fft_complex_factorize (const size_t n,
  28                            size_t *nf,
  29                            size_t factors[])
  30 {
  31   const size_t complex_subtransforms[] =
  32   {7, 6, 5, 4, 3, 2, 0};
  33
  34   /* other factors can be added here if their transform modules are
  35      implemented. The end of the list is marked by 0. */
  36
  37   int status = fft_factorize (n, complex_subtransforms, nf, factors);
  38   return status;
  39 }
  40
  41 static int
  42 fft_halfcomplex_factorize (const size_t n,
  43                                size_t *nf,
  44                                size_t factors[])
  45 {
  46   const size_t halfcomplex_subtransforms[] =
  47   {5, 4, 3, 2, 0};
  48
  49   int status = fft_factorize (n, halfcomplex_subtransforms, nf, factors);
  50   return status;
  51 }
  52
  53 static int
  54 fft_real_factorize (const size_t n,
  55                         size_t *nf,
  56                         size_t factors[])
  57 {
  58   const size_t real_subtransforms[] =
  59   {5, 4, 3, 2, 0};
  60
  61   int status = fft_factorize (n, real_subtransforms, nf, factors);
  62   return status;
  63 }
  64
  65
  66 static int
  67 fft_factorize (const size_t n,
  68                    const size_t implemented_subtransforms[],
  69                    size_t *n_factors,
  70                    size_t factors[])
  71
  72 {
  73   size_t nf = 0;
  74   size_t ntest = n;
  75   size_t factor;
  76   size_t i = 0;
  77
  78   if (n == 0)
  79     {
  80       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  81     }
  82
  83   if (n == 1)
  84     {
  85       factors[0] = 1;
  86       *n_factors = 1;
  87       return 0;
  88     }
  89
  90   /* deal with the implemented factors first */
  91
  92   while (implemented_subtransforms[i] && ntest != 1)
  93     {
  94       factor = implemented_subtransforms[i];
  95       while ((ntest % factor) == 0)
  96         {
  97           ntest = ntest / factor;
  98           factors[nf] = factor;
  99           nf++;
 100         }
 101       i++;
 102     }
 103
 104   /* deal with any other even prime factors (there is only one) */
 105
 106   factor = 2;
 107
 108   while ((ntest % factor) == 0 && (ntest != 1))
 109     {
 110       ntest = ntest / factor;
 111       factors[nf] = factor;
 112       nf++;
 113     }
 114
 115   /* deal with any other odd prime factors */
 116
 117   factor = 3;
 118
 119   while (ntest != 1)
 120     {
 121       while ((ntest % factor) != 0)
 122         {
 123           factor += 2;
 124         }
 125       ntest = ntest / factor;
 126       factors[nf] = factor;
 127       nf++;
 128     }
 129
 130   /* check that the factorization is correct */
 131   {
 132     size_t product = 1;
 133
 134     for (i = 0; i < nf; i++)
 135       {
 136         product *= factors[i];
 137       }
 138
 139     if (product != n)
 140       {
 141         GSL_ERROR ("factorization failed", GSL_ESANITY);
 142       }
 143   }
 144
 145   *n_factors = nf;
 146
 147   return 0;
 148 }
 149
 150
 151 static int
 152 fft_binary_logn (const size_t n)
 153 {
 154   size_t ntest ;
 155   size_t binary_logn = 0 ;
 156   size_t k = 1;
 157
 158   while (k < n)
 159     {
 160       k *= 2;
 161       binary_logn++;
 162     }
 163
 164   ntest = (1 << binary_logn) ;
 165
 166   if (n != ntest )
 167     {
 168       return -1 ; /* n is not a power of 2 */
 169     }
 170
 171   return binary_logn;
 172 }
 173
 174
 175
 176