2017年早稲田数学をPythonで解く

## はじめに
この記事は、[Python Advent Calendar 2017](https://qiita.com/advent-calendar/2017/python_python) 12日目の記事です。

## ここでやること
PyCon JP 2017の会場でもとてもお世話になった早稲田大学の数学入試問題をPythonで扱います。

※ この記事で扱う問題は、[代々木ゼミナール（予備校） | 早稲田大学 基幹理工・創造理工・先進理工の入試問題と解答例（2016年解答速報）](http://sokuho.yozemi.ac.jp/sokuho/s_mondaitokaitou/1/mondai/mondai/1279671_4424.html)から引用しています。

過去の記事
* [2015年センター試験数学ⅡBをPythonで解く](https://qiita.com/massa142/items/b224a07c1b33e9a39cd8)
* [2016年東大数学をPythonで解く](https://qiita.com/massa142/items/3285a1b1ed12a6ab0299)

## 環境
+ Python 3.6.3

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

## 第５問
![waseda_0216_sugaku_mon5_01.jpg](https://qiita-image-store.s3.amazonaws.com/0/19111/6926716d-a18d-7c5b-cebc-94015e6a3a2c.jpeg "waseda_0216_sugaku_mon5_01.jpg")

### (1)
+ 変数p, q, αをsymbols()を用いて変数定義
+ 関数f(x)とg(x)を定義

```python
In [1]: import sympy as sy
In [2]: p, q, α = sy.symbols('p q α', real=True)
In [3]: def f(x):
            return x**3+x**2+p*x+q
In [4]: def g(x):
            return -1/(x+ 1)
```

> f(x)=0 の任意の解αに足して g(α）も f(x)=0 の解である

とあるので、

+ f(α) = 0
+ f(g(α)) = 0
+ f(g(g(α))) = 0

の3つの式が成り立つことがわかる。

```python
In [5]: c1 = f(α)
In [6]: c2 = f(g(α))
In [7]: c3 = f(g(g(α)))
```

solve()を用いて、この3つの式からp, q, αの値を求める。

```python
In [8]: sy.solve({c1, c2, c3}, {p, q, α})
Out[8]: [{p: -2,
		  q: -1,
		  α: 7**(2/3)*(-2*2**(1/3)/3 + (-4 + 2**(2/3)*(-1 + sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3))*(1 - sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3)/84)/((1 - sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))},
		 {p: -2,
		  q: -1,
		  α: -7**(2/3)*(2*2**(1/3)/3 + (1 + sqrt(3)*I)*(4 - 2**(2/3)*(-1 - sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3))*(7 + 21*sqrt(3)*I)**(1/3)/84)/((1 + sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))},
		 {p: -2,
		  q: -1,
		  α: -1/3 + 2**(1/3)*7**(2/3)/(3*(1 + 3*sqrt(3)*I)**(1/3)) + 28**(1/3)*(1 + 3*sqrt(3)*I)**(1/3)/6}]
```

となるので、(1)の答えは **p: -2, q: -1** となる。

これを元にf(x)を再定義しておく。

```python
In [9]: def f(x):
            return x**3+x**2-2*x-1
```
### (2)

(2)の証明はグラフによる図示で済ませればよいので、f(x)をmatplotlibを用いて描いてみる。

```python
In [10]: import matplotlib.pyplot as plt
In [11]: import numpy as np
In [12]: x = np.linspace(-2, 2, 100);
In [13]: plt.plot(x, f(x))
Out[13]: [<matplotlib.lines.Line2D at 0x10a735e80>]
```
![ダウンロード.png](https://qiita-image-store.s3.amazonaws.com/0/19111/5f0a642c-361e-333c-a4ee-3549c720dc40.png "ダウンロード.png")

-2<x<2の間にy=0となる点が3つあるので、

**f(x)=0は-2<x<2の範囲に3つの実数解を持つことがわかる**

### (3)(4)

とここまではすんなりいけたけど、(3)(4)は正直うまく解けなかった(´・ω・｀)

実際の試験だと

```math
2 \cos{\left (2 θ \right )} = 4 \cos^{2}{\left (θ \right )} - 2　
```
```math
2 \cos{\left (3 θ \right )} = 8 \cos^{3}{\left (θ \right )} - 6\cos{\left (θ \right )}
```

であることを頭に入れて、逆算から式をゴニョっていくのが速い気がする。

ここではせっかくPythonでやるんだから愚直なやり方で解こうと思って、f(x), f(g(x)), f(g(g(x)))に 2*sy.cos(θ) を突っ込んでみたら

```python
In [14]: θ = sy.symbols('θ', positive=True)
In [15]: sy.solve({f(2*sy.cos(θ)), f(g(2*sy.cos(θ))), f(g(g(2*sy.cos(θ))))})
Out[15]: [{θ: -acos(7**(2/3)*(-56*2**(1/3) + 7**(1/3)*(-4 - 2**(2/3)*7**(1/3)*(1 - sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))*(1 - sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))/(168*(1 - sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))) + 2*pi},
 {θ: -acos(7**(2/3)*(-56*2**(1/3) + 7**(1/3)*(-4 - 2**(2/3)*7**(1/3)*(1 + sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))*(1 + sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))/(168*(1 + sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3))) + 2*pi},
 {θ: -acos(-1/6 + 2**(1/3)*7**(2/3)/(6*(1 + 3*sqrt(3)*I)**(1/3)) + 28**(1/3)*(1 + 3*sqrt(3)*I)**(1/3)/12) + 2*pi},
 {θ: acos(7**(2/3)*(-2**(1/3)/3 + (-4 + 2**(2/3)*(-1 + sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3))*(1 - sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3)/168)/((1 - sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3)))},
 {θ: acos(-7**(2/3)*(2**(1/3)/3 + (1 + sqrt(3)*I)*(4 - 2**(2/3)*(-1 - sqrt(3)*I)*(7 + 21*sqrt(3)*I)**(1/3))*(7 + 21*sqrt(3)*I)**(1/3)/168)/((1 + sqrt(3)*I)*(1 + 3*sqrt(3)*I)**(1/3)))},
 {θ: acos(-1/6 + 2**(1/3)*7**(2/3)/(6*(1 + 3*sqrt(3)*I)**(1/3)) + 28**(1/3)*(1 + 3*sqrt(3)*I)**(1/3)/12)}]
```

となって、うまい感じの出力にはならなかった。

うーん困ったということで、[駿台の解答例](http://www2.sundai.ac.jp/yobi/sv/sundai/sokuhou/sokuhou_D/1337447054754.html)をカンニングしてみる。

すると、

```math
z = i \sin{\left (θ \right )} + \cos{\left (θ \right )}
```
```math
2 \cos{\left (θ \right )} = z + \frac{1}{z}
```

を利用すればよいことがわかる。

これを元に解いてみると、

```python
In [16]: z = sy.symbols('i z')
In [17]: sy.simplify(f(z+1/z))
Out[17]: (z**6 + z**5 + z**4 + z**3 + z**2 + z + 1)/z**3
In [18]: sy.solve(_)
Out[18]: [-cos(pi/7) - I*sin(pi/7),
 -cos(pi/7) + I*sin(pi/7),
 cos(2*pi/7) - I*sin(2*pi/7),
 cos(2*pi/7) + I*sin(2*pi/7),
 -cos(3*pi/7) - I*sin(3*pi/7),
 -cos(3*pi/7) + I*sin(3*pi/7)]
``` 

zに対して6つの解が得られる。

この中から

```math
z = i \sin{\left (θ \right )} + \cos{\left (θ \right )}
```

というzの形式に沿うものは、

* `cos(2*pi/7) + I*sin(2*pi/7)`
* `cos(4*pi/7) + I*sin(4*pi/7)` <= `-cos(3*pi/7) + I*sin(3*pi/7)`
* `cos(6*pi/7) + I*sin(6*pi/7)` <= `-cos(pi/7) + I*sin(pi/7)`

の3つに絞られる。

よって、(4)の答えは **θ=2/7π, 4/7π, 6/7π** である。

またこの回答から(3)の問の証明もなされる。

## おわりに

(3)のような問題は(4)を解くための誘導の意味合いが大きく、人間に最適化されているのでプログラミングで解くのは難しかった！