TIM Labs
とても簡単なリストの実験をしてみよう。
リストで繰り返しがあるとき、*で済ませると楽ができる。 とても簡単な例。
In [54]: A = [1,2,1,2,1,2]

In [55]: B = [1,2]*3

In [56]: print(A,B)
[1, 2, 1, 2, 1, 2] [1, 2, 1, 2, 1, 2]

In [57]: A==B
Out[57]: True

In [58]: A[3]+=5

In [59]: B[3]+=5

In [60]: print(A,B)
[1, 2, 1, 7, 1, 2] [1, 2, 1, 7, 1, 2]

In [61]: A==B
Out[61]: True
リストA,Bは当然同じ動きをする。 では、リストの要素を数字ではなく、リストに変えてみよう。 つまり、リストのリストの場合、どうなるだろうか?
In [80]: A = [[],[],[],[],[],[]]

In [81]: B = [[]]*6

In [82]: print(A); print(B); print(A==B)
[[], [], [], [], [], []]
[[], [], [], [], [], []]
True

In [83]: A[3].append(1)

In [84]: B[3].append(1)

In [85]: print(A); print(B); print(A==B)
[[], [], [], [1], [], []]
[[1], [1], [1], [1], [1], [1]]
False
リストAとBを2つの方法で作り、一致していることを確認した。 次に、A[3]とB[3]に1をappendした。 AとBは同じで、同じ操作をしたのだから、結果も同じになるはずなのだが、なっていない。
*6 により[]を6個にしても、実体は1つのままのようだ。

このでは、[]で実験したが、オブジェクトだったら全て同様の結果になるであろう。
オブジェクトを指している変数を、別の変数に代入しても同じオブジェクトを指すようになるだけというのが良く説明されているが、こういう場合にも同じ現象が当然起きる。

実験では、[]が6個だから、全部書いてもたいしたことはない。
もし、[]が100個だったら、n個だったらどうすればオブジェクトとして別々の[]を指定個数並べられるのだろうか?
from_zero_deep_learning (283x400).jpgだいぶ前に紹介した『ゼロからつくるDeep Learning』のコメントで、パラメータの勾配を求めるのに数値微分の方法が紹介され、ソースコードも載っている。

第4章ニューラルネットワークの学習
4.5 学習アルゴリズムの実装 (112ページ〜)

とても簡単な2層のニューラルネットワークを作って試そうという訳である。
学習に使うデータは、手書き数字のデータとして有名なMNISTデータセットを利用する。

以下が、今回試してみようという関数 train_neuralnet.py の一部である。

for i in range(iters_num):
    batch_mask = np.random.choice(train_size, batch_size)
    x_batch = x_train[batch_mask]
    t_batch = t_train[batch_mask]
    
    # 勾配の計算
    t0 = time.time()
    grad = network.numerical_gradient(x_batch, t_batch)    # 遅い
    #grad = network.gradient(x_batch, t_batch)             # 速い
    te = time.time()-t0
    
    # パラメータの更新
    for key in ('W1', 'b1', 'W2', 'b2'):
        network.params[key] -= learning_rate * grad[key]
    
    loss = network.loss(x_batch, t_batch)
    train_loss_list.append(loss)
    
    if i % iter_per_epoch == 0:
        train_acc = network.accuracy(x_train, t_train)
        test_acc = network.accuracy(x_test, t_test)
        train_acc_list.append(train_acc)
        test_acc_list.append(test_acc)
        print("{0:d}\t{1:f}\t{2:f}\t{3:f}".format(i,te,train_acc,test_acc))
この中の勾配の計算の部分の前後で時間を計測し、勾配の計算時間を求める。
結果の表示は、1エポック(データ全体を1回学習に使用すると1エポック)毎に表示している。
表示は、エポック番号、勾配計算時間、訓練精度、テスト精度になっている。

 117ページのヒント「鳥の絵」には、数値微分はとても遅いので、数値微分ではなくBackpropagation(誤差逆伝搬法)のプログラムも用意しているので、そちらを使って実験するようにと書かれていた。
コード中には、2つの勾配計算方法の呼び出しが併記されていた。

―――となると、つい、どのくらいの速度差が出るのか試してみたくなって、やってみた。

grad = network.gradient(x_batch, t_batch)

In [13]: run train_nueralnet.py
0	0.001145	0.104417	0.102800
1	0.001112	0.099150	0.100900
2	0.001019	0.112367	0.113500
3	0.001067	0.112367	0.113500
4	0.001076	0.112367	0.113500
5	0.001022	0.112367	0.113500
6	0.001126	0.104417	0.102800
7	0.001148	0.104417	0.102800
8	0.001018	0.089583	0.093700
9	0.001133	0.112367	0.113500


grad = network.numerical_gradient(x_batch, t_batch)

In [18]: run train_nueralnet.py
0	39.711862	0.098717	0.098000
1	40.883988	0.098717	0.098000
2	42.274781	0.098717	0.098000
3	40.617868	0.098717	0.098000
4	40.499471	0.098717	0.098000
5	40.557489	0.098717	0.098000
6	41.622460	0.112367	0.113500
7	40.826519	0.116450	0.118000
8	40.555552	0.112367	0.113500
9	40.555557	0.112367	0.113500
これを見ると、前半の誤差逆伝搬(高速)版では、1エポックが約1ミリ秒に対し、後半の数値微分(低速)版では40秒程度になっている。 つまり、約40,000倍速度向上が実現できている。 確かに、これは画期的である。
誤差逆伝搬法の説明は長くなるので、ここでは行わない。この本では直感的でわかりやすい計算グラフを使って説明している。
 この方法で計算を行えるようになり Deep Learning が使い物になった。
C言語では、
num = 12;
str = num > 10 ? "Big" : "Small";
という風に三項演算子が使えて便利だった。

Pythonで、これを
>>> num = 12
>>> if num > 10:
...     str = 'Big'
... else:
...     str = 'Small'
... 
>>> str
'Big'
と書くのは面倒だ。この4行のif文を1行にできないものだろうか?
Pythonの本を見ても、if の説明はこの4行タイプばかりである。

実は、これを1行で書くことが可能なのだ。
>>> num = 12
>>> str = 'Big' if num > 10 else 'Small'
>>> str
'Big'
>>> num = 8
>>> str = 'Big' if num > 10 else 'Small'
>>> str
'Small'
もちろん、代入文にせず、式だけでも書ける。
>>> num = 8
>>> 'Big' if num > 10 else 'Small'
'Small'
条件式の書式を書くと、次のようになる。
〈式t〉 if   〈条件式〉    else    〈式f〉
この書き方は、使いこなせれば便利なことは多い。それは、C言語での3項演算子(? :)と同じである。

以下は、良い子は読まない方が良いと思うが書いておく。

実は、これ以外にも1行で書く方法がある。といっても、かなりアレな方法である。
>>> num = 12
>>> ['Small','Big'][num>10]
'Big'
>>> num = 8
>>> ['Small','Big'][num>10]
'Small'
この書式は次のようになっている。
[〈式f〉, 〈式t〉] [〈条件式〉]   
条件式がTrueのときに1,Falseのときに0として扱われるので、この書き方が成立する。

この書き方でコンパクトになるが、非推奨である。
Pythonで円周率100万桁の計算をしてきたが、ちょっとだけ誤差がでてしまった。計算の桁数を十分に余裕を持って100万100桁まで計算して100万桁だけ採用すれば、きっと正しい値になるだろう。
でも、これでは面白くない。
あまり余分な桁数まで精度をあげることなく、つまり無駄なく計算をしたいものである。

さて、では、どういう風にやるのが良いだろうか。

各項を項番号nとしたとき、nの式で表した。
そして、n=0,1,2,3,....とやって、項が想定していた誤差epsより小さくなるまで足してきた。

しかし、このやり方は、とても「マズイ」。超高精度計算ではやってはいけない方法だ。
これをちょっとだけ説明しておこう。

ちょっとこんな実験をやってみた。
lamda関数をつい使ってみたが、lambda関数の説明は今回は省略する。とても便利なものだということだけ覚えておいてほしい。
In [30]: term = lambda n:(-1)**n/2**(10*n+6) * ( 
                - 2**5/(4*n+1) - 1/(4*n+3) + 2**8/(10*n+1) 
                - 2**6/(10*n+3) - 2**2/(10*n+5) - 2**2/(10*n+7) + 1/(10*n+9)  )

In [31]: term(0)
Out[31]: 3.1417658730158724
term(n)という関数を作って、初項を計算してみたが、すでに3.14...とかなり円周率に近い値になっている。 足し合わせた合計は、ますます円周率に近づくはずである。そして、各項は、どんどん小さくなっていく。
今のプログラムは横着をして、どの項数も同じ精度、有効数字で計算しているのだが、総和は3.14....で、足す数は非常に小さい数である。
足す数の有効数字が100万桁あっても、足された結果は3.14...としたときの有効桁数100万桁になるから、足す数は有効桁数が100万桁でも、桁合わせして足されるため、実際に足した後に反映される桁は一部である。
nが大きくなるにしたがって、実際の有効桁数はどんどん減ってしまう。
これが積もっていくと、誤差が何桁にもなってしまう。

ということで、足し方を変えよう。 前から足してダメなら、後ろから順番に足してみよう。
以下がそのプログラムである。
def FabriceBellard(len):
    var('a:z')
    term = (-1)**n/2**(10*n+6) * ( 
                - 2**5/(4*n+1) - 1/(4*n+3) + 2**8/(10*n+1) 
                - 2**6/(10*n+3) - 2**2/(10*n+5) - 2**2/(10*n+7) + 1/(10*n+9)  )
    eps  = 10**(-n).evalf(2,subs={n:len+2})

    nmax = 0;
    while True:
        t = term.evalf(10,subs={n:nmax})
        if abs(t)
最初に、何項になったらepsより小さくなるかを求めて、そこから初項まで逆順に加えている。
これで実行した結果を示す。
In [23]: len = 1000001

In [24]: st = time.time()

In [25]: p = FabriceBellard(len)

In [26]: time.time()-st
Out[26]: 41192.76380801201

In [27]: str(p)[-50:]
Out[27]: '56787961303311646283996346460422090106105779458151'

In [28]: str(sympy.pi.evalf(len))[-50:]
Out[28]: '56787961303311646283996346460422090106105779458151'
今回は、項の計算を2回行っているのだが、全体の計算時間は1回しかやらなかった前回よりも短くなっている。こういうところにも、後ろから足した効果が出ている。
計算結果は、sytmpy.pi の100万桁の値とちゃんと一致している。

この計算、もっと桁数をサボりながら100万桁計算することができるが、それくらいは各自でやってみよう。
桁数をどんどん増やしたいなら、はやりPythonでは限界がある。そういう場合は、C言語とかで頑張って欲しい。その場合でも、今回使った円周率の次の公式は役立つはずだ。 fb_pi.png 今から40年近く前、8ビットパソコン(マイコン)で1万桁を求めるのが流行ったりしたが、コンピュータの性能はどんどん上昇し、いまや個人でも1兆桁に挑戦できるくらいコンピュータの性能が向上している。

ということで、円周率の話はここで終わりとする。
CodingTheMatrixJP.gif行列プログラマー
――Pythonプログラムで学ぶ線形代数

Philip N. Klein 著
松田 晃一、弓林 司、脇本 佑紀、中田 洋、齋藤 大吾 訳
2016年10月 発行
688ページ、5400円(本体)
オライリー・ジャパン
ISBN978-4-87311-777-5




副題にPythonとあるので、Pythonで行列の計算処理をしっかりやるための本と受け取るのが普通と思うが、本書はそういう本ではない。
Pythonプログラミングの本というより、数学書に相当近い。

行列といえば、当然NumPyを使い倒すPythonプログラマ養成本と考えてしまうのだが、本書ではNumPyについては書かれていない。

数学書といったように、コードは含まれているのだが、それよりも数学、行列をしっかり教える本で、通常の行列計算以外にも、行列はいろりろ使い道があるなどの話題が多い。

逆に、普通のPythonの本ではあまり説明の無い内包についてはしつこく出てくる。
普通だったら、for文を使って、その後インデントしてという感じになるようなことを、内包を使って、ちょろっとコーディングしてしまおうという方針のようだ。そして、内包というのは、数学の表記と同じなので、とてもなじみやすいと書いてある。
といっても、それは数学者にとって、日常的な書き方がそのままPythonにも使えて便利ということで、普通の数学が苦手なプログラマには逆効果かも知れない。

本書の課題では、内包を使ってプログラムせよというのが多い。これをやっていると、次第に、def、:(コロン)、字下げなどが無いプログラムを書くようになり、Pythonらしくなくなる。

lambda関数も詳しく説明しているかと思ったら、こちらは入っていなかった。残念である。
内包とlambda関数を使い倒すと面白いと思うのだが。

ということで、本書は、かなり毛色の変わったPythonの本で、線形代数学にプログラミングが付加された本である。

プログラミングの普通でない側面を勉強したい場合には役に立つが、普通にPythonのプログラムを勉強しよう、Deep Learning の準備としてPythonの本を読もうと言う場合は、絶対に他書を読むべし。

―――と書いたが、まだ4割くらいしか読んでいない。

せっかく早く収束する円周率の公式を使って計算をはじめたのだが、永久ループに入ってしまったようだ。
永久ループから抜ける判断をしているのは、次の個所だけである。
        if abs(t)<10**(-len):
            break

これから考えられるのは、この条件式がいつまで経っても成立しないことだ。 tは高精度な計算をしているので、たぶん問題ないだろう。 限界を決めるのに、10のべき乗を利用しているのだが、もしかしてこれが問題だろうか。 それで、こんな実験をしてみた。
In [26]: 10**(-400)==0.0
Out[26]: True

In [27]: 10**(-300)==0.0
Out[27]: False
10**(-400)は、0.0と等しくなってしまったようだ。 ということは、
        if abs(t)<0.0:
という絶対に成立しないことをしていたようだ。 高精度の計算を行うには、もっと微小な値をちゃんと計算し、比較しないといけない。 どうすればよいだろうか?

ということで、微小な下限値を、tと同様に、まず式で求め、evalfで桁数を指定して評価するように変更してみた。
def FabriceBellard(len):
    var('a:z')
    term = (-1)**n/2**(10*n+6) * (-2**5/(4*n+1) - 1/(4*n+3) + 2**8/(10*n+1) 
        - 2**6/(10*n+3) - 2**2/(10*n+5) - 2**2/(10*n+7) + 1/(10*n+9))
    eps  = 10**(-n).evalf(2,subs={n:len+2})
    sigma = 0
    m = 0
    while True:
        t = term.evalf(len+2,subs={n:m})
        if abs(t)<eps:
            break
        sigma += t
        m += 1
    return sigma.evalf(len)
epsの値だが、余裕をもって2桁小さくした。
これで、1001桁で実行してみた。
In [51]: len = 1001

In [51]: str(FabriceBellard(len))[-50:]
Out[51]: '18577805321712268066130019278766111959092164201989'

In [52]: str(sympy.pi.evalf(len))[-50:]
Out[52]: '18577805321712268066130019278766111959092164201989'
大丈夫なようだ。
時間計測を含めて、ちょっと100万桁求めてみた。

In [68]: len = 1000001

In [69]: st = time.time()

In [70]: p = FabriceBellard(len)

In [71]: time.time()-st
Out[71]: 43364.538419008255

In [72]: str(p)[-50:]
Out[72]: '56787961303311646283996346460422090106105779458148'

In [73]: str(sympy.pi.evalf(len))[-50:]
Out[73]: '56787961303311646283996346460422090106105779458151'

43364秒かかった。約12時間。 一番最後の桁が3だけ違った。
まあいいか。
でもないかな。
1違うのは仕方がないにしても、3違うのは許容範囲を超える原因がありそうだ。

考えてみよう。
円周率を延々と求めることが延々と行われてきた。そして、効率的に計算するための公式が次々と考案されてきた。桁数を上げようとすると、急激に計算量が増えてしまうのだった。
そのあたりを説明するときりがないので、詳しくは円周率(Wikipedia)を見よ。

しかし、どうやら最近は非常に効率のよい公式が見つかっているようだ。たとえば、これだ。

fb_pi.pngこの収束は、2の10n乗分の1のペースのようだ。つまり、桁数nが増えても桁数に比例する程度の時間でしか計算量が増えないように思われる、なかなか画期的な公式だ。
この公式の導き方を知りたい場合は、次を読もう。
A new formula to compute the n'th binary digit of pi, Fabrice Bellard, January 20, 1997

この公式を利用して、円周率を求めるPythonの関数を作ってみた。

def FabriceBellard(len):
    var('a:z')
    term = (-1)**n/2**(10*n+6) * (-2**5/(4*n+1) - 1/(4*n+3) + 2**8/(10*n+1) 
               - 2**6/(10*n+3) - 2**2/(10*n+5) - 2**2/(10*n+7) + 1/(10*n+9))
    sigma = 0
    m = 0
    while True:
        t = term.evalf(len,subs={n:m})
        if abs(t)<10**(-len):
            break
        sigma += t
        m += 1
    return sigma

桁数をパラメータlen で指定する。
項の加算を打ち切る限界は、桁数から決めた。

まず、とりあえず小数点以下50桁(全体で51桁)を求めてみた。
In [11]: FabriceBellard(51)
Out[11]: 3.14159265358979323846264338327950288419716939937511

In [12]: sympy.pi.evalf(51)
Out[12]: 3.14159265358979323846264338327950288419716939937511

大丈夫なようだ。もっと桁数を増やしてみよう。
In [19]: str(FabriceBellard(301))[-50:]
Out[19]: '45648566923460348610454326648213393607260249141274'

In [20]: str(sympy.pi.evalf(301))[-50:]
Out[20]: '45648566923460348610454326648213393607260249141274'
300桁も問題なく動いているようだ。 さらに続けて、400桁も確かめてみよう。
In [22]: str(FabriceBellard(401))[-50:]
^C---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
しかし、400桁まで計算してみようとしたら、いつまで経っても、うんともすんとも言わないのでCtrl-Cで強制終了した。
何が悪いんだろう。
前回は、円周率を100万桁、とても効率的というか、インチキというか、そういう方法で求めた。

今回は、真面目に円周率を計算する方法を考えてみよう。

円周率は、随分前から計算するための公式が分かっており、wikipediaにたくさん載っている。

昔から有名なのは、次のマチンの公式であろう。

pi_matin.png
これに基づいて計算すると、この程度の精度は出る。
In [397]: 4 * (4*atan(1/5) - atan(1/239))
Out[397]: 3.1415926535897936

しかし、これでは単に計算しただけに過ぎない。もっと精度を上げる工夫をしよう。
In [417]: expr = 4*(4 * atan(1/5) - atan(1/239))

In [418]: expr.evalf(50)
Out[418]: 3.1415926535897935600871733186068013310432434082031
これで50桁表示できたようだ。前回のpiと比較してみよう。
In [419]: sympy.pi.evalf(50)
Out[419]: 3.1415926535897932384626433832795028841971693993751
ということで、全然一致していない。 どうやら、数式を計算してから、evalfで高精度にしても無駄のようだ。 なので、数式のままにしておいて、evalfの中で値を指定してみよう。
In [420]: expr = 4 * (4*atan(1/x) - atan(1/y))

In [421]: expr.evalf(50,subs={x:5,y:239})
Out[421]: 3.1415926535897932384626433832795028841971693993751
これで一応50桁は一致したようだ。 もうちょっと先の桁まで比較してみよう。 小数点以下100万桁に挑戦。
In [430]: str(expr.evalf(1000001,subs={x:5,y:239}))[-50:]
Out[430]: '56787961303311646283996346460422090106105779458151'

In [431]: str(sympy.pi.evalf(1000001))[-50:]
Out[431]: '56787961303311646283996346460422090106105779458151'
これで100万桁、ちゃんと計算できたようだ。めでたし、めでたし。
と思いたいところだが、単にatanを使って、桁数指定して計算しただけなので、やっぱりインチキ感は否めない。どうしよう。
前回、SymPyでの連分数を作る関数list_to_frac(A,B)を定義した。
これを利用して計算してみよう。
連分数は、有理数だけではなく、無理数、超越数なども表せる。

まず練習として、2の平方根を求めてみよう。

cf_root_2a.pngこれから、2つのリストを[1,2,...]と[1,1,...]を作ることで2の平方根が求まる。
In [91]: pa = [1]+[2]*20

In [92]: pa
Out[92]: [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

In [93]: pb = [1]*20

In [94]: pb
Out[94]: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

In [95]: list_to_frac(pa,pb)
Out[95]: 1.414213562373095

では、円周率を求めてみよう。Wikipediaに載っていた式をπについて整理すると次式を得る。

cf_pi.pngこれを元に円周率は簡単に計算できる。
n = 10

In [97]: pa = [0]+list(range(1,2*n,2))

In [98]: pa
Out[98]: [0, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19]

In [99]: pb = [4]+[x**2 for x in range(1,n+1)]

In [100]: pb
Out[100]: [4, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

In [101]: list_to_frac(pa,pb)
Out[101]: 3.14159254044654
円周率は、
In [103]: math.pi
Out[103]: 3.141592653589793
なので、精度がちょっと足りないので、25段にしてみよう。
In [104]: n = 25

In [105]: pa = [0]+list(range(1,2*n,2))

In [106]: pb = [4]+[x**2 for x in range(1,n+1)]

In [107]: list_to_frac(pa,pb)
Out[107]: 3.141592653589793
ちゃんと計算できているようだ。 円周率というと、延々と1万桁、100万桁、、、、と求められているが、どうやったら求まるのだろうか? 効率の良い計算方法がいろいろ発表されているが、連分数では効率が良くないので、これ以上深入りしない。

でも、それでは不満だろうから、ちょっと円周率を100万桁求められる方法を書いて終わりにする。
In [108]: import sympy

In [109]: sympy.pi.evalf(20)
Out[109]: 3.1415926535897932385
これで、evalf()で桁数を指定すると、その桁まで表示される。 なお、整数部(3.)の部分も含まれるので、小数点以下何桁というのと照合するときは桁数を1加える必要がある。

これを利用すると、1万桁でも、100万桁でも簡単に求まる。
以下では、最後の50桁だけ示した。
In [115]: str(sympy.pi.evalf(10001))[-50:]
Out[115]: '46101264836999892256959688159205600101655256375678'

In [116]: str(sympy.pi.evalf(1000001))[-50:]
Out[116]: '56787961303311646283996346460422090106105779458151'

円周率100万桁までの本としては、以下の本が一部では非常に有名です。


円周ankoku_pi.jpg率1,000,000桁表
牧野貴樹 著
暗黒通信団 (2016/03) 発行
314円(本体)
 ISBN-13: 978-4873100371
連分数は数学の中でもとりわけ美しい魅惑的な分野だ。
しかし、分数の計算は面倒というのが小学校以来ある。
それを、SymPyで解決してしまおう。ContinuedFraction4.png

まず、連分数を最初に示す。
こんな数式を扱いたいのである。

右のpprint(f)のfは連分数の式で、それをpprint(f)で美しくprintした。

この形の再帰形式になっている分数を連分数という。
b1,b2,.... が全て1の場合を正則というそうだ。

aとbの数列を与えて連分数を生成する関数は以下のようになる。
分母のAの方が、分子のBより1つだけ長い。

import math
from sympy import *

def list_to_frac(A,B):
    expr = 0
    zip(A[1:],B)
    for (a,b) in reversed(list(zip(A[1:],B))):
        expr = b/(a+expr)
    return A[0]+expr
In [22]: symsa = symbols('a0:5')

In [23]: symsa
Out[23]: (a0, a1, a2, a3, a4)

In [24]: symsb = symbols('b1:5')

In [25]: symsb
Out[25]: (b1, b2, b3, b4)

In [26]: f = list_to_frac(symsa,symsb)

In [27]: f
Out[27]: a0 + b1/(a1 + b2/(a2 + b3/(a3 + b4/a4)))
symbols()を使うことで、まとめて数列の添字付変数を作っている。

list_to_frac()に与える分母、分子の数列であるが、これはリストでもタプルでも構わない。pprint()を用いると美しい分数形式で表示されるのだが、単にprint()とかの場合には、カッコだらけのベタッとした感じに表示される。

連分数についてはネット上にもたくさん情報はあるが、BLUE BACKSに以下の入門書がある。

BLUEBACKS_contfraction.jpgBLUE BACKS - B1770
連分数の不思議
    無理数の発見から超越数まで

木村 俊一 著

新書版、328ページ
2012/05/21 発売
1,100円(税別)
講談社
ISBN-978-4-06-257770-0

最近のコメント