中学・高校数学でオイラーのγの近似値を求めるプログラミング：数学×Pythonプログラミング入門（1/2 ページ）

調和数列を使ってオイラーのγ（ガンマ）の近似値を求める問題を通して、Σを使った総和の計算をプログラミングする。また、Pythonのリストや繰り返し処理についてまとめ、身近な事例を使って、繰り返し処理の制御変数とリストのインデックスをどう対応させるかという問題について考える。

» 2021年08月19日 05時00分公開

[羽山博，著]

この記事は会員限定です。会員登録（無料）すると全てご覧いただけます。

連載目次

　前回は、ウオーミングアップの回として、フェルマーの小定理やアンドリカの予想、超球の体積といったトピックを通してPythonプログラミングの初歩を振り返りました。

　今回から、数式などの「数学の言葉」で表された問題をプログラムとして記述していくことに本格的に取り組んでいきましょう。

　まず、調和数列を使ってオイラーのγ（ガンマ）と呼ばれる値の近似値を求める問題を例に、Σ（総和）の計算をプログラミングするために使われるお決まりのパターンを見ていきます。さらに、繰り返し処理を行うときに最もつまずきやすいポイント――繰り返し処理の制御変数とリストのインデックスをどう対応させるかという問題――について、可処分所得の対前年度比の計算などの身近な例を使って考えていきます。

　練習問題としては、ネイピア数（自然対数の底e）の近似値を求める例、2つの変数の共分散を求める例、株価のn日移動平均を求める例を取り上げます。

連載：

『数学×Pythonプログラミング入門 ― 中学・高校数学で学ぶ』

この連載では、中学や高校で学んだ数学を題材にして、Pythonによるプログラミングを学びます。といっても、数学の教科書に載っている定理や公式だけに限らず、興味深い数式の例やAI／機械学習の基本となる例を取り上げながら、数学的な考え方を背景としてプログラミングを学ぶお話にしていこうと思います。

羽山博

筆者紹介： IT系ライター、大学教員（非常勤）。書道、絵画を経て、ピアノとバイオリンを独学で始めるも学習曲線は常に平坦。趣味の献血は、最近脈拍が多く98回で一旦中断。

目標：オイラーのγの近似値を求める

　オイラーという数学者はあまりにも有名なので、そんな名前は聞いたことがない、という人はほとんどいないのではないでしょうか。世界で最も美しいと言われるオイラーの関係式（e^iπ＝－1）にも登場しますし、身近なところでは「一筆書きできる道」はオイラー閉路と呼ばれます。しかし、オイラーのγ（ガンマ）（またはオイラーの定数）と呼ばれる値はあまり知られていないかもしれません。

　オイラーのγは以下のような値ですが、実は、この値の性質はまだよく分かっていません。しかも、無理数なのかどうかもいまだに分かっていません（ちなみに、前回取り上げたガンマ関数のx＝1における微分係数Γ'(1) と－γが等しくなります。数学は不思議だらけですね）。

　この値を求めるための数式は以下の通りで、調和数列と自然対数が使われています。

　オイラーのγは、AI／機械学習ですぐに使うような値ではありませんが、数式をプログラミングしていくための基本的なパターンを含んでいるので、このプログラミングに取り組むことは、AI／機械学習だけでなく、どのような計算をする場合にも、そのための基礎を養うのに役立ちます。

　さて、この数式をプログラムとして表すことを考えていくわけですが、Σが総和を表すことが分かっても、式の中に書かれている∞（無限大）の扱いに困りますね。プログラムで無限に足し算を行うことはできないからです（永遠に答えが求められませんから）。そこで、近似値を求めるためにmに好きな値（ある程度大きな値）を指定するものと考えてください。

　では、コードを書いてみてください。……と言いたいところですが、いきなりそう言われてスラスラとできる人は多くはないかと思います（いや、これぐらいなら楽勝だよ、と思った方はぜひ独力で取り組んでみてください。答え合わせはこちら）。もしかすると「ホントに中学／高校数学の範囲なの？」と思われる方もいるかもしれませんね。しかし、間違いなく中学／高校の知識だけで理解できる数式です。そこで、この数式を理解するために、以下の内容についておさらいしておきましょう。（1）式を見て、ちょっと「ビビった」という人もご心配なく。一歩ずつゆっくり進めるので必ず理解できます。

数列 …… 等差数列や等比数列
級数 …… 数列の和の公式
自然対数 …… 底がe（＝2.71828...）の対数

　なお、これらの内容について、既に十分に理解しているのでおさらいなんて必要ないという方は、2章まで読み飛ばしてもらって構いません。では、数列から見ていきましょう。

1. 数列のおさらいをかねたサンプルプログラム

　数列とは、a₁,a₂,...,a_nという数の並びです。これを簡単に{a_n}と表すこともあります。

（例）

1,4,3,6,1 …… サイコロを5回振って出た目の並び
0.11, －0.7, －0.05, －0.49, 0.26 …… 2021年8月5日～11日のAppleの株価（終値）の前日比（単位：ドル）

　これらの例では、値がどう並んでいるかというルールはよく分かりませんが、数が並んでいるわけですから、レッキとした数列です*1。ただし、高校数学では、等差数列や等比数列などのルールが決まったものだけを数列と呼んでいるようです。それらのルールは以下の通りです。

等差数列： 初項a₁に一定の値（公差d）を足して次々に各項が求められる
等比数列： 初項a₁に一定の値（公比r）を掛けて次々に各項が求められる

*1　サイコロの目の並びや株価のデータなどは、数列としてではなく、ベクトルとして表すことの方が一般的です。例えば株価が上下するパターンがどのような要因と連動しているかを調べたい場合などにベクトルを使った計算が便利だからです。練習問題で取り扱っている気温やビールの売り上げなども数列と言えますが、ベクトルとして扱うことにより、さまざまな分析ができます。

　等差数列と等比数列の例を見ておきましょう。

1, 3, 5, 7, 9 …… 初項1、公差2の等差数列（奇数の自然数）
100, 101, 102.01, 103.0301 …… 初項100、公比1.01の等比数列（元金100円に対する1％複利の利息）

　等差数列や等比数列にはルールがあるので、上のように数字を並べるのではなく、ルールを数式で表しておいた方が正確で、応用の幅も広がります。等差数列の第n項の値つまりa_nは、初項をa₁、公差をdとすると、以下のように表されます。

　念のため確認しておきましょう。以下、穴埋めになっているのでちょっと考えてみてください。オレンジ色の部分をクリックすると答えが表示されます。例えば、上の例で3番目の奇数にあたる値は、初項a₁＝1、公差d＝2、n＝3なので、

　　a₃ ＝　1　＋ (　3　－1)× 　2　＝ 5

となり、ちゃんと第3項の5が求められます。

　一方、等比数列の第n項は公比をrとすると、以下のようになります。

　こちらも、上の例で確認しておきましょう。3番目の項は、初項a₁＝100、公比r＝1.01、n＝3なので、

　　a₃ ＝　100　× 　1.01　^3　－1 ＝ 102.01

となり、ちゃんと第2項の102.01が求められました（初項に預金額を、公比に1＋利率を、nに利払い回数を指定すれば、定期預金の元利合計の計算ができますよ）。

　では、ルールが分かったところで、軽くウオーミングアップといきましょう。奇数の自然数を1から10個分、リストとして返す関数を作ってみてください。なお、これ以降のサンプルプログラムについては作成から実行までの流れを動画で見られるようにしてあります。以下の例を含めて、プログラミングに不慣れな方に向けての動画はこちらにまとめてあるのでぜひ参考にしてください。

動画1　初心者向けのプログラミングの基本（等差数列、総和、平均、個数の求め方）

　「初項をaという変数で表して1を代入する」「公差をdという変数で表して2を代入する」といった具合に、数式で表されている内容を一つ一つ丁寧にコードとして書いていきます。以下のリスト1では、できるだけ公式がそのままの形で表されるようにコードを書いてみました。空欄になっている箇所（_____の部分）を埋めてみてください。

　各行の末尾にある#以降の文章はコメントなので、プログラムの実行には影響を及ぼしません。実際にコードを入力するときには、コメントの入力は省略してもらっても構いません（が、必要に応じて書いておいた方が、コードが分かりやすくなります）。

def makeodds(m):
numbers = [] # 結果を入れるための空のリストを用意する
a = 1 # 初項
d = 2 # 公差
for n in range(1, m+1): # 1～第m項まで
numbers.append(_____) # 第n項を求めてリストに追加する
return numbers

makeodds(10)
# 出力例：[1, 3, 5, 7, 9, 11, 13, 15, 17, 19]

リスト1　等差数列のリスト（1以上の奇数の数列）を作る
makeodds関数の仮引数mには、作成する数列の項の個数を与える。最初にnumbersという名前の空のリストを作っておく。for文によってnの値を1からmまで変えながら、繰り返しを行う。繰り返しの中では、等差数列の第n項を求める公式を使って計算した値をnumbersというリストに追加する。appendはリストに値を追加するメソッド。

*　内容とは関係ありませんが、ここでは「リスト」という用語が混同しやすいので、気をつけてください。タイトルに記してある「リスト1」というのはプログラムのコードを意味します。一方「等差数列のリスト」や「空のリスト」の「リスト」は、Pythonで取り扱われるデータのことを表します。Pythonのリストとは上の出力例のように[]で囲まれた幾つかの値の並びのことです。

　では、答えを記しておきましょう。

（答え）　　a + (n-1)*d　

　公式のままですね。

【まとめ】for文による繰り返し、range関数による値の並びの作成

　for文の書き方やrange関数の使い方にはちょっと注意が必要です。既にfor文やrange関数になじんでいる方はこのまとめは読み飛ばしてもらっても構いません（でも、意外な発見があるかもしれないので、斜め読み程度には目を通しておいてもらうのがいいと思います）。

　例によって、ここでは必要最小限のまとめと、数学との関連、留意点のみに絞って説明しているので、Python文法のより詳細な内容については、『Python入門』をご参照ください。

　では、for文の書き方を簡単な例で図解してみましょう。

図1　for文の書き方
for文ではinの後に書かれた値の並びから値を1つ取り出して変数に代入し、文を実行する。続いてまた値の並びから値を1つ取り出して変数に代入し、文を実行する。このような処理を繰り返し、値の並びから取り出せるものがなくなったらfor文を抜け、次に進む。for文の1行目の最後の:を忘れないように。また、繰り返して実行したい文はインデント（字下げ）してfor文の範囲内にあることを表す。

　続いて、range関数の使い方です。range(a, b)とすると、a以上、b未満の整数を表す値の並び（rangeオブジェクト）が得られます*2。b以下ではなくb未満であることに注意してください。なお、aを省略してrange(b)とすると、0以上、b未満となります。

*2　rangeオブジェクトには実際の値の並びではなく、数列のルールだけが記憶されています。値が必要になれば、ルールを基に計算を行って値を返してくれます。一方、リストは実際の値を記憶しています。従って、ルールが決まっている数列であれば、リストにするよりもrangeオブジェクトにした方がメモリの消費が少なくて済みます。

　図1のコードを実行するとどういう結果になるでしょうか。range(0, 5)は、0以上5未満なので、iには順に0, 1, 2, 3, 4が代入され、それらの2乗つまり、0, 1, 4, 9, 16が順に出力されます。

　ところで、range関数は、range(a, b, c)のように指定すると「cずつ増やす」という意味になります。従って、リスト1の奇数のリストは以下の1行のコードで作成することもできます。

list(range(1, 20, 2))
# 出力例：[1, 3, 5, 7, 9, 11, 13, 15, 17, 19]

リスト2　等差数列のリスト（1以上の奇数の数列）を作る（違ったやり方）
range関数だけで等差数列が作れる。2番目の引数には、数列の要素数ではなく、最後の項の値＋1を指定していることに注意（19までを作成したいので、19＋1＝20未満とした）。全ての値を表示するにはlist関数を使ってリストに変換するのが手っ取り早い。

　実は、range関数はまさに等差数列を作るための関数なのです。ただし、作成できるのは整数の等差数列に限ります。従って、どの引数にも0.5などの小数部がある数値は指定できません。

　なお、for文の後にはrange関数だけでなく、リストや文字列などを指定することもできます。以下のコードを実行すれば、最初に見たAppleの株価が2021年8月5日～11日の間に結局どれだけ上下したかが計算できます。

diff = [0.11, -0.7, -0.05, -0.49, 0.26]
rise = 0
for x in diff: # リストの値を順に処理する
rise += x

rise
#出力例:-0.8699999999999999

リスト3　リストの値を総合計する
この例ではfor文のinの後にリストが指定されている。従って、最初はxに0.11が代入され、riseに足される。次はxに-0.7が代入され、riseに足される。以下同様に、リストの値がなくなるまで繰り返される。+=はそれまでの左辺の値に右辺の値を足した値を左辺の値とする演算子で、左辺の値に右辺の値を「足し込む」といった働きになる。この例では、5日間に株価が0.87ドル下がったことが分かる(-0.86999...となっているのは浮動小数点数の計算に伴う丸め誤差)。

　リスト3は合計を求めるためのお決まりのパターンです。最初に合計を求めるための変数に0を入れておき、繰り返し処理の中で値を次々と足していくというわけです。これについては、後でもう一度説明します。

コラム　数学の文字とプログラミングの変数名の微妙な違い

　数学では、何番目かの項を表すのにnという文字を使って「第n項」のように表すことがよくあります。また、a₁,a₂,...などの添え字は1から始まる数を使うことがよくあります。つまり1番目の値、2番目の値と考えるわけです。一方、プログラミングでは、nはデータの個数を表すのに使うことが多いようです。また、繰り返し処理の制御に使う変数やリストのインデックスを表す変数としてはi,j,kといった文字を使うことが一般的です。さらに、リストのインデックスは0番から始まります。従って、プログラミングに慣れた人にとって、リスト1は以下のように書いた方が自然に感じられるかもしれません。太字部分が変更箇所です。その場合、数式に書かれている文字とコード内の変数がそのまま対応しなくなることもあるので、注意が必要です。

def makeodds(n): # n個の奇数を作る関数
numbers = []
a = 1 # 初項
d = 2 # 公差
for i in range(0, n): # 0～第n-1項まで
numbers.append(a + i*d) # 公式を少し書き換えないといけない
return numbers

makeodds(10)
# 出力例：[1, 3, 5, 7, 9, 11, 13, 15, 17, 19]

リスト4　データの個数をnとし、初項が0番目の要素であることを意識したコード
前掲のリスト1と比較してみると、makeodds関数の仮引数がmからnに書き換わっていることが分かる（プログラミングでは「個数」を表すのに、nを使うことが多いので）。リストのインデックスを表すための変数にはiやjなどを使うことが多いので、for文の制御変数をnからiに変更してある。また、リストのインデックスは0から始まるので、range関数の第1引数を1ではなく0にした。そのため、第2引数がm+1からnに書き換わっている。append関数に指定していた等差数列の公式に沿ったa + (n-1)*dという記述は、上のように書き換えたことに伴い、a + i*dに変わっている。

　コードをどのように表現するかということに関しては、状況に依存するところもありますが、数式を素直に表現し、その意味がコードから分かるようにするならリスト1の書き方になりそうです。一方、問題やアルゴリズムがプログラミング言語を意識した表現で示されているなら、リスト4のような書き方の方が自然かと思われます。いずれにしても、変数名や変数の意味などは、明確にしておく必要があります。

2. 調和数列の和を求めるプログラム

動画2　調和数列とオイラーのγの近似値

　続いて、数列の和を求めるプログラムを見ていきます。と言っても、等差数列や等比数列の第n項までの和を求める場合は公式が使えるので、単純にそれらの公式を使って計算すれば済みます。ここではこれ以上詳しく取り扱いませんが、一応、公式だけ記しておきます。なお、数列の和のことを級数と呼ぶこともあります。

等差数列の第n項までの和S_n

　［A］では、1章で見たa_nを求める公式を代入しています。

等比数列の第n項までの和S_n

　さて、いよいよここからが本題です。大事なのでもう一回、言っておきましょう。ここからが本題です。目標のところに示した数式の一部は調和数列と呼ばれる数列の和になっています。

　調和数列とは、等差数列の逆数を並べた数列です。例えば、以下のような数列が調和数列です。

　この例は、1以上の奇数の数列（初項1、公差2の等差数列）の逆数の数列です*3。

*3　調和数列なんてどこに現れるんだと思われる向きもあるかもしれません。実は、音楽の世界に調和数列が登場します。ピタゴラス音律や純正律のド・ソ・ドという和音では、それらの間隔が調和数列になっています。この場合、下のドの弦の長さを1としたとき、ソは2/3、上のドは1/2となります。それらの逆数を取ると（周波数の比を求めると）、1, 3/2, 2、つまり2/2, 3/2, 4/2となります。これは初項1、公差1/2の等差数列になっています。次の動画はド・ソ・ドが調和数列になっているヴェルクマイスターIIIと呼ばれる音律でその和音を弾いてみた例です（が、著者には音感がないので平均律との違いが全く分かりません）。

動画3　調和数列の和音

　さて、（1）式の中に含まれている調和数列の和は、自然数（初項1、公差1の等差数列）の逆数の数列の和になっていることが分かりますね。その部分だけを取り出して、どのような計算をしているのかが分かるように足し算で表してみましょう。

　実は等差数列や等比数列と違って、調和数列には和を求める公式がありません。そのため、最初から各項を順に足していくしかありません。しかし、このことは、Σによる総和の計算をプログラミングするということにほかなりません（と、大げさに言いましたが、基本中の基本です）。ここでは、初項から第100項までを足してみましょう。もちろん、+演算子で100個の項をつなぐわけではなく、繰り返し処理を使います。

def sumharmonic(m):
sum = 0
for k in range(1, m+1):
sum += 1 / k
return sum

sumharmonic(100)
# 出力例：5.187377517639621

リスト5　調和数列の和を求める
総合計（累計）を求めるお決まりのパターン。最初に、合計を求めるための変数sumに0を代入しておき、繰り返し処理の中で1つずつ値をsumに足していく。繰り返しが終わった時点でsumの値が全ての合計になっている。

コラム　個数を求めるためのお決まりのパターン

　合計を求める方法は、個数を求める場合にも使えます。対象となるデータが見つかったときに、値そのものではなく1を足していけばいいのです。例えば、アンケート調査などで、0～4までの5つの選択肢のうち、どの番号を選んだ人が何人いるかを集計するものとします（つまり、頻度表を作ることになります）。以下のようなコードを書くといいでしょう。

data = [0, 3, 2, 1, 2, 4, 1, 3, 3, 2] # 回答のリスト
freq = [0]*5 # 頻度表。最初は全ての値を0とする
for sel in data: # 回答を順に処理
freq[sel] += 1 # 頻度表[回答]を増やす

freq
# 出力例：[1, 2, 3, 3, 1]

リスト6　頻度表を作成する
dataは収集された回答データ。最初の人は0を選び、次の人は3を選び……といった具合に並んでいる。freqは頻度表。freq[0]には0を選んだ人の人数を、freq[1]には1を選んだ人の人数を……という具合に記録するが、最初は全て0を入れておく。[0]*5は、0という値が5つ並んだリストという意味。

　for文の処理は、10人の人が順に番号のついたボールを5つの箱に入れるのを想像してもらえれば容易に理解できると思います。最初の人は⓪のボールを0番の箱に入れます。次の人は③のボールを3番の箱に入れます。同様に最後の人まで繰り返すというわけです。最後に、箱の中のボールの数を数えればどの番号を何人の人が選んだかという頻度表が作成できますね。

　この例では、選択肢の番号がリストのインデックスとうまく対応しているので、簡潔なコードが書けました。しかし、ルールがもう少し複雑な場合には、if文などを使って場合分けをしたり、計算によってインデックスを求めたりする必要があります。が、if文については回を改めて見ていこうと思います。計算によってインデックスを求める例は後で見ていきます。

3.　オイラーのγの近似値を求めるプログラム

　2章で見た調和数列の和は、mの値をどんどん大きくすると、無限大に発散することが分かっています。しかし、最初に見た以下の式は、オイラーのγ（＝0.5772156...）に収束します。

　前述したように、プログラムで無限に足し算を行うことはできないので、近似値を求めるためにmに好きな値を指定して計算できるようにしましょう。従って、limの部分（mの値を無限大に近づけていく）は気にしなくてもいいということになります。式は以下のようになりますね。

　なお、lnというのはlog_enのことです。eを底とする対数（自然対数）はよく使われるので、lnと略して書くことがよくあります（lnはlog naturalの略です）。また、（2）式の≃は「ニアリーイコール」と読み、ほぼ等しいということを表します。ともあれ、これならリスト5を少し書き換えるだけでできますね。

import math

def eulersgamma(m):
sum = 0
for k in range(1, m+1):
sum += 1 / k
return sum - math.log(m)

eulersgamma(10000)
# 出力例：0.5772656640681646