数学問題-二重根号の簡約（Python 版）

数学問題-3 で，「二重根号の簡約は，直接的に解を求める事が出来ない（simplify では解が得られない）」と書いたが，解を求めるプログラムを Python で書いてみた。エラーチェックや，数式を引数にするという書式にしたので，ちょっと長くなった。

使用例

>>> double_root("sqrt(6-sqrt(32))")
('2-sqrt(2)', 0.5857864376269049)
>>> double_root("sqrt(6-2*sqrt(8))")
('2-sqrt(2)', 0.5857864376269049)
>>> double_root("sqrt(106+2*sqrt(1288))")
('sqrt(92)+sqrt(14)', 13.33332043339938)
>>> double_root("sqrt(106+sqrt(5152))")
('sqrt(92)+sqrt(14)', 13.33332043339938)
>>> double_root("sqrt(106+4*sqrt(322))")
('sqrt(92)+sqrt(14)', 13.33332043339938)
>>> double_root("sqrt(106+sqrt(368))")
"can't simplify"
>>> double_root("sqrt(-16+sqrt(368))")
"can't simplify"
>>> double_root("sqrt(16-sqrt(-368))")
'square root of negative value is not allowed (inner sqrt)'
>>> double_root("sqrt(16+sqrt(-368))")
'square root of negative value is not allowed (inner sqrt)'
>>> double_root("sqrt(16-sqrt(368))")
'square root of negative value is not allowed (outer sqrt)'
>>> double_root("sqrt(18+2*sqrt(81))")
('6', 6)
>>> double_root("sqrt(18-2*sqrt(81))")
('0', 0)
>>> double_root("sqrt(13+sqrt(144))")
('5', 5)
>>> double_root("sqrt(13-sqrt(144))")
('1', 1)

プログラム

>>> import re
>>> import numpy as np
>>> from math import sqrt
>>> 
>>> def double_root(equation):
...     def solve(a_plus_b, a_mult_b):
...         for b in range(1, int(a_mult_b**0.5 + 1)):
...             if a_mult_b % b == 0 and a_plus_b == b + a_mult_b // b:
...                 return (a_mult_b // b, b)
...         return (np.nan, np.nan)
...     def simplify(a):
...         int_root = int(a**0.5)
...         if int_root**2 == a:
...             return str(int_root)
...         else:
...             return "sqrt(%d)" % a
...     inner = re.sub("sqrt\(", "", equation, 1)
...     inner = re.sub("\)", "", inner, 1)
...     match = re.search("sqrt\(-[0-9]*\)", equation)
...     if match is not None:
...         return "square root of negative value is not allowed (inner sqrt)"
...     elif eval(inner) < 0:
...         return "square root of negative value is not allowed (outer sqrt)"
...     result1 = eval(equation)
...     sign = '+' if re.search('-', equation) is None else '-'
...     equation = re.sub("\)", "", equation)
...     equation = re.sub("sqrt\(", "", equation)
...     equation = re.sub("[-+*]", ",", equation)
...     # print(equation)
...     if equation[0] == ',':
...         return "can't simplify"
...     parsed = list(map(int, equation.split(",")))
...     # print(parsed)
...     if len(parsed) == 3:
...         parsed[1] = parsed[1]**2 * parsed[2]
...     parsed[1] = parsed[1] / 4
...     if parsed[1] != int(parsed[1]):
...         return "can't simplify"
...     # print(parsed[1], parsed[2], "\n")
...     res = solve(parsed[0], int(parsed[1]))
...     if np.isnan(res[0]):
...         return 'can\'t simplify'
...     res_str = "%s%s%s" % (simplify(res[0]), sign, simplify(res[1]))
...     result2 = eval(res_str)
...     if np.allclose(result1, result2):
...         if int(result2) == result2:
...             res0, res1 = sqrt(res[0]), sqrt(res[1])
...             res_str = str(int(res0 + res1)) if sign == '+' else str(int(res0 - res1))
...         return res_str, eval(res_str)
...     return "error: original = %g, simplified = %g" % (result1, result2)