fourn.cpp

此程序为VC++常用数值算法这本书附赠的光盘中包含了本书中全部的源代码

void fourn(double data[],double nn[],int &ndim,int &isign)
{
	int ntot,idim,i1,i2,i3,i2rev,i3rev,k1,k2,nprev,n;
	double nrem,ip3,tempr,tempi,ibit,ip1,ip2,wi,wr,wpi,wpr;
	double ifp1,ifp2,wtemp,theta;
    ntot = 1;
    for (idim = 1; idim<=ndim; idim++)
	{
        ntot = ntot * nn[idim];
	}
    nprev = 1;
    for (idim = 1; idim<=ndim; idim++)
	{
        n = nn[idim];
        nrem = ntot / (n * nprev);
        ip1 = 2 * nprev;
        ip2 = ip1 * n;
        ip3 = ip2 * nrem;
        i2rev = 1;
        for (i2 = 1; i2<=ip2; i2+=ip1)
		{
            if( i2 < i2rev)
			{
                for (i1 = i2; i1<=i2 + ip1 - 2; i1+=2)
				{
                    for (i3 = i1;i3<=ip3;i3+=ip2)
					{
                        i3rev = i2rev + i3 - i2;
                        tempr = data[i3];
                        tempi = data[i3 + 1];
                        data[i3] = data[i3rev];
                        data[i3 + 1] = data[i3rev + 1];
                        data[i3rev] = tempr;
                        data[i3rev + 1] = tempi;
                    }
                }
            }
            ibit = ip2 / 2;
yi:         if( ibit >= ip1 && i2rev > ibit)
			{
                i2rev = i2rev - ibit;
                ibit = ibit / 2;
                goto yi;
			}
            i2rev = i2rev + ibit;
        }
        ifp1 = ip1;
er:     if (ifp1 < ip2)
		{
            ifp2 = 2 * ifp1;
            theta = isign * 6.28318530717959 / (ifp2 / ip1);
            wpr = -2.0 * (sin(0.5 * theta))*(sin(0.5 * theta));
            wpi = sin(theta);
            wr = 1.0;
            wi = 0.0;
            for (i3 = 1; i3<=ifp1; i3+=ip1)
			{
              for (i1 = i3;i1<=i3 + ip1 - 2;i1+=2)
			  {
                for(i2 = i1; i2<=ip3; i2+=ifp2)
				{
                  k1 = i2;
                  k2 = k1 + ifp1;
                  tempr=float(wr)*data[k2]-float(wi)*data[k2+1];
                  tempi=float(wr)*data[k2+1]+float(wi)*data[k2];
                  data[k2] = data[k1] - tempr;
                  data[k2 + 1] = data[k1 + 1] - tempi;
                  data[k1] = data[k1] + tempr;
                  data[k1 + 1] = data[k1 + 1] + tempi;
                }
              }
              wtemp = wr;
              wr = wr * wpr - wi * wpi + wr;
              wi = wi * wpr + wtemp * wpi + wi;
            }
            ifp1 = ifp2;
            goto er;
		}
        nprev = n * nprev;
	}
}