数学問題-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)