dpfafft.c

该程序是用vc开发的对动态数组进行管理的DLL

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
PFAFFT - Functions to perform Prime Factor (PFA) FFT's, in place

npfa_d		return valid n for fcomplex-to-fcomplex PFA
npfar_d		return valid n for real-to-fcomplex/fcomplex-to-real PFA
npfao_d		return optimal n for fcomplex-to-fcomplex PFA
npfaro_d		return optimal n for real-to-fcomplex/fcomplex-to-real PFA
pfacc_d		1D PFA fcomplex to fcomplex
pfacr_d		1D PFA fcomplex to real
pfarc_d		1D PFA real to fcomplex
pfamcc_d		multiple PFA fcomplex to real
pfa2cc_d		2D PFA fcomplex to fcomplex
pfa2cr_d		2D PFA fcomplex to real
pfa2rc_d		2D PFA real to fcomplex

*****************************************************************************
Function Prototypes:
int npfa_d (int nmin);
int npfao_d (int nmin, int nmax);
int npfar_d (int nmin);
int npfaro_d (int nmin, int nmax);
void pfacc_d (int isign, int n, real_complex z[]);
void pfacr_d (int isign, int n, real_complex cz[], real rz[]);
void pfarc_d (int isign, int n, real rz[], real_complex cz[]);
void pfamcc_d (int isign, int n, int nt, int k, int kt, real_complex z[]);
void pfa2cc_d (int isign, int idim, int n1, int n2, real_complex z[]);
void pfa2cr_d (int isign, int idim, int n1, int n2, real_complex cz[], real rz[]);
void pfa2rc_d (int isign, int idim, int n1, int n2, real rz[], real_complex cz[]);

*****************************************************************************
npfa_d:
Input:
nmin		lower bound on returned value (see notes below)

Returned:	valid n for prime factor fft

******************************************************************************
npfao_d
Input:
nmin		lower bound on returned value (see notes below)
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for prime factor fft

******************************************************************************
npfar_d
Input:
nmin		lower bound on returned value

Returned:	valid n for real-to-real_complex/real_complex-to-real prime factor fft

*****************************************************************************
npfaro_d:
Input:
nmin		lower bound on returned value
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for real-to-real_complex/real_complex-to-real prime factor fft

******************************************************************************
pfacc_d:
Input:
isign		sign of isign is the sign of exponent in fourier kernel
n		length of transform (see notes below)
z		array[n] of real_complex numbers to be transformed in place

Output:
z		array[n] of real_complex numbers transformed

******************************************************************************
pfacr_d:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform (see notes below)
cz          array[n/2+1] of real_complex values (may be equivalenced to rz)

Output:
rz          array[n] of real values (may be equivalenced to cz)

******************************************************************************
pfarc_d:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform; must be even (see notes below)
rz          array[n] of real values (may be equivalenced to cz)

Output:
cz          array[n/2+1] of real_complex values (may be equivalenced to rz)

******************************************************************************
pfamcc_d:
Input:
isign       	sign of isign is the sign of exponent in fourier kernel
n           	number of real_complex elements per transform (see notes below)
nt          	number of transforms
k           	stride in real_complex elements within transforms
kt          	stride in real_complex elements between transforms
z           	array of real_complex elements to be transformed in place

Output:
z		array of real_complex elements transformed

******************************************************************************
pfa2cc_d:
Input:
isign       	sign of isign is the sign of exponent in fourier kernel
idim        	dimension to transform, either 1 or 2 (see notes)
n1          	1st (fast) dimension of array to be transformed (see notes)
n2          	2nd (slow) dimension of array to be transformed (see notes)
z           	array[n2][n1] of real_complex elements to be transformed in place

Output:
z		array[n2][n1] of real_complex elements transformed

******************************************************************************
pfa2cr_d:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
idim        dimension to transform, which must be either 1 or 2 (see notes)
n1          1st (fast) dimension of array to be transformed (see notes)
n2          2nd (slow) dimension of array to be transformed (see notes)
cz          array of real_complex values (may be equivalenced to rz)

Output:
rz          array of real values (may be equivalenced to cz)

******************************************************************************
pfa2rc_d:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
idim        dimension to transform, which must be either 1 or 2 (see notes)
n1          1st (fast) dimension of array to be transformed (see notes)
n2          2nd (slow) dimension of array to be transformed (see notes)
rz          array of real values (may be equivalenced to cz)

Output:
cz          array of real_complex values (may be equivalenced to rz)

******************************************************************************
Notes:
Table of valid n and cost for prime factor fft.  For each n, cost
was estimated to be the inverse of the number of ffts done in 1 sec
on an IBM RISC System/6000 Model 320H, by Dave Hale, 08/04/91.
(Redone by Jack Cohen for 15 sec to rebuild NTAB table on advice of
David and Gregory Chudnovsky, 05/03/94).
Cost estimates are least accurate for very small n.  An alternative method
for estimating cost would be to count multiplies and adds, but this
method fails to account for the overlapping of multiplies and adds
that is possible on some computers, such as the IBM RS/6000 family.

npfa_d:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.

npfao_d:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.
If nmin does not exceed 720720, then the returned n will not be 
less than nmin.  The optimal n is chosen to minimize the estimated
cost of performing the fft, while satisfying the constraint, if
possible, that n not exceed nmax.

npfar and npfaro:
Current implemenations of real-to-real_complex and real_complex-to-real prime 
factor ffts require that the transform length n be even and that n/2 
be a valid length for a real_complex-to-real_complex prime factor fft.  The 
value returned by npfar satisfies these conditions.  Also, see notes 
for npfa.

pfacc:
n must be factorable into mutually prime factors taken 
from the set {2,3,4,5,7,8,9,11,13,16}.  in other words,
	n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
where
	0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
is required for pfa to yield meaningful results.  this
restriction implies that n is restricted to the range
	1 <= n <= 720720 (= 5*7*9*11*13*16)

pfacr:
Because pfacr uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfarc:
Because pfarc uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfamcc:
To perform a two-dimensional transform of an n1 by n2 real_complex array 
(assuming that both n1 and n2 are valid "n"), stored with n1 fast 
and n2 slow:
    pfamcc(isign,n1,n2,1,n1,z); (to transform 1st dimension)
    pfamcc(isign,n2,n1,n1,1,z); (to transform 2nd dimension)

pfa2cc:
Only one (either the 1st or 2nd) dimension of the 2-D array is transformed.

If idim equals 1, then n2 transforms of n1 real_complex elements are performed; 
else, if idim equals 2, then n1 transforms of n2 real_complex elements are 
performed.

Although z appears in the argument list as a one-dimensional array,
z may be viewed as an n1 by n2 two-dimensional array:  z[n2][n1].

Valid n is computed via the "np" subroutines.

To perform a two-dimensional transform of an n1 by n2 real_complex array 
(assuming that both n1 and n2 are valid "n"), stored with n1 fast 
and n2 slow:  pfa2cc(isign,1,n1,n2,z);  pfa2cc(isign,2,n1,n2,z);

pfa2cr:
If idim equals 1, then n2 transforms of n1/2+1 real_complex elements to n1 real 
elements are performed; else, if idim equals 2, then n1 transforms of n2/2+1 
real_complex elements to n2 real elements are performed.

Although rz appears in the argument list as a one-dimensional array,
rz may be viewed as an n1 by n2 two-dimensional array:  rz[n2][n1].  
Likewise, depending on idim, cz may be viewed as either an n1/2+1 by 
n2 or an n1 by n2/2+1 two-dimensional array of real_complex elements.

Let n denote the transform length, either n1 or n2, depending on idim.
Because pfa2rc uses pfa2cc to do most of the work, n must be even 
and n/2 must be a valid length for pfa2cc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfa2rc:
If idim equals 1, then n2 transforms of n1 real elements to n1/2+1 real_complex 
elements are performed; else, if idim equals 2, then n1 transforms of n2 
real elements to n2/2+1 real_complex elements are performed.

Although rz appears in the argument list as a one-dimensional array,
rz may be viewed as an n1 by n2 two-dimensional array:  rz[n2][n1].  
Likewise, depending on idim, cz may be viewed as either an n1/2+1 by 
n2 or an n1 by n2/2+1 two-dimensional array of real_complex elements.

Let n denote the transform length, either n1 or n2, depending on idim.
Because pfa2rc uses pfa2cc to do most of the work, n must be even 
and n/2 must be a valid length for pfa2cc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

******************************************************************************
References:  
Temperton, C., 1985, Implementation of a self-sorting
in-place prime factor fft algorithm:  Journal of
Computational Physics, v. 58, p. 283-299.

Temperton, C., 1988, A new set of minimum-add rotated
rotated dft modules: Journal of Computational Physics,
v. 75, p. 190-198.

Press et al, 1988, Numerical Recipes in C, p. 417.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 04/27/89
Revised by Baoniu Han to handle double precision. 12/14/98
*****************************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

typedef double real;
typedef dcomplex real_complex;


#define NTAB 240
static struct {
int n;  real c;
} nctab[NTAB] = {
{       1, 0.000052 },
{       2, 0.000061 },
{       3, 0.000030 },
{       4, 0.000053 },
{       5, 0.000066 },
{       6, 0.000067 },
{       7, 0.000071 },
{       8, 0.000062 },
{       9, 0.000079 },
{      10, 0.000080 },
{      11, 0.000052 },
{      12, 0.000069 },
{      13, 0.000103 },
{      14, 0.000123 },
{      15, 0.000050 },
{      16, 0.000086 },
{      18, 0.000108 },
{      20, 0.000101 },
{      21, 0.000098 },
{      22, 0.000135 },
{      24, 0.000090 },
{      26, 0.000165 },
{      28, 0.000084 },
{      30, 0.000132 },
{      33, 0.000158 },
{      35, 0.000138 },
{      36, 0.000147 },
{      39, 0.000207 },
{      40, 0.000156 },
{      42, 0.000158 },
{      44, 0.000176 },
{      45, 0.000171 },
{      48, 0.000185 },
{      52, 0.000227 },
{      55, 0.000242 },
{      56, 0.000194 },
{      60, 0.000215 },
{      63, 0.000233 },
{      65, 0.000288 },
{      66, 0.000271 },
{      70, 0.000248 },
{      72, 0.000247 },
{      77, 0.000285 },
{      78, 0.000395 },
{      80, 0.000285 },
{      84, 0.000209 },
{      88, 0.000332 },
{      90, 0.000321 },
{      91, 0.000372 },
{      99, 0.000400 },
{     104, 0.000391 },
{     105, 0.000358 },
{     110, 0.000440 },
{     112, 0.000367 },
{     117, 0.000494 },
{     120, 0.000413 },
{     126, 0.000424 },
{     130, 0.000549 },
{     132, 0.000480 },
{     140, 0.000450 },
{     143, 0.000637 },
{     144, 0.000497 },
{     154, 0.000590 },
{     156, 0.000626 },
{     165, 0.000654 },
{     168, 0.000536 },
{     176, 0.000656 },
{     180, 0.000611 },
{     182, 0.000730 },
{     195, 0.000839 },
{     198, 0.000786 },
{     208, 0.000835 },
{     210, 0.000751 },
{     220, 0.000826 },
{     231, 0.000926 },
{     234, 0.000991 },
{     240, 0.000852 },
{     252, 0.000820 },
{     260, 0.001053 },
{     264, 0.000987 },
{     273, 0.001152 },
{     280, 0.000952 },
{     286, 0.001299 },
{     308, 0.001155 },
{     312, 0.001270 },
{     315, 0.001156 },
{     330, 0.001397 },
{     336, 0.001173 },
{     360, 0.001259 },
{     364, 0.001471 },
{     385, 0.001569 },
{     390, 0.001767 },
{     396, 0.001552 },
{     420, 0.001516 },
{     429, 0.002015 },
{     440, 0.001748 },
{     455, 0.001988 },
{     462, 0.001921 },
{     468, 0.001956 },
{     495, 0.002106 },
{     504, 0.001769 },
{     520, 0.002196 },
{     528, 0.002127 },
{     546, 0.002454 },
{     560, 0.002099 },
{     572, 0.002632 },
{     585, 0.002665 },
{     616, 0.002397 },
{     624, 0.002711 },
{     630, 0.002496 },
{     660, 0.002812 },
{     693, 0.002949 },
{     715, 0.003571 },
{     720, 0.002783 },
{     728, 0.003060 },
{     770, 0.003392 },
{     780, 0.003553 },
{     792, 0.003198 },
{     819, 0.003726 },
{     840, 0.003234 },
{     858, 0.004354 },
{     880, 0.003800 },
{     910, 0.004304 },
{     924, 0.003975 },
{     936, 0.004123 },
{     990, 0.004517 },
{    1001, 0.005066 },
{    1008, 0.003902 },
{    1040, 0.004785 },
{    1092, 0.005017 },
{    1144, 0.005599 },
{    1155, 0.005380 },
{    1170, 0.005730 },
{    1232, 0.005323 },
{    1260, 0.005112 },
{    1287, 0.006658 },
{    1320, 0.005974 },
{    1365, 0.006781 },
{    1386, 0.006413 },
{    1430, 0.007622 },
{    1456, 0.006679 },
{    1540, 0.007032 },
{    1560, 0.007538 },
{    1584, 0.007126 },
{    1638, 0.007979 },
{    1680, 0.007225 },
{    1716, 0.008961 },
{    1820, 0.008818 },
{    1848, 0.008427 },
{    1872, 0.009004 },
{    1980, 0.009398 },
{    2002, 0.010830 },
{    2145, 0.012010 },
{    2184, 0.010586 },
{    2288, 0.012058 },
{    2310, 0.011673 },
{    2340, 0.011700 },
{    2520, 0.011062 },
{    2574, 0.014313 },
{    2640, 0.013021 },
{    2730, 0.014606 },
{    2772, 0.013216 },
{    2860, 0.015789 },
{    3003, 0.016988 },
{    3080, 0.014911 },
{    3120, 0.016393 },
{    3276, 0.016741 },
{    3432, 0.018821 },
{    3465, 0.018138 },
{    3640, 0.018892 },
{    3696, 0.018634 },
{    3960, 0.020216 },
{    4004, 0.022455 },
{    4095, 0.022523 },
{    4290, 0.026087 },
{    4368, 0.023474 },
{    4620, 0.024590 },
{    4680, 0.025641 },
{    5005, 0.030303 },
{    5040, 0.025253 },
{    5148, 0.030364 },
{    5460, 0.031250 },
{    5544, 0.029412 },