sum演算子

sum演算子について
繰り返しループとして使う

sum演算子について

gnuplot4.5のある時期以降、和を取るためのsum演算子が導入されました。この演算子は以下のように用います。

sum [i = 1:10] i

ここで、sumのすぐ後の大かっこ内にあるiはダミー変数で、その後の数字の範囲（上の例では1から10の間）を1ずつ変化しながら動きます。そして、それぞれのダミー変数の値に対して、大かっこの後の表現が計算され、最後にその和が返されます。つまり上のコマンドでは以下の数式を計算していることになります。

例えば、以下のようにすると、確かに1から10までの和である55が計算できていることが分かります。

print sum [i = 1:10] i

sum演算子は関数定義の中で使ったり、plotのusingの中で使ったりすることもできます。まず、関数定義の中で使う例として、級数展開の項で説明した矩形波のFourier展開

を示します。1行目でsum演算子を使って関数step_fourier(x,n)を定義しています。これは、上のFourier級数のnまでの部分和を返す関数です。このスクリプトを実行すると、下の図のようになります。なお、再帰定義を使った級数展開では和をとる繰り返し回数に制限があり、n=100程度が上限でしたが、こちらの例ではn=1000でも（やや時間はかかりますが）問題なく計算できている点にも注目して下さい。

step_fourier(x, n) = sum [i=1:n] 2.0/pi*(1-(-1)**i)*sin(i*x)/i
set samples 5000
plot step_fourier(x, 10), step_fourier(x, 100), step_fourier(x, 1000)

また、usingの中でsum演算子を使う例として、テストの平均点のプロットを紹介します。ある授業の小テストの得点が回ごとに並んだ以下のようなデータファイルtest_result.datがあったとします。

#A君 B君 C君 … Z君
 50  35  59 … 100
 80  20  49 …  99
 90  30  29 …  87

このとき、以下のようなスクリプトを用いれば、各回の平均値をプロットすることができます。

N = 26 # 生徒の総数
plot "test_result.dat" using 0:(sum [i=1:N] column(i))*(1.0/N) w lp

このusingの中身はusing 0:($1+$2+$3+...+$26)/26.0と同じなのですがsumをつかうとよりすっきり分かりやすく記述できます。また、生徒数が変わった場合などの変更も容易です。

他に、数値積分も和で表せますので、sum演算子を使って実現できます。

繰り返しループとして使う

sumを累次評価演算子(,)や三項演算子と組み合わせれば、和以外の一般的な繰り返しループとしても使うことができます。

そのための基本的なやり方は以下の通りです。

sum [i=M:N] (何かやりたいこと… , (i==N ? 評価したい値 : 0))

つまり、累次評価演算子を用いて、繰り返し処理をしたいことを行っていき、累次評価演算子の最後に和を取るべき値をsumに与えるのですが、ここでループの最後（つまりi==N）以外は0を与えてやるわけです。すると、sumの結果は、ループの最後に評価した値に一致します。

まずは、簡単な例を以下に示します。

a = 1.0
print sum [i=1:10] (a=a*i, (i==10 ? a : 0))

この例では、iが1から10まで1ずつ増えていく間、変数aの値をi倍していきます。しかし、sumで和を取る値は、iが10以外ならば常にゼロで、iが10のときだけaの値になります。結局、sumが返すのはi==10のときのaの値、つまり10! = 3628800となります。

次に、もう少し複雑な例をいくつか紹介します。

まず、正弦波の無限積展開

の例を示します。

f(x,n) = x* sum [i=1:n] ( (i==1 ? a=(1.0-(x/pi/i)**2) : a=a*(1.0-(x/pi/i)**2)) , (i==n ? a : 0.0) ) 
set samples 1000
plot [-50:50][-5:5] f(x,10), f(x,100), f(x,1000), sin(x)

スクリプトの実行例からわかるように、上の式の積の上限を大きくしていくと、確かに正弦関数に近付いていくことが分かります。

次に、ある整数n以下の素数の数をカウントする関数の例を紹介します。

# int(n)が素数の場合は1を返し、そうでない場合は0を返す関数
prime(n) = sum [i=1:int(sqrt(n))] \
    ( (i==1 ? _result = 0 : _result = _result + (int(n)%i==0 ? 1 : 0)) , \
    ( i==int(sqrt(n)) && _result < 1 ? 1 : 0 ))

# n以下の素数の数を返す関数
primes(n) = sum [i=2:int(n)] prime(i)

set xrange [0:1000]
set samples 1001
plot primes(x)

上記のprime(n)関数はいわゆる「試し割り法」を用いて、nが素数かどうか判定します。この関数は引数nが素数の場合は1、そうでない場合は0を返すので、単純にprime(1)+prime(2)+....+prime(n)という和を取ってやることで、n以下の素数の数を返す関数primes(n)を定義できます。

上記スクリプトを実行すると、以下のようになります。例えば、グラフの右端では、ちゃんと1000以下の素数の数である168を返しています。

なお、このようなある整数n以下の素数の数を表す関数はπ(n)と書かれ、「素数定理」によると十分大きいnに対して以下の近似式が成り立ちます。

これを確かめるために、以下のようなデータファイルintegers.datを準備します。

そして、以下のようなスクリプトを走らせます。

# int(n)が素数の場合は1を返し、そうでない場合は0を返す関数
prime(n) = sum [i=1:int(sqrt(n))] ( (i==1 ? _result = 0 : _result = _result + (int(n)%i==0 ? 1 : 0)) , ( i==int(sqrt(n)) && _result < 1 ? 1 : 0 ))

# n以下の素数の数を返す関数
primes(n) = sum [i=2:int(n)] prime(i)

set logscale xy
plot "integers.dat" using 1:(primes($1)), x/log(x)

実行にはかなりの時間がかかります（私の環境では15分くらい?）が、結果は以下のようになります。確かに、素数定理が成り立っていることが分かります。