2016年東大数学をPythonで解く

## はじめに
この記事は、[Python Advent Calendar 2016](http://qiita.com/advent-calendar/2016/python) 11日目の記事です。

## ここでやること
[去年のPython Adevent Calendar](http://qiita.com/massa142/items/b224a07c1b33e9a39cd8)ではセンター試験の数学ⅡBを扱いましたが、今年は東大数学(文系)を扱っていきます。

※ この記事で扱う問題は、[河合塾（総合教育機関・予備校）/ 2016年度国公立大二次試験・私立大入試解答速報](http://kaisoku.kawai-juku.ac.jp/nyushi/honshi/16/t01.html)から引用しています。

## 環境
+ Python 3.5.0
+ Anaconda 3-2.4.0

## Sympyの基本的な使い方
### Sympyとは？
Pythonの記号計算ライブラリ
公式ドキュメント : http://www.sympy.org/en/index.html
日本語資料 : http://www.turbare.net/transl/scipy-lecture-notes/packages/sympy.html

### Symbols - 変数定義
```python
In [1]: from sympy import *
In [2]: x + 1
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-4cf92658b648> in <module>()
----> 1 x + 1

NameError: name 'x' is not defined

In [3]: x = symbols('x')
In [4]: x + 1
Out[4]: x + 1
```

### expand - 展開
```python
In [5]: expand((x + 1)**2)
Out[5]: x**2 + 2*x + 1
```

### factor - 因数分解
```python
In [6]: factor(x**4 - 3*x**2 + 1)
Out[6]: (1 + x - x**2)*(1 - x - x**2)
```

### simplify - 簡約
```python
In [7]: simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
Out[7]: x - 1
```

### limit - 極限
```python
In [8]: limit(x, x, oo)
Out[8]: oo
```

### diff - 微分
```python
In [9]: diff(cos(x), x)
Out[9]: -sin(x)
In [10]: diff(x**3 + x**2 - x - 1, x)
Out[10]: 3*x**2 + 2*x - 1
```

### integrate - 積分
```python
In [11]: integrate(cos(x), x)
Out[11]: sin(x)
In [12]: integrate(x**3 + x**2 - x - 1, x)
Out[12]: x**4/4 + x**3/3 - x**2/2 - x
```

### Matrix - 行列
```python
In [13]: Matrix([[1, 2, 3], [-2, 0, 4]])
Out[13]:
Matrix([
[ 1, 2, 3],
[-2, 0, 4]])
```

### solve - 式を解く
```python
In [14]: solve(x**2 - 1, x)
Out[14]: [-1, 1]
```

## 第1問
![1.gif](https://qiita-image-store.s3.amazonaws.com/0/19111/a216cbf5-1e49-c4e0-f80b-4713741cdaee.gif "1.gif")

+ 変数xyはSymbolsを用いて変数定義
+ 各座標PQRをMatrixを用いてベクトルに落とし込む

```python
In [1]: import sympy as sy
In [2]: x, y = sy.symbols('x y')
In [3]: P = sy.Matrix([x, y])
In [4]: Q = sy.Matrix([-x, -y])
In [5]: R = sy.Matrix([1, 0])
```

+ △PQRが鋭角三角形という条件をベクトルの内積を使って表現する
 	+ ∠QPRが90°以下の場合、$\vec{PQ}$・$\vec{PR}$ > 0
	+ ∠PQRが90°以下の場合、$\vec{QP}$・$\vec{QR}$ > 0 
	+ ∠PRQが90°以下の場合、$\vec{RP}$・$\vec{RQ}$ > 0
+ 内積と角度に関する詳しい内容を知りたい方には[「内積が見えると統計学も見える」第5回 プログラマのための数学勉強会 発表資料](http://www.slideshare.net/matsukenbook/5-55326268)がおすすめ
	+ $\vec{PQ}$・$\vec{PR}$ = PQ.norm() x PR.norm() x cos∠QPR
	+ `norm`はベクトルの長さ

```python
In [6]: sy.simplify((Q - P).dot(R - P) > 0)
Out[6]: 2*x*(x - 1) + 2*y**2 > 0
In [7]: sy.simplify((P - Q).dot(R - Q) > 0)
Out[7]: 2*x*(x + 1) + 2*y**2 > 0
In [8]: sy.simplify((P - R).dot(Q - R) > 0)
Out[8]: -x**2 - y**2 + 1 > 0
```

+ ここから得られた3つの式を、円の条件として見やすいようにリファクタリング

```math
2x(x-1)+2y^2 > 0 <=> (x-\frac{1}{2})^2+y^2 > \frac{1}{4}
```
```math
2x(x+1)+2y^2 > 0 <=> (x+\frac{1}{2})^2+y^2 > \frac{1}{4}
```
```math
-x^2 - y^2 + 1 > 0 <=> x^2+y^2 < 1
```

+ この条件を満たすグラフを描画

```python
In [9]: import matplotlib.pyplot as plt 
In [10]: fig = plt.figure()
In [11]: ax = plt.gca()
In [12]: ax.add_patch(plt.Circle((0,0),1,fc="#770000"))
Out[12]: <matplotlib.patches.Circle at 0x109689518>
In [13]: ax.add_patch(plt.Circle((0.5,0),0.5, fc="#FFFFFF"))
Out[13]: <matplotlib.patches.Circle at 0x109689f28>
In [14]: ax.add_patch(plt.Circle((-0.5,0),0.5, fc="#FFFFFF"))
Out[14]: <matplotlib.patches.Circle at 0x109696710>
```
+ X軸とY軸のアスペクト比を揃える
+ 軸の上限・下限を設定

```python
In [15]: ax.set_aspect('equal')
In [16]: plt.xlim([2, 2])
Out[16]: (-2, 2)
In [17]: plt.ylim([-2, 2])
Out[17]: (-2, 2)
```

+ 解答となるグラフを描画

```python
In [18]: plt.show()
```

![figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/19111/87afbc38-574d-bc4c-928a-328c418517ac.png "figure_1.png")

**解答: 点P(x, y)の範囲はグラフの赤色部分である**


## 参考書籍

[Pythonからはじめる数学入門](http://amzn.to/2gz5dyq)

今回紹介したSympyやmatplotlibについても紹介されてます。
数学プログラミングの取っ掛かりにおすすめです。

## さいごに

こんな感じでSympyを使うことで大学入試レベルだと簡単に解くことができます。
東大数学がこんな易しいレベルでいいのかという危惧もありますが...

ここでは第1問のみ扱いましたが、興味がある方はぜひ他の数学の問題もPythonで解いてみてください！