FFTPACKPaul N. Swarztrauber
勝手な訳 桂田 祐史
1996 年 10 月 15 日, 2016 年 12 月 31 日
この文書の誤りを見つけた方は、勝手な訳者まで連絡して頂けると幸いです。連絡先: katurada@meiji.ac.jp

目次
1

勝手な訳注

2

2

コピーライトと内容一覧

33

各関数の一覧
3.1 rﬀti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 rﬀtf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 rﬀtb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 ezﬀti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 ezﬀtf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 ezﬀtb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 sinti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8 sint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 costi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.10 cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.11 sinqi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 sinqf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 sinqb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.14 cosqi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.15 cosqf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.16 cosqb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.17 cﬀti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.18 cﬀtf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.19 cﬀtb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
4
6
7
7
9
10
11
12
12
13
14
15
16
17
18
19
19
20

A dﬀtpack の C 言語プログラムからの利用
A.1 筆者が利用している環境の説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 dﬀtpack の入手とインストール . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 dﬀtpack の C 言語プログラムからの利用 . . . . . . . . . . . . . . . . . . . . . . . .
A.4 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
22
22
23
23

1

A.5 おまけ: サンプル・プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.5.1 周期境界条件の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.5.2 Dirichlet 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.5.3 Neumann 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
24
27
29

勝手な訳注

1

この文書は、有名な FFT ライブラリィ・パッケージである

FFTPACKby Paul N. Swarztrauberの INDEX ファイル1 を翻訳したものに、訳者が C プログラマー向けの注を補ったものです。
元々のプログラムは、FORTRAN の伝統的なプログラミング・スタイルで書かれていて、例えば配列の添字の下限は 1 であると仮定してあります。FFT を学んだ人は、時には添字の下限を 0 にした方が、すっきりした表現になることを知っているでしょう。C 言語の場合は、FORTRAN とは異なり、配列の添字は 0 から始まるので、自然にプログラムを書けば、すっきりとした表現になるわけです。そうした場合、元の文書に書いてある説明とは食い違いが出て来ますので、そこのところは訳者が説明を補っています。なお、一部では、配列の添字の下限を 1 からにした方が自然な場
合もあります。訳者の方針は、C 言語からの利用の仕方を説明するところでは、 1 要素程度のメモ
リーの無駄は無視して、なるべく自然な表現になるようにすることでした。なお、FORTRAN 77 では、配列の添字の下限を好きなところから始めるように宣言できますから、無駄がなく、しかも自然な書き方をすることが可能です。その場合にも、 C 言語向けの注は参考になると思います。

2016 年 12 月追記メモ
ほぼ放置していたのですが、昔書いたプログラムを動かすことになって、20 年ぶりに ﬀtpack を使うことになりました。そのため、この文書に手を入れる必要性が発生しました。
ただし、今回実際に使うのは、dﬀtpack2 です。199x 年当時は、オリジナルの FFTPACK を自分で倍精度化したライブラリィを作り、それを C 言語のプログラムから利用していました。その際、サブルーチンの名前をオリジナルと同じものにしていましたが、後で発見した dﬀtpack では、単精度版とは異なる名前がつけられています。
dﬀapack のサブルーチン名の規則



原則は、rﬀti → dﬀti, ezﬀti → dzﬀti, sinti → dsinti のように、先頭の文字を “d” にする。
例外は複素数バージョンで cﬀti → zﬀti のように、先頭の文字を “z” にする

もちろん、このようにする方がまともだと (今の私は) 考えます。そこで、今回はこの dﬀtpack に乗り換えて、独自倍精度化したものは捨てることにしました。
今回、C 言語の単精度浮動小数点型 float としてあるところを、この文書の以前の版では REALと書いてあったのは、実体が倍精度浮動小数点型 double であったからですが、今回は素直に floatと書くことにしました。
今回も実際の応用プログラムは C 言語で書くことにしたので、インターフェースを作る必要がありますが、今回はそれを fftpack-c-interface として (身の回り以外で使う人がいるとも思えないが) 公開することにしました。その使い方については付録 ??に説明します。
1

dﬀtpack の doc に相当する。

2

2

2

コピーライトと内容一覧* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *version 4

april 1985

a package of fortran subprograms for the fast fouriertransform of periodic and other symmetric sequencesbypaul n swarztrauber
national center for atmospheric researchboulder,colorado 80307

which is sponsored by the national science foundation* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *このパッケージは、以下に列挙するような、実および複素周期数列やある種の対称数列に対する
高速フーリエ変換のプログラムから構成されています (カッコ内は dﬀtpack のサブルーチン名)。

rﬀti (dﬀti) rﬀtf と rﬀtb を初期化する
rﬀtf (dﬀtf) 実周期数列の順変換
rﬀtb (dﬀtb) 実係数配列の逆変換
ezﬀti (dzﬀti) ezﬀtf と ezﬀtb を初期化する
ezﬀtf (dzﬀtf) 単純化された実周期の順変換
ezﬀtb (dzﬀtb) 単純化された実周期の逆変換
sinti (dsinti) sint を初期化する
sint (dsint) 実奇関数列のサイン変換
costi (dcosti) cost を初期化する
cost (dcost) 実偶関数列のコサイン変換
sinqi (dsinqi) sinqf と sinqb を初期化する
sinqf (dsinqf) 奇波数の順サイン変換
sinqb (dsinqb) sinqf の非正規化逆変換
cosqi (dsinqi) cosqf と cosqb の初期化
cosqf (dsinqf) 奇波数の順コサイン変換
cosqb (dsinqb) cosqf の非正規化逆変換
cﬀti (zﬀti) cﬀtf と cﬀtb の初期化
cﬀtf (zﬀtf) 複素周期数列の順変換
cﬀtb (zﬀtb) 複素周期数列の非正規化逆変換

3

3

各関数の一覧

3.1

rﬀti

機能と引数並びFORTRAN の場合
subroutine rﬀti(n, wsave)integer nreal wsave(*)
C の場合
void rﬀti(int n, float wsave[]);サブルーチン rﬀti は rﬀtf と rﬀtb の両者で使われる配列 wsave を初期化する。n の素因数分解と三角関数の表が計算され、wsave に格納される。
入力パラメーター

n 変換する数列の長さ
出力パラメーター

wsave 少なくとも長さ 2n + 15 の作業用配列。 n が変更されない限り、rﬀtf と rﬀtbの双方で同じ作業用配列を使うことが出来る。n の異なる値に対しては異なる配列
wsave が必要になる。wsave の内容は rﬀtf もしくは rﬀtb の呼び出しまで変更してはならない。

3.2

rﬀtf

機能と引数並びFORTRAN の場合
subroutine rﬀtf(n,r,wsave)integer nreal r(*),wsave(*)
C の場合
void rﬀtf(int n, float r[], float wsave[]);サブルーチン rﬀtf は実周期数列の Fourier 係数を計算する (Fourier 解析3 )。この変換は後述の出力パラメーター r のところで定義されるものである。

3

ここでは、周期関数の関数値の列から Fourier 係数を求めることを Fourier 解析 (Fourier Analysis)」と呼んでいる。

4

入力パラメーターn 変換される配列 r の長さ。その方法は n が小さな素数の積である時にとても効率的である。異なる作業用配列を与えるならば n を変えてもよい。
r 変換される数列を含む長さ n の実数型配列。
wsave rﬀtf を呼ぶプログラムにおける少なくとも長さ 2n+15 の作業用配列。配列 wsave
はサブルーチン rﬀti(n,wsave) を呼ぶことにより初期化されねばならない。異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化は nが変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。rﬀtf と rﬀtb では同じ配列 wsave が使える。
出力パラメーターr Fourier 係数が入ります。
FORTRAN の場合 まずn∑r(1) =r(i).

i=1

もしも n が偶数の場合 ℓ = n/2, n が奇数の場合 ℓ = (n + 1)/2 とおく。
k = 2, . . . , l に対して

r(2k − 2) =n∑(r(i) cosi=1r(2k − 1) =n∑

2π(k − 1)(i − 1)n(−r(i) sini=1

もしも n が偶数ならばr(n) =n∑)2π(k − 1)(i − 1)
n

,
)
.

(−1)i−1 r(i).

i=1

C の場合 まず
r[0] =

n−1
∑

r[i].

i=0

もしも n が偶数の場合 ℓ = n/2 − 1, n が奇数の場合 ℓ = (n − 1)/2 とおく。
k = 1, . . . , ℓ に対して
n−1
∑

2πki
r[i] cos,r[2k − 1] =ni=0r[2k] = −

n−1
∑
i=0

r[i] sin

2πki.
n

もしも n が偶数ならば

r[n − 1] =n−1∑(−1)i r[i].

i=0wsave rﬀtf もしくは rﬀtb を呼び出すまで壊してはならない計算結果を含む。
注意 この変換は正規化されていないので、rﬀtf に続いて rﬀtb を呼び出すと、入力列を n 倍することになる。

5

3.3

rﬀtb

機能と引数並びFORTRAN の場合
subroutine rﬀtb(n,r,wsave)integer nreal r(*), wsave(*)
C の場合
void rﬀtb(int n, float r[], float wsave[])
サブルーチン rﬀtb は、実周期数列をその Fourier 係数から計算する (Fourier synthesis)。この変換は後述の出力パラメーター r のところで定義される。
入力パラメーターn 変換される配列 r の長さ。その方法は n が小さな素数の積である時にとても効率的である。異なる作業用配列を与えるならば n を変えてもよい。
r 変換される数列を含む長さ n の実数配列。
wsave rﬀtb を呼ぶプログラムにおける少なくとも長さ 2n + 15 の作業用配列。配列
wsave はサブルーチン rﬀti(n,wsave) を呼ぶことにより初期化されねばならない。
異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーターr 関数値の列が入る。
FORTRAN の場合 偶数 n に対しては、i = 1, · · · , n に対して、
r(i) = r(1) + (−1)i−1 r(n)
)
n/2 (
∑
2π(k − 1)(i − 1)2π(k − 1)(i − 1)+
2r(2k − 2) cos− 2r(2k − 1) sin
.
n
n
k=2
奇数 n に対しては、i = 1, · · · , n に対して、(n+1)/2 (r(i) = r(1)+∑k=22π(k − 1)(i − 1)2π(k − 1)(i − 1)
2r(2k − 2) cos− 2r(2k − 1) sinnnC の場合 偶数 n に対しては、i = 0, · · · , n − 1 に対して、
n/2−1 (

r[i] = r[0] + (−1) r[n − 1] +
i

∑
k=1

2πki
2πki
− 2r[2k] sin
2r[2k − 1] cosnn

奇数 n に対しては、i = 0, · · · , n − 1 に対して、(n−1)/2 (r[i] = r[0] +

∑
k=1

2πki
2πki
− 2r[2k] sin
2r[2k − 1] cos
n
n

)
.

wsave rﬀtf もしくは rﬀtb を呼び出すまで壊してはならない計算結果を含む。
6

)
.

)
.

注意 この変換は正規化されていないので、rﬀtf に続いて rﬀtb を呼び出すと、入力列を n 倍することになる4 。

3.4

ezﬀti

機能と引数並びFORTRAN の場合
subroutine ezﬀti(n,wsave)integer nreal wsave(*)
C の場合
void ezﬀti(int n, float wsave[])
サブルーチン ezﬀti は ezﬀtf と ezﬀtb の両方で使われる配列 wsave を初期化する。n の素因数分解と三角関数の表が計算されて wsave に格納される。
入力パラメーター

n 変換する数列の長さ
出力パラメーター

wsave 少なくとも長さ 3n + 15 の作業用配列。n が変更されない限り ezﬀtf と ezﬀtbの双方で同じ作業用配列が使える。異なる n の値に対しては異なる配列 wsave が要求される。

3.5

ezﬀtf

機能と引数並びFORTRAN の場合
subroutine ezﬀtf(n,r,azero,a,b,wsave)integer nreal r(*),azero,a(*),b(*),wsave(*)
C の場合
void ezﬀtf(int n, float r[], float *azero, float a[], float b[], float wsave[])
サブルーチン ezﬀtf は実周期数列の Fourier 変換を計算する (Fourier 解析)。この変換は後述の
出力パラメーター azero, a, b のところで定義してある。 ezﬀtf は rﬀtf の、簡単だが、よりのろまなバージョンである。

4

これは誤植ではないかな。多分 n/2 倍ではないかと思う。

7

C プログラマーへのコメント 配列 r[], a[] については、0 から要素を詰めて格納するのが自然だ
が、sin の項の Fourier 係数を納める b[] だけは、添字は 1 から始めた方が自然だと思われる。そこで、以下では、
float r[MAXN], a[MAXN/2+1], b[MAXN/2+1], wsave[3*MAXN+15];のように宣言しておいて、

ezﬀtf(n, r, &a[0], a+1, b+1, wsave);あるいは

ezﬀtf(n, r, a, a+1, b+1, wsave);と呼び出して使うことを念頭に説明してある。b[0] は未使用で無駄になっているが、これが嫌ならば、

ezﬀtf(n, r, &a[0], a+1, b, wsave);のようにすればよい (宣言も b[MAXN/2] で良くなる)。実は

ezﬀtf(n, r, a, b, wsave);のように呼び出せるようにインターフェイスを作ることも可能だが、 FORTRAN の場合の呼び出し方との互換性を重視してやめにした。
入力パラメーターn 変換される配列 r の長さ。その方法は n が小さな素数の積である時にとても効率的である。
r 変換される数列を含む長さ n の実数型配列。r は破壊されない。
wsave ezﬀtf を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン ezﬀti(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。 ezﬀtf と ezﬀtb で同じ配列 wsave が使える。
出力パラメーター

azero FORTRAN の場合は

n∑r(i)/n. C の場合は

i=1

n−1
∑

r[i]/n.

i=0a,b Fourier 係数が入る。
FORTRAN の場合 偶数 n に対しては b(n/2) = 0 で a(n/2) =n∑(−1)i−1 r(i)/n.

i=1

偶数 n に対しては kmax = n/2 − 1, 奇数 n に対しては kmax = (n − 1)/2 と定義し、k = 1, · · · , kmax に対して
n∑22πk(i − 1)a(k) =
r(i) cos,nni=1b(k) =n∑2i=1

C の場合 偶数 n に対しては b[n/2] = 0 で a[n/2] =nr(i) sinn−1∑2πk(i − 1).
n(−1)i r[i]/n. 偶数 n に

i=0

対しては kmax = n/2 − 1, 奇数 n に対しては kmax = (n − 1)/2 と定義し、
k = 1, · · · , kmax に対して

2∑2πkia[k] =
r[i] cos
,
n i=0
n
n−1

8

2πki
2∑
r[i] sinb[k] =
.
n i=0
n
n−1

3.6

ezﬀtb

機能と引数並びFORTRAN の場合
subroutine ezﬀtb(n,r,azero,a,b,wsave)integer nreal r(*),azero,a(*),b(*),wsave(*)
C の場合
void ezﬀtb(int n, float r[], float *azero, float a[], float b[], float wsave[])
サブルーチン ezﬀtb は実周期数列の Fourier 変換を計算する (Fourier 同調)。この変換は後述の出力パラメーター r のところで定義してある。 ezﬀtb は rﬀtb の、簡単だが、よりのろまなバージョンである。

C プログラマーへのコメント まず ezﬀtf に対してのコメントはここでも有効である。そこで、以下では、
float r[MAXN], a[MAXN/2+1], b[MAXN/2+1], wsave[3*MAXN+15];のように宣言しておいて、

ezﬀtb(n, r, &a[0], a+1, b+1, wsave);あるいは

ezﬀtb(n, r, a, a+1, b+1, wsave);と呼び出して使うことを念頭に説明してある。b[0] は未使用で無駄になっているが、これが嫌ならば、

ezﬀtb(n, r, &a[0], a+1, b, wsave);のようにすればよい (宣言も b[MAXN/2] で良くなる)。なお、ezﬀtf とは異なり、a, b は入力であるから、&a[0] でなく、a[0] を引数として渡すようにすることも可能であったが (というか、C のプ
ログラムとしては、その方が自然)、ezﬀtf との対称性を重視して、このままにしてある。
入力パラメーターn 変換される配列 r の長さ。その方法は n が小さな素数の積である時にとても効率的である。
azero 定数 Fourier 係数a, b 残りの Fourier 係数を含む配列。これらの配列は破壊されない。これらの配列の長さは n が偶数か奇数かによる。 n が偶数の場合、n/2 だけの場所が必要である。 nが奇数の場合、(n − 1)/2 だけの場所が必要である。
wsave ezﬀtb を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン ezﬀti(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。 ezﬀtf と ezﬀtb で同じ配列 wsave が使える。
9

出力パラメーターr 関数値の列が入る。
FORTRAN の場合 n が偶数ならば kmax = n/2, n が奇数ならば kmax = (n −1)/2 と定義し、i = 1, · · · , n に対して
)
kmax
∑(
2πk(i − 1)2πk(i − 1)r(i) = azero ++ b(k) sin.
a(k) cosnnk=1
複素数を用いて書くと、j = 1, · · · , n に対して
kmax∑

r(j) =

c(k) exp

k=−kmax

2πik(j − 1)
,
n

ここで

a(k) − ib(k),2c(0) = azero, ,√ただし i = −1.

c(−k) = c(k) (k = 1, · · · , kmax),c(k) =*************** 振幅 – 位相による記述 ***********************FORTRAN の場合 i = 1, 2, · · · , n に対して
r(i) = azero +

kmax
∑

(
αk cosk=1

ここで)2πk(i − 1)
+ βk ,n√a2k + b2k ,αk =cos βk = a(k)/αk ,sin βk = −b(k)/αk .

C の場合 i = 0, 1, · · · , n − 1 に対して
r[i] = azero +

kmax
∑

(
αk cos

k=1

ここで

√
αk =

3.7

a2k + b2k ,cos βk = a(k)/αk ,)2πki+ βk ,
nsin βk = −b(k)/αk .

sinti

機能と引数並びFORTRAN の場合
subroutine sinti(n,wsave)integer nreal wsave(*)
C の場合
void sinti(int n, float wsave[])
サブルーチン sinti は sint で使われる配列 wsave を初期化する。 n の素因数分解と三角関数の表が wsave に格納される。

10

入力パラメーターn 変換される配列 r の長さ。その方法は n + 1 が小さな素数の積である時にとても効率的である。
出力パラメーター

wsave 少なくとも int(2.5n + 15) を配置するだけの長さの作業用配列。異なる n の値に対しては、異なる配列 wsave が要求される。wsave の内容は sint の呼び出しまでは変更してはならない。

3.8

sint

機能と引数並びFORTRAN の場合
subroutine sint(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void sint(int n, float x[], float wsave[])
サブルーチン sint は奇関数列 x(i) の離散 Fourier サイン変換を計算する。この変換は後述の出力パラメーター x のところで定義される。
sint は自分自身の非正規化逆変換である。なぜならば、sint 呼び出しに続けてもう一度 sint を呼び出すと、入力数列 x は 2(n + 1) 倍されるからである。
サブルーチン sint で使われる配列 wsave はサブルーチン sinti(n,wsave) 呼び出しで初期化せねばならない。

C プログラマーへのコメント 配列 x[] には、sin に対応する Fourier 係数が格納されるので、添
字は 1 から始まるとした方が自然である。そこで、以下ではfloat x[MAXN+1];のように宣言しておいて、

sint(n, x+1, wsave);のように呼び出すことを仮定して説明をした。こうした場合、FORTRAN と C で、式の定義に違いはなくなる (ただ r(·) を r[·] として読み代えるだけ)。
入力パラメーターn 変換される配列 r の長さ。その方法は n + 1 が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。
wsave sint を呼ぶプログラムにおける少なくとも長さ int(2.5n + 15) の作業用配列。配列 wsave はサブルーチン sinti(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
11

出力パラメーターx i = 1, · · · , n に対して
x(i) =

n∑2x(k) sin

k=1

kiπ
.
n+1

sint 呼び出しに引き続き、もう一度 sint を呼び出すと数列 x は 2(n + 1) 倍される。それゆえ sint は自分自身の非正規化逆変換であると言える。
wsave sint の呼び出しまで壊してはならない、初期化計算の結果を含む。

3.9

costi

機能と引数並びFORTRAN の場合
subroutine costi(n,wsave)integer nreal wsave(*)
C の場合
void costi(int n, float wsave[])
サブルーチン costi は cost で使われる配列 wsave を初期化する。 n の素因数分解と三角関数の表が wsave に格納される。
入力パラメーターn 変換される配列 r の長さ。その方法は n − 1 が小さな素数の積である時にとても効率的である。
出力パラメーター

wsave 少なくとも長さ 3n + 15 の作業用配列。異なる n の値に対しては、異なる配列
wsave が要求される。wsave の内容は cost の呼び出しまでは変更してはならない。

3.10

cost

機能と引数並びFORTRAN の場合
subroutine cost(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void cost(int n, float x[], float wsave[])12

サブルーチン cost は偶関数列 x(i) の離散 Fourier コサイン変換を計算する。この変換は後述の出力パラメーター x のところで定義される。
cost は自分自身の非正規化逆変換である。なぜならば、cost 呼び出しに続けてもう一度 cost を呼び出すと、入力数列 x は 2(n − 1) 倍されるからである。
サブルーチン cost で使われる配列 wsave はサブルーチン costi(n,wsave) 呼び出しで初期化せねばならない。
入力パラメーターn 変換される配列 r の長さ。n は 1 より大きくなければならない。その方法は n − 1 が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。
wsave cost を呼ぶプログラムにおける少なくとも長さ 3n+15 の作業用配列。配列 wsave
はサブルーチン costi(n,wsave) を呼ぶことにより初期化されねばならず、異なる
n の値に対しては、異なる配列 wsave が使われねばならない。この初期化は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーターx 変換された数列。
FORTRAN の場合 i = 1, · · · , n に対して
i−1x(i) = x(1) + (−1)x(n) +

n−1∑

2x(k) cosk=2(k − 1)(i − 1)π.
n−1

cost 呼び出しに引き続き、もう一度 cost を呼び出すと数列 x は 2(n − 1) 倍される。それゆえ cost は自分自身の非正規化逆変換であると言える。
C の場合 i = 0, · · · , n − 1 に対して
x[i] = x[0] + (−1) x[n − 1] +
i

n−2
∑
k=1

2x[k] cos

kiπ
.
n−1

cost 呼び出しに引き続き、もう一度 cost を呼び出すと数列 x は 2(n − 1) 倍される。それゆえ cost は自分自身の非正規化逆変換であると言える。
wsave cost の呼び出しまで壊してはならない、初期化計算の結果を含む。

3.11

sinqi

機能と引数並びFORTRAN の場合
subroutine sinqi(n,wsave)integer nreal wsave(*)
C の場合
void sinqi(int n, float wsave[])
サブルーチン sinqi は sinqf と sinqb の双方で使われる配列 wsave を初期化する。 n の素因数分解と三角関数の表が wsave に格納される。

13

入力パラメーターn 変換する数列の長さ。その方法は n が小さな素数の積である時にとても効率的である。
出力パラメーター

wsave 少なくとも長さ 3n + 15 の作業用配列。 n が変更されない限り、sinqf と sinqbの双方で同じ作業用配列が使える。異なる n の値に対しては、異なる配列 wsaveが使われねばならない。wsave の内容は sinqf あるいは sinqb の呼び出しまでは
変更してはならない。

3.12

sinqf

機能と引数並びFORTRAN の場合
subroutine sinqf(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void sinqf(int n, float x[], float wsave[])
サブルーチン sinqf は quarter wave data の高速 Fourier 変換を計算する。すなわち、sinqf は奇波数を持ったサイン級数表現から数列を計算する。この変換は後述の出力パラメーター x のところで定義される。
sinqb は sinqf の非正規化逆変換である。なぜならば、sinqb 呼び出しに続けて sinqf を呼び出
すと、数列 x は 4n 倍されるからである。
サブルーチン sinqf で使われる配列 wsave はサブルーチン sinqi(n,wsave) 呼び出しで初期化せねばならない。

C プログラマーへのコメント 配列 x[] には、sin に対応する Fourier 係数が格納されるので、添
字は 1 から始まるとした方が自然である。そこで、以下ではfloat x[MAXN+1];のように宣言しておいて、

sinqf(n, x+1, wsave);のように呼び出すことを仮定して説明をした。こうした場合、FORTRAN と C で、式の定義に違いはなくなる (ただ r(·) を r[·] として読み代えるだけ)。
入力パラメーターn 変換される配列 x の長さ。その方法は n が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。

14

wsave sinqf を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン sinqi(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーター

x i = 1, · · · , n に対して
i−1x(i) = (−1)x(n) +

n−1∑

2x(k) sink=1

(2i − 1)kπ2nsinqf 呼び出しに続いて sinqb を呼び出すと、数列 x を 4n 倍することになる。それゆえ sinqb は sinqf の非正規化逆変換である。
wsave sinqf や sinqb の呼び出しまで壊してはならない、初期化計算の結果を含む。

3.13

sinqb

機能と引数並びFORTRAN の場合
subroutine sinqb(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void sinqb(int n, float x[], float wsave[])
サブルーチン sinqb は quarter wave data の高速 Fourier 変換を計算する。すなわち、sinqb は奇波数を持ったサイン級数表現から数列を計算する。この変換は後述の出力パラメーター x のところで定義される。
sinqf は sinqb の非正規化逆変換である。なぜならば、sinqf 呼び出しに続けて sinqb を呼び出
すと、数列 x は 4n 倍されるからである。
サブルーチン sinqb で使われる配列 wsave はサブルーチン sinqi(n,wsave) 呼び出しで初期化せねばならない。

C プログラマーへのコメント 配列 x[] には、sin に対応する Fourier 係数が格納されるので、添
字は 1 から始まるとした方が自然である。そこで、以下ではfloat x[MAXN+1];のように宣言しておいて、

sinqb(n, x+1, wsave);のように呼び出すことを仮定して説明をした。こうした場合、FORTRAN と C で、式の定義に違いはなくなる (ただ r(·) を r[·] として読み代えるだけ)。

15

入力パラメーターn 変換される配列 x の長さ。その方法は n が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。
wsave sinqb を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン sinqi(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーター

x i = 1, · · · , n に対して
x(i) =

n∑4x(k) sink=1

(2k − 1)iπ.
2nsinqb 呼び出しに続いて sinqf を呼び出すと、数列 x を 4n 倍することになる。それゆえ sinqf は sinqb の非正規化逆変換である。
wsave sinqf や sinqb の呼び出しまで壊してはならない、初期化計算の結果を含む。

3.14

cosqi

機能と引数並びFORTRAN の場合
subroutine cosqi(n,wsave)integer nreal wsave(*)
C の場合
void cosqi(int n, float wsave[])
サブルーチン cosqi は cosqf と cosqb の双方で使われる配列 wsave を初期化する。 n の素因数分解と三角関数の表が wsave に格納される。
入力パラメーター

n 変換する配列の長さ。その方法は n が小さな素数の積である時にとても効率的である。
出力パラメーター

wsave 少なくとも長さ 3n + 15 の作業用配列。n が変更されない限り、cosqf と cosqbの双方で同じ作業用配列が使える。異なる n の値に対しては、異なる配列 wsaveが使われねばならない。wsave の内容は cosqf あるいは cosqb の呼び出しまでは
変更してはならない。

16

3.15

cosqf

機能と引数並びFORTRAN の場合
subroutine cosqf(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void cosqf(int n, float x[], float wsave[])
サブルーチン cosqf は quarter wave data の高速 Fourier 変換を計算する。すなわち、cosqb は奇波数を持ったコサイン級数表現から係数を計算する。この変換は後述の出力パラメーター x のところで定義される。
cosqf は cosqb の非正規化逆変換である。なぜならば、cosqf 呼び出しに続けて cosqb を呼び出すと、数列 x は 4n 倍されるからである。
サブルーチン cosqf で使われる配列 wsave はサブルーチン cosqi(n,wsave) 呼び出しで初期化せねばならない。
入力パラメーターn 変換される配列 x の長さ。その方法は n が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。
wsave cosqf を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン bf cosqi(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーター

x Fourier 係数を格納する。
FORTRAN の場合 i = 1, 2, · · · , n に対して
x(i) = x(1) +n∑2x(k) cosk=2(2i − 1)(k − 1)π.
2ncosqf 呼び出しに引き続き cosqb を呼び出すと数列 x は 4n 倍される。それゆえ cosqb は cosqf の非正規化逆変換である。
C の場合 i = 0, 1, · · · , n − 1 に対して
x[i] = x[0] +

n−1
∑
k=1

2x[k] cos

(2i + 1)kπ.
2ncosqf 呼び出しに引き続き cosqb を呼び出すと数列 x は 4n 倍される。それゆえ cosqb は cosqf の非正規化逆変換である。
wsave cosqf や cosqb の呼び出しまで壊してはならない、初期化計算の結果を含む。
17

3.16

cosqb

機能と引数並びFORTRAN の場合
subroutine cosqb(n,x,wsave)integer nreal x(*),wsave(*)
C の場合
void cosqb(int n, float x[], float wsave[])
サブルーチン cosqb は quarter wave data の高速 Fourier 変換を計算する。すなわち、cosqb は奇波数を持ったコサイン級数表現から数列を計算する。この変換は後述の出力パラメーター x のところで定義される。
cosqb は cosqf の非正規化逆変換である。なぜならば cosqb 呼び出しに続けて cosqf を呼び出
すと数列 x を 4n 倍することになるから。
サブルーチン cosqb で使われる配列 wsave はサブルーチン cosqi(n,wsave) の呼び出しで初期化せねばならない。
入力パラメーターn 変換される配列 x の長さ。その方法は n が小さな素数の積である時にとても効率的である。
x 変換される数列を含む長さ n の配列。
wsave cosqb を呼ぶプログラムにおける少なくとも長さ 3n + 15 の作業用配列。配列
wsave はサブルーチン cosqi(n,wsave) を呼ぶことにより初期化されねばならず、異なる n の値に対しては、異なる配列 wsave が使われねばならない。この初期化
は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーターx 関数値の列が格納される。
FORTRAN の場合 i = 1, · · · , n に対して
x(i) =

n∑4x(k) cosk=1

(2k − 1)(i − 1)π.
2ncosqb 呼び出しに続けて cosqf を呼び出すと数列 x を 4n 倍することになる。
それゆえ cosqf は cosqb の非正規化逆変換である。
C の場合 i = 0, 1, · · · , n − 1 に対して
x[i] =

n−1∑

4x[k] cos

k=0

(2k + 1)iπ.
2ncosqb 呼び出しに続けて cosqf を呼び出すと数列 x を 4n 倍することになる。
それゆえ cosqf は cosqb の非正規化逆変換である。
wsave cosqf や cosqb の呼び出しまで壊してはならない、初期化計算の結果を含む。
18

3.17

cﬀti

機能と引数並びFORTRAN の場合
subroutine cﬀti(n,wsave)integer nreal wsave(*)
C の場合
void cﬀti(int n, float wsave[])
サブルーチン cﬀti は cﬀtf や cﬀtb の双方で使われる配列 wsave を初期化する。n の素因数分解と三角関数の表が wsave に格納される。
入力パラメーターn 変換される数列の長さ。
出力パラメーター

wsave 少なくとも長さ 4n + 15 の作業用配列。n が変更されない限り、cﬀtf と cﬀtb の双方で同じ作業用配列が使える。異なる n の値に対しては、異なる配列 wsave が使われねばならない。wsave の内容は cﬀtf あるいは cﬀtb の呼び出しまでは変更
してはならない。

3.18

cﬀtf

機能と引数並びFORTRAN の場合
subroutine cﬀtf(n,c,wsave)integer nreal wsave(*)
complex c(*)
C の場合
void cﬀtf(int n, COMPLEX c[], float wsave[])
サブルーチン cﬀtf は複素離散 Fourier 変換を計算する (Fourier analysis)、すなわち cﬀtf 複素周期数列の Fourier 係数を計算する。この変換は後述の出力パラメーター c のところで定義される。
この変換は正規化されていない。正規化した変換を得る為には出力を n で割らなければならない。
そうしないで、cﬀtf の呼び出しに続けて cﬀtb の呼び出しをすると、数列を n 倍することになる。
サブルーチン cﬀtf で使われる配列 wsave はサブルーチン cﬀti(n,wsave) の呼び出しによって初期化せねばならない。

19

入力パラメーターn 複素数列 c の長さ。その方法は n が小さな素数の積である時にとても効率的である。
c 数列を含む長さ n の複素数配列。
wsave cﬀtf を呼ぶプログラムにおける少なくとも長さ 4n+15 の作業用配列。配列 wsave
はサブルーチン cﬀti(n,wsave) を呼ぶことにより初期化されねばならず、異なる
n の値に対しては、異なる配列 wsave が使われねばならない。この初期化は n が変更されない限り、繰り返す必要はないので、引き続く変換は最初よりも速く得られる。
出力パラメーターc Fourier 係数が格納される。
FORTRAN の場合 j = 1, · · · , n に対して
c(j) =

n∑c(k) exp

k=1

ここで i =√

−2πi(j − 1)(k − 1)
,
n

−1.

C の場合 j = 0, 1, · · · , n − 1 に対して
c[j] =n−1∑
k=0

ここで i =√

c[k] exp−2πijk,n−1.

wsave サブルーチン cﬀtf もしくは cﬀtb の呼び出しまでの間に破壊してはならない初期化の計算結果を含む。

3.19

cﬀtb

機能と引数並びFORTRAN の場合
subroutine cﬀtb(n,c,wsave)integer nreal wsave(*)
complex c(*)
C の場合
void cﬀtb(int n, COMPLEX c[], float wsave[])
サブルーチン cﬀtb は複素離散 Fourier 逆変換を計算する (Fourier 同調)、すなわち cﬀtb 複素周期数列をその Fourier 係数から計算する。この変換は後述の出力パラメーター c のところで定義される。
cﬀtf の呼び出しに続けて cﬀtb の呼び出しをすると、数列を n 倍することになる。
サブルーチン cﬀtb で使われる配列 wsave はサブルーチン cﬀti(n,wsave) の呼び出しによって初期化せねばならない。

20

入力パラメーターn 複素数列 c の長さ。その方法は n が小さな素数の積である時にとても効率的である。
c 数列を含む長さ n の複素数配列。
wsave cﬀtb を呼ぶプログラムにおける少なくとも長さ 4n + 15 の実作業用配列。配列
wsave はサブルーチン cﬀti(n,wsave) を呼び出すことによって初期化しなければならず、異なる n の値に対しては異なる配列 wsave を使わねばならない。cﬀtf とcﬀtb で同じ配列 wsave が使える。
出力パラメーターc 関数値の列が納められる。
FORTRAN の場合 j = 1, · · · , n に対して
c(j) =

n∑c(k) exp

k=1

ここで i =√

i(j − 1)(k − 1)2π
,
n

−1.

C の場合 j = 0, 1, · · · , n − 1 に対して
c[j] =n−1∑
k=0

ここで i =√

c[k] exp

2πijk
,
n

−1.

wsave サブルーチン cﬀtf もしくは cﬀtb の呼び出しまでの間に破壊してはならない初期化の計算結果を含む。
[“send index for vﬀtpk” describes a vectorized version of ﬀtpack]

21

付録

A

dﬀtpack の C 言語プログラムからの利用

A.1

筆者が利用している環境の説明

Apple の Mac (OS は 10.11) に、Apple が提供している Xcode と、MacPorts の GCC 5 (gfortranを含む) をインストールしている。

A.2

dﬀtpack の入手とインストール

http://www.netlib.org/fftpack/dp.tgz を入手して展開する。以下では curl というコマンド
でダウンロードしている (使っている Mac OS に、curl が標準で備わっている)。
curl -O http://www.netlib.org/fftpack/dp.tgz
tar xzf dp.tgz
cd dfftpack

ここで Makefile を適当に編集する。筆者の場合は、FORTRAN コンパイラーを指定するマクロ FC
を適当なものに置き換えるくらいで済んだ。例えば gfortran-mp-5 とするには、
Makefile の書き換え

#FC=g77
FC=gfortran-mp-5


のように書き換える (FC=g77 というのをコメント・アウトして、FC=gfortran-mp-5 という行を書き加えた)。後は

make
make test

これで各サブルーチンのテストの結果、誤差が小さいことを確認する。問題がなければインストー
ルする。




sudo make install

これで /usr/local/lib/libdfftpack.a, /usr/local/include/dfftpack.h がインストールされる。

オリジナルの単精度版をインストールするにはオリジナルの ﬀtpack の入手・インストールも大体同様であるが、dp.tar.gz のようなアーカイ
ブ・ファイルがないので、http://www.netlib.org/fftpack/ にあるファイルを全部入手することになる。WWW ブラウザーでやるのは面倒なので、私は wget5 を用いてwget --mirror --no-parent -P . -nH http://www.netlib.org/fftpack/
cd fftpackとした (fftpack/ の / は重要)。この後は上と同様である。
5

これは Mac OS には標準で入っていないが、色々便利なので導入しても良いだろう。Mac Ports が使えるならば、
sudo port install wget とすればインストールできる。

22

A.3

dﬀtpack の C 言語プログラムからの利用

UNIX の伝統で、C コンパイラーと FOTRAN コンパイラーのオブジェクト・ファイルに互換性があるので、簡単に C 言語のプログラムから、FORTRAN でコンパイルして作った FFTPACK のラ
イブラリィーが呼び出せる。
大したものではないけれど、http://nalab.mind.meiji.ac.jp/~mk/fftpack-c-interface.tar.
gz を公開する。それを使うと以下のようにしてインストール出来る (dﬀtpack の入手から通して記述してある)。


curl -O http://www.netlib.org/fftpack/dp.tgz
tar xzf dp.tgz
cd dfftpackcurl -O http://nalab.mind.meiji.ac.jp/~mk/fftpack-c-interface.tar.gz
tar xzf fftpack-c-interface.tar.gz
mv fftpack-c-interface/* .
mv Makefile Makefile.org
cp Makefile-dfftpack Makefile

こうしてから Makefile を編集する。書き換え方は上と同じである。その後も上と同じ。







make
make test
sudo make install




A.4

まとめ


倍精度版



curl -O http://www.netlib.org/fftpack/dp.tgz
tar xzf dp.tgz
cd dfftpackcurl -O http://nalab.mind.meiji.ac.jp/~mk/fftpack-c-interface.tar.gz
tar xzf fftpack-c-interface.tar.gz
mv fftpack-c-interface/* .
mv Makefile Makefile.org
cp Makefile-dfftpack Makefile
makemake test
sudo make install

これで /usr/local/lib/libdfftpack.a, /usr/local/include/dfftpack.h がインストールされ
る。

23


単精度版



wget --mirror --no-parent -P . -nH http://www.netlib.org/fftpack/
cd fftpackcurl -O http://nalab.mind.meiji.ac.jp/~mk/fftpack-c-interface.tar.gz
tar xzf fftpack-c-interface.tar.gz
mv fftpack-c-interface/* .
mv Makefile Makefile.org
cp Makefile-fftpack Makefile
makemake test
sudo make install

これで /usr/local/lib/libfftpack.a, /usr/local/include/fftpack.h がインストールされる。

A.5

おまけ: サンプル・プログラム(工事中、まだバグがあると思う。)
スペクトル法というと大げさだけれど、熱方程式の初期値境界値問題を FFT を用いて解くプログ
ラムを紹介する。

A.5.1

周期境界条件の場合

まずは周期境界条件の場合から。

(2)

∂u
∂ 2u
(x, t) ((x, t) ∈ (0, 2π) × (0, ∞)),
(x, t) =∂t∂x2
u(0, t) = u(2π, t) = 0, ux (0, t) = ux (2π, t) (t ∈ (0, ∞)),

(3)u(x, 0) = f (x) (x ∈ [0, 2π]).

(1)

解 u が x について周期 2π であることは、比較的簡単に分かる。それゆえ、各時刻 t において、u(·, t) は Fourier 級数展開できる。すなわち

u(x, t) =∞a0 (t) ∑+(ak (t) cos kx + bk (t) sin kx)
2
k=1の形に展開できる。この ak , bk については、

a′k (t) = −k 2 ak (t),ak (0) = Ak ,b′k (t) = −k 2 bk (t),bk (0) = Bkが成り立つことが分かる。ただし、Ak , Bk は f の Fourier 係数である:
∫
∫
1 2π
1 2π
Ak =
f (x) cos kx dx (k = 0, 1, · · · ), Bk =f (x) sin kx dx
π 0
π 0(k = 1, 2, · · · ).

そこで関数 f の値を計算してから、離散 Fourier 変換して Ak , Bk (の近似値) を求め、

ak (t) = Ak e−k t ,2bk (t) = Bk e−k t2を計算してから、離散共役 Fourier 変換して u(x, t) (の近似値) を計算する、というアルゴリズムが考えられる。
FFTPACK は、実数値周期関数を扱うために、

24

(a) dzffti(), dzfftf(), dzfftb()(b) dffti(), dfftf(), dfftb()という 2 系統の関数群を用意してある。
まず (a) から。

アルゴリズム 1 と dzffti(), dzfftf(), dzfftb() の利用



double x[N+1], u[N+1], a[N/2+1], b[N/2+1], at[N/2+1], bt[N/2+1];(中略)
h=2*pi/N;// t=0 での ufor (i = 0; i <= N; i++) {x[i] = i * h;u[i] = f(x[i]);}// t=0 での u の Fourier 係数を離散 Fourier 変換で求める
dzffti(N, work);
dzfftf(N, u, a, a+1, b+1, work);// if (N % 2 == 0) b[N/2] = 0.0;
(中略)// 時刻 t における ak(t), bk(t) を求める
at[0]=a[0];for (k = 1; k <= N/2; k++) {factor = exp(- k * k * t);at[k] = a[k] * factor;
bt[k] = b[k] * factor;}// 時刻 t での u を離散共役 Fourier 変換で求める
dzfftb(N, u, at, at+1, bt+1, work);// 例えばグラフでも描いてみる
u[N] = u[0];g_move(x[0], u[0]);for (i = 1; i <= N; i++)g_plot(x[i], u[i]);

N の値が int 型の変数 N に記憶されているとする。
f は周期関数なので、u[0]= f (x0 ) と u[N]= f (xN ) の値は等しい。
dzfftf() に u[i]= f (xi ) (i = 0, 1, · · · , N − 1) という N 個の数値を入力して、離散 Fourier 係
数 A0 , A1 , B1 , · · · (N 個) を求めている。詳しく述べると、• N が偶数のとき、dzfftf() から A0 , A1 , B1 , . . . , AN/2−1 , BN/2−1 , AN/2 という N 個の数が返される。(このとき実は b[N/2] に値は設定されないので、k= N/2 に対して、bt[k] = b[k]* factor; を実行するのは厳密にはおかしい。事前に b[N/2]=0.0; のようにしておくのが無難かもしれない。)• N が奇数のとき、N/2 という式の値は (N − 1)/2. dzfftf() から A0 , A1 , B1 , . . . , A(N −1)/2 ,B(N −1)/2 という N 個の値が返される。
dzfftb() に a0 (t), a1 (t), b1 (t), · · · (N 個) の数値を入力して、u[i]= u(xi , t) (i = 0, 1, · · · , N − 1)という N 個の数値を求めている。
上のプログラムでは、b[0], bt[0] は未使用である。節約するようにプログラムを書くことも出来るが、この方が分かりやすいであろう。
次に (b) について述べる。数学的には (a) と等価であることに注意しよう。

25


アルゴリズム 1 と dffti(), dfftf(), dfftb() の利用



double x[N+1], r[N+1], rt[N+1];(中略)
h=2*pi/N;// t=0 での ufor (i = 0; i < N; i++) {x[i] = i * h;r[i] = f(x[i]);}// 時刻 t=0 での Fourier 係数を離散 Fourier 変換で求める
dffti(N, work);
dfftf(N, r, work);// if (N % 2 == 0) r[N] = 0.0;
for (k = 0; k < N; k++) r[k] /= N;(中略)// 時刻 t における ak(t), bk(t) を求める
rt[0]=r[0];for (i = 1; i < N; i += 2) {k = (i + 1) / 2;
factor = exp(- k * k * t);rt[i]= r[i]* factor;
rt[i+1] = r[i+1] * factor;}// 時刻 t での u を離散共役 Fourier 変換で求める
dfftb(N, rt, work);// 例えばグラフでも描いてみる
rt[N] = rt[0];g_move(x[0], rt[0]);for (i = 1; i <= N; i++)g_plot(x[i], rt[i]);f は周期関数なので、r[0]= f (x0 ) と r[N]= f (xN ) の値は等しい。
dfftf() に r[i]= f (xi ) (i = 0, 1, · · · , N − 1) という N 個の数値を入力して、離散 Fourier 係数
A0 , A1 , B1 , · · · (N 個) を求め、r[k] (k = 0, 1, · · · , N − 1) に返される。詳しく述べると、• N が偶数のとき、dfftf() から A0 , A1 , B1 , . . . , AN/2−1 , BN/2−1 , AN/2 という N 個の数が返される。(r[N] に値が設定されないが、rt[N] に値を設定するのは厳密にはおかしいので、
r[N]=0.0; のようにしておくのが無難かもしれない。)• N が奇数のとき、N/2 という式の値は (N − 1)/2. dfftf() から A0 , A1 , B1 , . . . , A(N −1)/2 ,B(N −1)/2 という N 個の値が返される。
dzfftb() に a0 (t), a1 (t), b1 (t), · · · (N 個) の数値を入力して、rt[i]= u(xi , t) (i = 0, 1, · · · , N −1)という N 個の数値を求めている。
最後に、複素数バージョンの zffti(), zfftf(), zfftb() について述べる。

26


アルゴリズム 1 と zffti(), zfftf(), zfftb() の利用#include <complex.h>double x[N+1];
complex double c[N+1], w[N+1];(中略)
h=2*pi/N;// t=0 での ufor (i = 0; i <= N; i++) {x[i] = i * h;c[i] = f(x[i]);}// 時刻 t=0 での Fourier 係数を離散 Fourier 変換で求める
zffti(N, work);
zfftf(N, c, work);for (k = 0; k < N; k++) c[k] /= N;(中略)// 時刻 t における ck(t) を求めて w[k] にセット
w[0] = c[0];for (k = 1; k <= N / 2; k++) {factor = exp(- k * k * t);w[k]= c[k]* factor;w[N - k] = c[N - k] * factor;}// 時刻 t での u を離散共役 Fourier 変換で求める
zfftb(N, w, work);// 例えばグラフでも描いてみる (Re w=u である)
w[N] = w[0];g_move(x[0], creal(w[0]));for (i = 1; i <= N; i++)g_plot(x[i], creal(w[i]));

f は周期関数なので、最初の c について、c[0]= f (x0 ) と c[N]= f (xN ) の値は等しい。
zfftf() に c[i]= f (xi ) (i = 0, 1, · · · , N − 1) という N 個の数値を入力し、離散 Fourier 係数−1{Ck }N
k=0 を求めている。
2
zfftb() に w[k]= Ck e−ℓ t (k = 0, 1, · · · , N − 1) を入力して (ただし、ℓ = min{k, N − k})、w[i]= u(xi , t) (i = 0, 1, · · · , N − 1) を求めている。
−1
f が実数値関数であるから、Fourier 係数 {cn }n∈Z は c−n = cn を満たし、離散 Fourier 係数 {Ck }N
k=0⌊ N −1 ⌋は、C0 ∈ R, CN −k = Ck (k = 1, 2, · · · , 2 ) を満たす。つまり、この問題を解くために zffti(),zfftf(), zfftb() を使うのは無駄なことをしている。
• N が偶数のとき、k = N/2 に対して、c[k] と c[N − k] は同じものを指すことに注意せよ。
• N が奇数のとき、N/2 の値は N2−1 であり、k = N2−1 のとき、N − k = N2+1 である。
A.5.2

Dirichlet 境界条件

(5)

∂ 2u
∂u
(x, t) =(x, t) ((x, t) ∈ (0, π) × (0, ∞)),∂t∂x2
u(0, t) = u(π, t) = 0 (t ∈ (0, ∞)),

(6)u(x, 0) = f (x) (x ∈ [0, π]).

(4)

27

各時刻 t で関数 u(·, t) を Fourier 正弦展開する。すなわち

u(x, t) =∞∑bk (t) sin kx (x ∈ [0, π])k=1の形に展開できる。

b′k (t) = −k 2 bk (t),bk (0) = Bkが成り立つ。ただし Bk は f の Fourier 係数である。
∫
2 π
f (x) sin kx dx (k = 1, 2, · · · ).
Bk =π 0
そこで関数 f の xi (i = 1, 2, . . . , N − 1) での値を計算してから、離散 Fourier 正弦変換して Bk (の近似値) を求めて
2
bk (t) = Bk e−k t (k = 1, 2, · · · , N − 1)を計算してから、離散 Fourier 正弦変換して u(xi , t) (i = 1, 2, · · · , N − 1) の近似値を求める、というアルゴリズムが考えられる。

アルゴリズム 2 と dsinti(), dsint() の利用



double x[N+1], b[N+1], r[N+1];(中略)
h=2*pi/N;// t=0 での ufor (i = 1; i < N; i++) {x[i] = i * h;b[i] = f(x[i]);}x[0] = x[N] = 0;// 時刻 t=0 での Fourier 係数を離散 Fourier 変換で求める
dsinti(N - 1, work);
dsint(N - 1, b + 1, work);for (k = 1; k < N; k++) b[k] /= (2*N);(中略)// 時刻 t における bk(t) を求める
for (k = 1; k < N; k++)r[k] = b[k] * exp(- k * k * t);// 離散共役 Fourier 変換
dsint(N - 1, r + 1, work);// 例えばグラフでも描いてみる
r[0] = r[N] = 0.0;
g_move(x[0], r[0]);for (i = 1; i <= N; k++)g_plot(x[i], r[i]);dsint() に、b[i]= f (xi ) (i = 1, 2, · · · , N − 1) という N − 1 個の数値を入力し、離散 Fourier 正−1
弦係数 {Bk }N
k=1 を求めている。
2
dsint() に、r[k]= Bk e−k t (k = 1, 2, . . . , N − 1) という N − 1 個の数値を入力し、r[i]= u(xi , t)(i = 1, 2, · · · , N − 1) を求めている。

28

A.5.3

Neumann 境界条件

(8)

∂ 2u
∂u
(x, t) =(x, t) ((x, t) ∈ (0, π) × (0, ∞)),∂t∂x2
ux (0, t) = ux (π, t) = 0 (t ∈ (0, ∞)),

(9)u(x, 0) = f (x) (x ∈ [0, π]).

(7)∞a0 (t) ∑u(x, t) =+ak (t) cos kx
2
k=1の形に展開できる。

a′k (t) = −k 2 ak (t),ak (0) = Akが成り立つ。ただし Ak は f の Fourier 係数である。
∫
2 π
f (x) cos kx dx
Ak =
π 0

(k = 0, 1, · · · ).

そこで関数 f の xi (i = 0, 1, 2, . . . , N ) での値を計算してから、離散 Fourier 余弦変換して Ak (の近似値) を求め2
ak (t) = Ak e−k t (k = 0, 1, · · · , N )を計算してから、離散 Fourier 余弦変換して u(xi , t) (i = 0, 1, · · · , N ) の近似値を求める、というア
ルゴリズムが考えられる。

アルゴリズム 3 と dcosti(), dcost() の利用



double x[N+1], a[N+1], r[N+1];(中略)
h=2*pi/N;// t=0 での ufor (i = 0; i <= N; i++) {x[i] = i * h;a[i] = f(x[i]);}// 時刻 t=0 での Fourier 係数を離散 Fourier 変換で求める
dcosti(N + 1, work);
dcost(N + 1, a, work);for (k = 0; k <= N; k++) a[k] /= (2*N);(中略)// 時刻 t における ak(t) を求める
r[0] = a[0];for (k = 1; k <= N; k++)r[k] = a[k] * exp(- k * k * t);// 離散共役 Fourier 変換
dcost(N + 1, r, work);// 例えばグラフでも描いてみる
g_move(x[0], r[0]);for (i = 1; i <= N; i++)g_plot(x[i], r[i]);

dcost() に、a[i]= f (xi ) (i = 0, 1, · · · , N ) という N + 1 個の数値を入力し、離散 Fourier 余弦係数 a[k]= Ak (k = 0, 1, · · · , N ) を求めている。
2
dcost() に、r[k]= Ak e−k t (k = 0, 1, . . . , N ) という N + 1 個の数値を入力し、r[i]= u(xi , t)(i = 0, 1, · · · , N ) を求めている。

29