1のn乗根
正 $17$ 角形が定規とコンパスだけで作図可能であることを19歳のガウスが証明したのは有名な話です．また，正 $7$ 角形や正 $13$ 角形などは作図不可能であることが知られています．作図が可能とか可能でないなどの議論は体論やガロア理論によって解決されました．$2^n+1$ という形の素数をフェルマー素数というのですが，$p$ を素数とするとき，正 $p$ 角形が作図可能であることの必要十分条件は $p$ がフェルマー素数であることが知られています．現時点で知られているフェルマー素数は
$$
3, 5, 17, 257, 65537
$$
の5つだけです．ですから，頂点の数が素数個であるような作図可能正多角形はたった $5$ 種類しかないことになります．
数学的な細かい話を一切抜きにすると，一般に正 $n$ 角形が（目盛りのない）定規とコンパスで作図可能かどうかは，$1$ の $n$ 乗根（のすべて）が平方根だけを用いて表せるかどうか（ $\fallingdotseq$ 二次方程式を繰り返し解くという操作で得られるか）と同じ問題です．ここで，$1$ の $n$ 乗根とは $n$ 乗して $1$ になる数のことですから，方程式
$$x^n = 1 \Longleftrightarrow x^n - 1 = 0$$
の解です．例えば $n=3$ のときは
$$
x^3-1 = (x-1)(x^2+x+1)
$$
より
$$
x=1, \frac{-1\pm\sqrt{3}i}{2}
$$
が $1$ の$3$ 乗根です．$n=5$ とすると
$$
x^5 - 1 = (x-1)(x^4+x^3+x^2+x+1) = 0
$$
の解が $1$ の $5$ 乗根になります．$x^4+x^3+x^2+x+1=0$ を解くのは，高校生にとってはまあまあの問題です．ここでは $x^5-1=0$ をSymPyに解かせてみましょう．
from sympy import solve, symbols, init_printing
init_printing()
x = symbols('x')
solve(x**5-1)
$$\left [ 1, \quad - \frac{1}{4} + \frac{\sqrt{5}}{4} - i \sqrt{\frac{\sqrt{5}}{8} + \frac{5}{8}}, \quad - \frac{1}{4} + \frac{\sqrt{5}}{4} + i \sqrt{\frac{\sqrt{5}}{8} + \frac{5}{8}}, \quad - \frac{\sqrt{5}}{4} - \frac{1}{4} - i \sqrt{- \frac{\sqrt{5}}{8} + \frac{5}{8}}, \quad - \frac{\sqrt{5}}{4} - \frac{1}{4} + i \sqrt{- \frac{\sqrt{5}}{8} + \frac{5}{8}}\right ]$$
いずれも平方根を用いて表せている（$i=\sqrt{-1}$）ので，正 $5$ 角形は作図可能，という議論が成立します．
さて，高校数学の複素数平面の章で勉強することですが，一般に $1$ の $n$ 乗根は
$$
\cos \frac{2k\pi}{n} + i\sin \frac{2k\pi}{n}\quad(k=0, 1, 2, \ldots, n-1)
$$
で表せます．オイラーの公式
$$
e^{\theta i} = \cos\theta + i\sin\theta
$$
を用いれば
$$
e^{2k\pi i/n}\quad(k=0, 1, 2,\ldots, n-1)
$$
とも表せます．$e^{2\pi i} = 1$ であることと指数法則を考えれば，この数を $n$ 乗すれば $1$ になるのは明らかです．例として，$1$ の $13$ 乗根を計算してみましょう．
from sympy import exp, pi, I, N #I: 虚数単位 N: 数値計算
n = 13
for k in range(n):
rt = exp(2*k*pi*I/n)
print(N(rt))
1.00000000000000
0.88545602565321 + 0.464723172043769*I
0.568064746731156 + 0.822983865893656*I
0.120536680255323 + 0.992708874098054*I
-0.354604887042536 + 0.935016242685415*I
-0.748510748171101 + 0.663122658240795*I
-0.970941817426052 + 0.239315664287558*I
-0.970941817426052 - 0.239315664287558*I
-0.748510748171101 - 0.663122658240795*I
-0.354604887042536 - 0.935016242685415*I
0.120536680255323 - 0.992708874098054*I
0.568064746731156 - 0.822983865893656*I
0.88545602565321 - 0.464723172043769*I
これを $13$ 乗して $1$ になるか確認します．
from sympy import expand
for k in range(n):
rt = exp(2*k*pi*I/n)
print(expand(N(rt)**n))
1.00000000000000
0.999999999999999 + 1.77635683940025e-15*I
1.0 + 1.77635683940025e-15*I
1.0 + 7.45931094670027e-17*I
1.0 - 2.22044604925031e-15*I
1.0 - 3.5527136788005e-15*I
1.0 - 9.0205620750794e-16*I
1.0 + 9.0205620750794e-16*I
1.0 + 3.5527136788005e-15*I
1.0 + 2.22044604925031e-15*I
1.0 - 7.45931094670027e-17*I
1.0 - 1.77635683940025e-15*I
0.999999999999999 - 1.77635683940025e-15*I
いずれも虚部の部分は $10^{-15}$ 乗ほどで，ほぼ $0$ とみなせるので，全体はほぼ $1$ になっています．
実際問題，$\sqrt{2}$ などは循環しない小数，いわゆる無理数です．無理数を有限小数で近似することは厳密性を欠くことになりますので，数学的には，さっきの $1$ の $5$ 乗根のような，平方根などを用いた厳密な表記，すなわち真値を期待します．ガウスの偉業は，$1$ の $17$ 乗根を三角関数を用いずに平方根だけを用いて表すことができるというものです．$1$ の $17$ 乗根は $x^{17}-1 = 0$ の解ですから
$$(x-1)(x^{16}+x^{15}+\cdots + 1) = 0$$
の解を求めることになります．$x=1$ は明らかとして，$16$ 次方程式
$$x^{16} + x^{15} + \cdots + 1 = 0$$
の解法はそれなりに難しく，初等的な数論という位置付けにあるものの，体論やガロア理論を勉強すると理解が深まったりします．
実際に1の17乗根を計算してみよう
本当は細かく解説もしたいんですが，断念しました．コードだけを示します．
from sympy import symbols, expand, simplify, init_printing, N, sqrt, exp, pi, I, Integer, primitive_root
init_printing()
class LoopedTuple(tuple): # t=(1, 2, 3) だったら，t[3] で 1 が取り出せるタプル
def __getitem__(self, n):
return super(LoopedTuple, self).__getitem__(n % len(self))
def cluster_by_mod(ns, m): # 剰余類に分ける関数
return [tuple(ns[i::m]) for i in range(m)]
def divisor_sequence(n): # ガロア群の正規部分群の列の位数のリストを返す関数
from sympy import divisors
ds = divisors(n)
seq = [ds.pop(-1)]
for d in reversed(ds):
if seq[-1] % d == 0:
seq.append(d)
return seq
#\zetaの多項式の次数を下げる関数
def zreduce(expr, zeta, p, strict=True):
if expr.args == ():
return expr
else:
if expr.__class__.__name__ == 'Pow':
b, n = expr.args
if b == zeta and n.is_integer:
r = n % p
if r == 0:
return Integer(1)
elif strict and r == p-1:
return (sum([-zeta**i for i in range(p-1)]))
else:
return zeta**r
return expr.__class__(*[zreduce(arg, zeta, p, strict) for arg in expr.args])
# 部分体の元を部分体の基底で表す関数
def as_linear_combination(expr, zeta, p, to):
result = zreduce(expand(expr), zeta, p)
result2 = zreduce(expand(expr), zeta, p, strict=False)
if len(result2.args) <= len(result.args):
result = result2
for i in reversed(range(1, p)):
result = result.subs(zeta**i, to[zeta**i])
return result
# 二次方程式の解を部分体の基底で表す関数
def solve_quad_eq(a, b, zeta, p, to1, to2, square_expanding = False):
from sympy import re as real_part
from sympy import im as imaginary_part
Na, Nb = a, b
for i in reversed(range(1, p)):
Na = Na.subs(zeta**i, exp(2*pi*i*I/p))
Nb = Nb.subs(zeta**i, exp(2*pi*i*I/p))
alc = lambda expr, to: as_linear_combination(expr, zeta, p, to)
s = alc(a + b, to2)
prod = alc(a * b, to2)
if square_expanding == True:
s2 = alc(expand((a+b)**2), to2)
D = alc((s2-4*prod).subs(to1), to2)
else:
D = alc((s**2-4*prod).subs(to1), to2)
if real_part(Na) > real_part(Nb) or imaginary_part(Na) > imaginary_part(Nb):
return {alc(a, to1): (s+sqrt(D))/2, alc(b, to1): (s-sqrt(D))/2}
else:
return {alc(a, to1): (s-sqrt(D))/2, alc(b, to1): (s+sqrt(D))/2}
# 置換の繰り返しで真値を求める関数
def substitute(expr, zeta, p, to_indexed_basis, to_root):
result = expr
for i in reversed(range(1, p)):
result = result.subs(zeta**i, to_indexed_basis[zeta**i])
for to in to_root:
result = result.subs(to)
return result
def calc_pthroot_of_1(p, only_one = False): #p: フェルマー素数 理論上 3, 5, 17, 257, 65537 のみ入れられる．
pr = primitive_root(p) #原始根 ガロア群の生成元に関連
zeta = symbols(r'\zeta') # 1以外の1の17乗根の1つ（偏角が最小のもの）
# Q(\zeta)の基底
basis = LoopedTuple([zeta**(pr**i % p) for i in range(p-1)])
# Q(\zeta)の基底の添字表記
indexed_basis = LoopedTuple([symbols(r'\zeta_{{{i}}}'.format(i=i)) for i in range(p-1)])
# 添字の基底 -> 基底 の変換用
to_basis = dict(zip(indexed_basis, basis))
# 基底 -> 添字の基底 の変換用
to_indexed_basis = dict(zip(basis, indexed_basis))
# Q(\zeta)の部分体の基底
subbases = [LoopedTuple([sum(t) for t in cluster_by_mod(list(basis), d)]) for d in divisor_sequence(p-1)]
# 部分体の基底に用いる文字のリスト
letters = r'\zeta \alpha \beta \gamma \delta \varepsilon \eta \theta \iota \kappa \lambda \mu \nu \xi \omicron \pi \rho \sigma \tau \upsilon \phi \xi \psi \omega'.split()
# Q(\zeta)の部分体の添字基底
indexed_subbases =[LoopedTuple([symbols(letters[level]+r'_{'+str(i)+'}') for i in range(d)]) for level, d in enumerate(divisor_sequence(p-1))]
# 部分体の基底 <-> 部分体の添字基底 変換用
to_subbasis = {}
to_indexed_subbasis = []
for level, (subbasis, indexed_subbasis) in enumerate(zip(subbases, indexed_subbases)):
to_indexed_subbasis.append({})
for elem, indexed_elem in zip(subbasis, indexed_subbasis):
to_subbasis[indexed_elem] = elem
if level == 0:
to_indexed_subbasis[level][elem] = indexed_elem
else:
for arg in elem.args:
to_indexed_subbasis[level][arg] = indexed_elem - elem + arg
# 部分体の基底を実際に解く
to_root = []
for level, subbasis in enumerate(subbases[:-1]):
d = Integer(len(subbases[level]))/2
roots = {}
for i in range(d):
a = subbases[level][i]
b = subbases[level][d+i]
roots.update(solve_quad_eq(a, b, zeta, p, to_indexed_subbasis[level], to_indexed_subbasis[level+1], square_expanding=True))
to_root.append(roots)
if only_one == True: #\zetaのみreturn
return substitute(zeta, zeta, p, to_indexed_basis, to_root)
else: # 全ての真値
return [substitute(zeta**i, zeta, p, to_indexed_basis, to_root) for i in range(1, p)]
p = 3, 5 で試してみます．
p = 3
pthroots = calc_pthroot_of_1(p)
pthroots
$$\left [ - \frac{1}{2} + \frac{\sqrt{3} i}{2}, \quad - \frac{1}{2} - \frac{\sqrt{3} i}{2}\right ]$$
p = 5
pthroots = calc_pthroot_of_1(p)
pthroots
$$\left [ - \frac{1}{4} + \frac{\sqrt{5}}{4} + \frac{1}{2} \sqrt{- \frac{5}{2} - \frac{\sqrt{5}}{2}}, \quad - \frac{\sqrt{5}}{4} - \frac{1}{4} + \frac{1}{2} \sqrt{- \frac{5}{2} + \frac{\sqrt{5}}{2}}, \quad - \frac{\sqrt{5}}{4} - \frac{1}{4} - \frac{1}{2} \sqrt{- \frac{5}{2} + \frac{\sqrt{5}}{2}}, \quad - \frac{1}{4} + \frac{\sqrt{5}}{4} - \frac{1}{2} \sqrt{- \frac{5}{2} - \frac{\sqrt{5}}{2}}\right ]$$
よさそうです．次は p = 17 でいきます．
p = 17
pthroots = calc_pthroot_of_1(p)
pthroots[0] # pthroots全てを表示すると，とんでもなく長い数式が表示されるので注意
$$- \frac{1}{16} + \frac{\sqrt{17}}{16} + \frac{1}{8} \sqrt{- \frac{\sqrt{17}}{2} + \frac{17}{2}} + \frac{1}{4} \sqrt{- \sqrt{\frac{\sqrt{17}}{2} + \frac{17}{2}} - \frac{1}{2} \sqrt{- \frac{\sqrt{17}}{2} + \frac{17}{2}} + \frac{3 \sqrt{17}}{4} + \frac{17}{4}} + \frac{1}{2} \sqrt{- \frac{17}{8} - \frac{1}{4} \sqrt{- \frac{\sqrt{17}}{2} + \frac{17}{2}} + \frac{\sqrt{17}}{8} + \frac{1}{2} \sqrt{\frac{1}{2} \sqrt{- \frac{\sqrt{17}}{2} + \frac{17}{2}} + \frac{3 \sqrt{17}}{4} + \sqrt{\frac{\sqrt{17}}{2} + \frac{17}{2}} + \frac{17}{4}}}$$
pthrootsに $1$ 以外の $16$ 個の $1$ の $17$ 乗根が入っています．近似値を見てましょう．
for k in range(p-1):
print(N(pthroots[k]))
0.932472229404356 + 0.361241666187153*I
0.739008917220659 + 0.673695643646557*I
0.445738355776538 + 0.895163291355062*I
0.092268359463302 + 0.995734176295034*I
-0.273662990072083 + 0.961825643172819*I
-0.602634636379256 + 0.798017227280239*I
-0.850217135729614 + 0.526432162877356*I
-0.982973099683902 + 0.18374951781657*I
-0.982973099683902 - 0.18374951781657*I
-0.850217135729614 - 0.526432162877356*I
-0.602634636379256 - 0.798017227280239*I
-0.273662990072083 - 0.961825643172819*I
0.092268359463302 - 0.995734176295034*I
0.445738355776538 - 0.895163291355062*I
0.739008917220659 - 0.673695643646557*I
0.932472229404356 - 0.361241666187153*I
exp関数で求めたものと比較してみます．
for k in range(1, p):
print(N(exp(2*k*pi*I/p)))
0.932472229404356 + 0.361241666187153*I
0.739008917220659 + 0.673695643646557*I
0.445738355776538 + 0.895163291355062*I
0.092268359463302 + 0.995734176295034*I
-0.273662990072083 + 0.961825643172819*I
-0.602634636379256 + 0.798017227280239*I
-0.850217135729614 + 0.526432162877356*I
-0.982973099683902 + 0.18374951781657*I
-0.982973099683902 - 0.18374951781657*I
-0.850217135729614 - 0.526432162877356*I
-0.602634636379256 - 0.798017227280239*I
-0.273662990072083 - 0.961825643172819*I
0.092268359463302 - 0.995734176295034*I
0.445738355776538 - 0.895163291355062*I
0.739008917220659 - 0.673695643646557*I
0.932472229404356 - 0.361241666187153*I
あくまで近似値ですが，寸分も違わず一致しています．本当に $x^{16} + x^{15} + \cdots + 1 = 0$ の解かどうかは，最小多項式を見て確認するのも手です．
from sympy import minimal_polynomial
x = symbols('x')
minimal_polynomial(pthroots[0], x)
$$x^{16} + x^{15} + x^{14} + x^{13} + x^{12} + x^{11} + x^{10} + x^{9} + x^{8} + x^{7} + x^{6} + x^{5} + x^{4} + x^{3} + x^{2} + x + 1$$
間違いありませんでした．
1の257乗根を求めてみよう
第4のフェルマー素数 $257$ でも同様の議論は可能です．しかし，とんでもなく時間がかかるので，クラッシュしても責任は取れませんのご注意ください（下記はMac (Retina 5K, 27-inch, 2017; 4.2 GHz Intel Core i7)の結果）．
%time rt = calc_pthroot_of_1(257, only_one=True) #only_one=Trueオプションで1つだけ求める
CPU times: user 4min 11s, sys: 297 ms, total: 4min 12s
Wall time: 4min 12s
ちなみに手元のMacBookPro(Retina, 15-inch, Mid 2015; 2.2 GHz Intel Core i7)で%timeコマンドを付加して計算した結果です：
CPU times: user 7min 20s, sys: 379 ms, total: 7min 21s
Wall time: 7min 21s
このrtはinit_printingで表示させない方がいいと思います．
さて，近似値を計算して確認してみましょう．
N(rt)
$$0.999701157843094 + 0.0244457564247445 i$$
N(exp(2*pi*I/257))
$$0.999701157843094 + 0.0244457564247445 i$$
一致していました．感動です．
Mathematica でインポートするために，$\LaTeX$ 形式で保存したいのですが，これもまた凄まじいです．
from sympy import latex
src = latex(rt)
print(src)
※Qiitaに結果がプレビューされなかったので，スクリーンショットでお見せします（sqrtがスペルチェックに引っかかってます）．
現在のコンピュータでも計算が大変な対象を理論上計算できると頭脳だけで結論づけられるガウスはやはり偉大です．
機会があれば，各種計算の仕組みなども解説したいと考えていますが，数学的に内容が難しいので見送るかもしれません．