      SUBROUTINE FFT(g,p)
*     複素数の入力データの組 g の高速フーリエ変換
*     変換値は g に戻される
      IMPLICIT NONE
      REAL*8 Wp(2), angle, factor(2), pi, temp(2), g(2,*)
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER L, N, N2, del_i, i, ip, j, jmax, k, p, Itmp1_
      N = 2**p                              ! データ点の数
      N2 = N/2
*     ビット順反転の方法による入力データの並べ換え
      j = 1
      DO i = 1, N-1
*        g(i,1) は f((i-1)*del_t) の実数部分
*        g(i,2) は f((i-1)*del_t) の虚数部分
*        データが実数ならば，g(i,2) = 0
         IF(i .LT. j) THEN
            temp(1) = g(1,j)
            temp(2) = g(2,j)
            g(1,j) = g(1,i)
            g(2,j) = g(2,i)
            g(1,i) = temp(1)
            g(2,i) = temp(2)
         END IF
         k = N2
         DO WHILE(k .LT. j)
            j = j-k
            k = k/2
         END DO
         j = j + k
      END DO
*     ダニエルソン-ランチョスのアルゴリズムの開始
      jmax = 1
      DO L = 1, p
         del_i = 2*jmax
         Wp(1) = 1                          ! Wp を W^0 で初期化
         Wp(2) = 0
         angle = pi/jmax
         factor(1) = cos(angle)             ! W^p の新しい値と古い値の比    
         factor(2) = -sin(angle)
         DO j = 1, jmax
            DO i = j, N, del_i
*              長さ 2^L の変換を計算
               ip = i + jmax
               temp(1) = g(1,ip)*Wp(1) - g(2,ip)*Wp(2)
               temp(2) = g(1,ip)*Wp(2) + g(2,ip)*Wp(1)
               g(1,ip) = g(1,i) - temp(1)
               g(2,ip) = g(2,i) - temp(2)
               g(1,i) = g(1,i) + temp(1)
               g(2,i) = g(2,i) + temp(2)
            END DO
*           新しい W^p を計算
            temp(1) = Wp(1)*factor(1) - Wp(2)*factor(2)
            temp(2) = Wp(1)*factor(2) + Wp(2)*factor(1)
            DO Itmp1_ = 1, 2
               Wp(Itmp1_) = temp(Itmp1_)
            END DO
         END DO
         jmax = del_i
      END DO
      END
