一番重大なエラーはどれか(Python)

とある質問サイトに，以下の投稿があった。そのサイトでは，適切に投稿しないとインデントが反映されないので，インデントは直しておいたが，もしかしたら，質問者のインデントも不適切な部分があったのかもしれない。

import math

number = input("正の整数を入力してください: ")

def trial_division(target):
 dest = int(math.sqrt(target))
 
 for i in range(2,dest):
  if target % 1 == 0:
   print(str(target) + "は素数ではありません")
   return
 print(str(target) + "は素数です")

trial_division(number)

というプログラミングを作成したのですが、エラーが出ました。
6行目と最後の行が間違っているようですが、なにが間違っているのでしょうか。 

唯一の回答者は，
　targetという変数が意味不明ですな
　numberじゃね？ 
　ってより、numberをtergetにすりゃいいってはなしかなー 
などと与太を飛ばしております。

あなたならわかりますよね。

#############################

論理エラーでしょうね。range(2,dest) は range(2, dest + 1) としないとだめ（Python 特有の罠）
target % 1 の分母は 1 ではなく i だというのは，本か何かのプログラムの写し間違い？（変数名がおかしいしなあ）
6行目 のエラーは math.sqrt(target) で起きているんだろうね。入力した number は，文字列だから，数値演算できない。

最後の行 の間違いというのは，input と def の位置が悪いと言うこと。

全部ダメだ

import math

def trial_division(number):
    dest = int(math.sqrt(number))
    for i in range(2, dest + 1):
        if number % i == 0:
            print(str(number) + "は素数ではありません")
            return
    print(str(number) + "は素数です")

number = int(input("正の整数を入力してください: "))
trial_division(number)
trial_division(23)
trial_division(24)
trial_division(25)
trial_division(1000000007)

モンテカルロ法により円周率を求める

Qiita では，これも定期的に出てくるやつ。

長々しいプログラムが散見されるが，以下のように簡潔に書く。

import numpy as np
def sim(n):
    xy = np.random.random((n, 2))**2
    print(np.mean(xy.sum(axis=1) < 1)*4)

< 1 のところを <= 1 にしないといけないのではないかとか，np.random.random は 0 以上，1 未満 の一様乱数なので，0 以上，1 以下の一様乱数を使わないといけないのではないかとか，おかしなことをいう人もいるが，そのようなおかしな人を相手にする必要はない。

三次元空間における2点を通る直線

二次元における話を三次元に拡張し，sympy で解いてみた。

x-y-z の三次元空間で，x-y 平面（床面としよう）が鏡になっている。z 軸が正の方向（x-y 平面の上方向）からこの鏡にレーザー光を当てると反射する。

点 (a, b, c) = (5, 2, 3) からレーザー光を発射し，反射したレーザー光が点(d, e, f) = (-1, 5, 8) を通るようにするには，x-y 平面のどこを狙えばよいか（(a, b, c) などは，x, y, z の順の座標値）。

from sympy import *
var('a, b, c, d, e, f, x, y, z')

eq1 = Eq((x - a) / (d - a), (z - c) / (f - c))
eq2 = Eq((y - b) / (e - b), (z - c) / (f - c))
eq1 = eq1.subs([(a, 5), (c, -3), (d, -1), (f, 8), (z, 0)])
eq2 = eq2.subs([(b, 2), (c, -3), (e,  5), (f, 8), (z, 0)])
solve([eq1, eq2], [x, y])

数値計算の落とし穴（非常に希なケース？）その3

16677181699666574 >= 3**x を満たす最大の整数 x を求めよ

Python の sympy は，綺麗な解を与える。

>>> from sympy import *
>>> var('a, n, x')
(a, n, x)
>>> a = Ge(n, 3**x)
>>> a
n >= 3**x

>>> solve(a.subs(n, 16677181699666574), x)
x <= log(16677181699666574)/log(3)

>>> solve(a.subs(n, 16677181699666574), x).evalf(30)
x <= 34.0000000000000002728995950932

>>> solve(a.subs(n, 16677181699666573), x).evalf(30)
x <= 34.0000000000000002183196760745

>>> solve(a.subs(n, 16677181699666572), x).evalf(30)
x <= 34.0000000000000001637397570559

>>> solve(a.subs(n, 16677181699666571), x).evalf(30)
x <= 34.0000000000000001091598380373

>>> solve(a.subs(n, 16677181699666570), x).evalf(30)
x <= 34.0000000000000000545799190186

>>> solve(a.subs(n, 16677181699666569), x).evalf(30) # ここが分かれ道
x <= 34.0

>>> solve(a.subs(n, 16677181699666568), x).evalf(30)
x <= 33.9999999999999999454200809814

二次方程式の解を求める - R の場合

「二次方程式の解を求めるプログラムのレビュー」で，実数解の場合に解の公式だけで解を求めてはならないと述べた。

前の記事で，Python の sympy がどうやっているのかレビューした結果を示した。

今回は，R において，sympy 相当の解を求める方法を示す。

sympy の Rational() に相当するのは，gmp ライブラリにおける，bigq クラスの小数である。しかし，許されている演算は，四則演算，比較演算，abs() 程度である。四則演算の中には ^ もあるが，平方根を取る '^0.5' はインプリメントされていない。そこで，以下のような Sqrt() 関数を自分で定義して使わなければならない。

> library(gmp)

> # ニュートン・ラフソン法で bigq クラスの小数の平方根を求める関数
> Sqrt = function(n) {
+   n = as.bigq(n)
+   x = as.bigq(1)
+   EPSILON = as.bigq(x, 1e20)
+   while (TRUE) {
+     x2 = (x + n / x) / 2
+     if (abs(x - x2) <= EPSILON) {
+       return(x2)
+     }
+     x = x2
+   }
+ }


ちゃんと動くか試しておこう
	
> a = Sqrt(1000)
> print(as.numeric(a*a))
[1] 1000

x^2 + 1.0000000000000001*x + 0.0000000000000001 = 0 を解く

> b = as.bigq("10000000000000001", "10000000000000000")
> c = as.bigq( "1", "10000000000000000")	
> x1 = (-b - Sqrt(b^2 - 4*c)) / 2
> print(x1)
Big Rational ('bigq') :
[1] -159999999999999984000000000000000800000000000000000000000000000001/159999999999999984000000000000000800000000000000000000000000000000
> print(as.numeric(x1))
[1] -1
	
> x2 = (-b + Sqrt(b^2 - 4*c)) / 2
> print(x2)
Big Rational ('bigq') :
[1] -15999999999999998400000000000000079999999999999999/159999999999999984000000000000000800000000000000000000000000000000
> print(as.numeric(x2))
[1] -0.0000000000000001
	
> x3 = c / x1
> print(x3)
Big Rational ('bigq') :
[1] -15999999999999998400000000000000080000000000000000/159999999999999984000000000000000800000000000000000000000000000001
> print(as.numeric(x3))
[1] -0.0000000000000001

解の公式だけで両方の実数解を正しく求めることができた。 当然ではあるが，大きい方の実数解を解の公式で求めて，小さい方の実数解は解と係数の関係から求めるというやり方でも正しく求める事が出来る。

二次方程式の解を求める - Python の sympy の場合

「二次方程式の解を求めるプログラムのレビュー」で，実数解の場合に解の公式だけで解を求めてはならないと述べたが，Python の sympy がどうやっているのかレビューした結果を示す。

以下の二次方程式を解け。
x**2 +1.0000000000000001*x + 0.0000000000000001 = 0

結論から述べると，正しいやり方は「二次方程式の係数を Rational() で与え，solve() で解を求める」ということである。

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

係数を Rational() で定義する。

>>> b = Rational(10000000000000001, 10000000000000000)
>>> c = Rational( 1, 10000000000000000)
>>> print(b, c)
10000000000000001/10000000000000000 1/10000000000000000

Rational() ならば精度は完全に保証される。

>>> print(b.evalf(30))
1.00000000000000010000000000000
>>> print(c.evalf(30))
1.00000000000000000000000000000e-16

solve() で解を求める

>>> ans = solve(x**2 + b * x + c)

>>> print(ans)
[-1, -1/10000000000000000]

つまり，

>>> print(ans[0].evalf(30))
-1.00000000000000000000000000000
>>> print(ans[1].evalf(30))
-1.00000000000000000000000000000e-16

以下は読まなくて良い。

ダメなやり方

python の float を Float() で変換しても，evalf() で有効数字を 30 桁表示させれば精度が足りないのは明らかである。

>>> b = Float(1.0000000000000001)
>>> c = Float(0.0000000000000001)
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000000000000000000 9.99999999999999979097786724035e-17

文字型で指定してやれば，精度は増す。

>>> b = Float('1.00000000000000010000000000000000000000000000')
>>> c = Float('0.00000000000000010000000000000000000000000000')
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000010000000000000 1.00000000000000000000000000000e-16

十分（？）な精度の b, c で二次方程式を解くと以下のようになる。

>>> ans = solve(x**2 + b * x + c, x)

デフォルトでの表示は正確な解が得られたと見まがう。

>>> print(ans)
[-1.00000000000000, -1.00000000000000e-16]

しかし，個別の解を .evalf(30) で表示させると，精度が足りないことが明白になる。

>>> print(ans[0].evalf(30))
-0.999999999999999888977697537484
>>> print(ans[1].evalf(30))
-1.00000000000000010235730316482e-16

解の公式だけで実数解を求める場合には，大小どちらの解も精度は十分であることがわかる。（何故 solve() の場合と違うの か？）

>>> x1 = (-b - sqrt(b**2 - 4 * c)) / 2
>>> print(x1.evalf(30))
-1.00000000000000000000000000000
>>> x2 = (-b + sqrt(b**2 - 4 * c)) / 2
>>> print(x2.evalf(30))
-1.00000000000000000000000000000e-16

数学問題-8

Python では sympy を使う。

数式の簡約において，simplify, radsimp, ratsimp などで，式の見た目がけっこう変わるので，まとめておく。

元の式
>>> a = exp(x) / (exp(2*x) + 1)

不定積分を求める
>>> b1 = integrate(a)
>>> b1
RootSum(4*_z**2 + 1, Lambda(_i, _i*log(2*_i + exp(x))))

[0, 1] で定積分を求める
>>> integrate(a, (x, 0, 1)).evalf()
0.432884741619829

見た目は全然違うが，diff() で元に戻ることを確認できる
>>> diff(b1) # --> a
exp(x)/(exp(2*x) + 1)

積分の結果 b1 を simplify()で簡約化する --> b2
>>> b2 = simplify(b1)
>>> b2
-I*log(exp(x) - I)/2 + I*log(exp(x) + I)/2

b2 は b1 と見た目が違う。diff(c) でも元の式 a と見た目が違う。
>>> c = diff(b2)
>>> c
I*exp(x)/(2*(exp(x) + I)) - I*exp(x)/(2*(exp(x) - I))

分数の和で表された式を一つの分数にまとめる ratsimp() で元の式 a に戻ることが確認できる。
>>> ratsimp(c)
exp(x)/(exp(2*x) + 1)

[0, 1] で定積分を求める
>>> integrate(ratsimp(c), (x, 0, 1)).evalf()
0.432884741619829

simplify() では別の形になる
>>> simplify(c)
1/(2*cosh(x))

式の形は違うが，定積分は同じになることが確認できる
>>> integrate(simplify(c), (x, 0, 1)).evalf()
0.432884741619829

分母を有理化する radsimp() もまた別の形になる
>>> radsimp(c)
(2*exp(3*x) + 2*exp(x))/(2*(exp(2*x) + 1)**2)

同じく式の形は違うが，定積分は同じになることが確認できる
>>> integrate(radsimp(c), (x, 0, 1)).evalf()
0.432884741619829

数学問題-7

Python では sympy を使う
WolframAlpha の「高等学校 数学」の例題を解いてみる

https://ja.wolframalpha.com/examples/mathematics/koukousugaku/math-iii/

問題 1　不定積分　∫(1/(x(log x)^2))dx

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

>>> integrate(1/(x*log(x)**2))
-1/log(x)

問題 2　定積分　∫(0からlog3まで) dx/(e^x+5e^(-x)-2)

>>> a = 1/(exp(x) + 5*exp(-x) - 2)
>>> b1 = integrate(a)
>>> b1
RootSum(16*_z**2 + 1, Lambda(_i, _i*log(8*_i + exp(x) - 1)))
>>> simplify(b1)
-I*log(exp(x) - 1 - 2*I)/4 + I*log(exp(x) - 1 + 2*I)/4

>>> integrate(a, (x, 0, log(3)))
RootSum(16*_z**2 + 1, Lambda(_i, _i*log(2/15 - 8*_i/5))) - RootSum(16*_z**2 + 1, Lambda(_i, _i*log(4/5 - 8*_i/5)))
>>> integrate(a, (x, 0, log(3))).evalf()
0.392699081698724
>>> (b1.subs(x, log(3)) - b1.subs(x, 0)).evalf()
0.392699081698724

sympy では得られないが，integrate(a) の別の形

>>> b2 = -atan((1-exp(x))/2)/2
>>> # radsimp(diff(b2)) = radsimp(a) ゆえ，integrate(diff(b2)) = b2 = integrate(a)
>>> radsimp(diff(b2))
exp(x)/(exp(2*x) - 2*exp(x) + 5)
>>> radsimp(a)
exp(x)/(exp(2*x) - 2*exp(x) + 5)
>>> b2.subs(x, log(3)).evalf()
0.392699081698724
>>> (b2.subs(x, log(3)) - b2.subs(x, 0)).evalf()
0.392699081698724

数学問題-6

Python では sympy を使う
WolframAlpha の「高等学校 数学」の例題を解いてみる

https://ja.wolframalpha.com/examples/mathematics/koukousugaku/math-iii/

問題　n が無限大に近付くときの √(n)(√ (n+2)-√n) の極限

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

>>> limit(sqrt(n) * (sqrt(n + 2) - sqrt(n)), n, oo) # 'oo' means ∞ i.e. inf
1

数学問題-5

Python では sympy を使う

WolframAlpha の「高等学校 数学」の例題を解いてみる

https://ja.wolframalpha.com/examples/mathematics/koukousugaku/math-iii/

問題　z^6+27=0 を解く

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

>>> f = z**6 + 27
>>> solve(f)
[-sqrt(3)*I, sqrt(3)*I, -3/2 - sqrt(3)*I/2, -3/2 + sqrt(3)*I/2, 3/2 - sqrt(3)*I/2, 3/2 + sqrt(3)*I/2]

因数分解して，一次式，二次式を解くと以下のようになる

>>> factor(f)
(z**2 + 3)*(z**2 - 3*z + 3)*(z**2 + 3*z + 3)
>>> solve(z**2-3*z+3)
[3/2 - sqrt(3)*I/2, 3/2 + sqrt(3)*I/2]
>>> solve(z**2+3*z+3)
[-3/2 - sqrt(3)*I/2, -3/2 + sqrt(3)*I/2]

数学問題-4

Python では sympy を使う

WolframAlpha の「高等学校 数学」の例題を解いてみる

https://ja.wolframalpha.com/examples/mathematics/koukousugaku/math-i/

問題　sin x + cos x = 1　を解く

>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

>>> f = sin(x) + cos(x) - 1
>>> ans = solve(f)
>>> ans
[0, pi/2]
>>> f.subs(x, ans[0]).evalf()
0
>>> f.subs(x, ans[1]).evalf()
0

以下の図で確認

>>> from sympy.plotting import plot
>>> plot(f, (x, -0.1, 1.7))

数学問題-二重根号の簡約（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)
	

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

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

使用例

> double.root("sqrt(6-sqrt(32))")
$equation
[1] "2 - sqrt(2)"

$expression
expression(2 - sqrt(2))

$eval
[1] 0.5857864

> double.root("sqrt(6-2*sqrt(8))")
$equation
[1] "2 - sqrt(2)"

$expression
expression(2 - sqrt(2))

$eval
[1] 0.5857864

> double.root("sqrt(106+2*sqrt(1288))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+sqrt(5152))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+4*sqrt(322))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+sqrt(368))")
[1] "can't simplify"
> double.root("sqrt(-16+sqrt(368))")
[1] "can't simplify"
> double.root("sqrt(16-sqrt(-368))")
[1] "square root of negative value is not allowed (inner sqrt)"
> double.root("sqrt(16+sqrt(-368))")
[1] "square root of negative value is not allowed (inner sqrt)"
> double.root("sqrt(16-sqrt(368))")
[1] "square root of negative value is not allowed (outer sqrt)"
> double.root("sqrt(18+2*sqrt(81))")
$equation
[1] "6"

$expression
expression(6)

$eval
[1] 6

> double.root("sqrt(13+sqrt(144))")
$equation
[1] "5"

$expression
expression(5)

$eval
[1] 5

プログラム

double.root = function(equation) {
  solve = function(a.plus.b, a.mult.b) {
    for (b in 1:sqrt(a.mult.b)) {
      if (a.mult.b %% b == 0 && a.plus.b == b + a.mult.b %/% b) {
        return(list(a = a.mult.b %/% b, b = b))
      }
    }
    return(NA)
  }
  simplify = function(a) {
    int.root = trunc(sqrt(a))
    if (int.root ^ 2 == a) {
      return(sprintf("%d", int.root))
    } else {
      return(sprintf("sqrt(%d)", a))
    }
  }
  inner = sub("\\)$", "", sub("sqrt\\(", "", equation))
  if (grepl("sqrt\\(-[0-9]*\\)", equation)) {
    return("square root of negative value is not allowed (inner sqrt)")
  } else if (eval(parse(text = inner)) < 0) {
    return("square root of negative value is not allowed (outer sqrt)")
  }
  result1 = eval(parse(text = equation))
  sign = ifelse(grepl("-", equation), "-", "+")
  equation = gsub("sqrt\\(", "", equation)
  equation = gsub("\\)", "", equation)
  # print(equation)
  parsed = as.integer(unlist(strsplit(equation, "[-+*]")))
  # print(parsed)
  if (is.na(parsed[1])) {
    return("can't simplify")
  }
  if (length(parsed) == 3) {
    parsed[2] = parsed[2] ^ 2 * parsed[3]
  }
  parsed[2] = parsed[2] / 4
  if (parsed[2] != trunc(parsed[2])) {
    return("can't simplify")
  }
  # cat(parsed[1], parsed[2], "\n")
  res = solve(parsed[1], parsed[2])
  if (length(res) != 2) {
    return("can't simplify")
  }
  res.str = paste(simplify(res$a), sign, simplify(res$b))
  result2 = eval(parse(text = res.str))
  if (all.equal(result1, result2)) {
      if (sign == "+") {
        res.str = as.character(sqrt(res$a) + sqrt(res$b))
      } else {
        res.str = as.character(sqrt(res$a) - sqrt(res$b))
      }
    result = parse(text = res.str)
    return(list(
      equation = res.str,
      expression = result,
      eval = eval(result)
    ))
  }
  return(sprintf("error: original = %g, simplified = %g", result1, result2))
}

数学問題-3

Python では sympy を使う
WolframAlpha の「高等学校 数学」の例題を解いてみる
https://ja.wolframalpha.com/examples/mathematics/koukousugaku/

 

二重根号の簡約は，直接的に解を求める事が出来ない（simplify では解が得られない）。
そこで，解法の基礎に帰って，それを sympy で解く。

 

問題　sqrt(124 - 6*sqrt(291)) を簡約する

 

>>> from sympy import *
>>> import math

>>> var('a, b, x, y')
(a, b, x, y)

>>> x = 124
>>> y = 10476 # = 36*291 = 4*2619
>>> simplify(sqrt(x-sqrt(y))) # math.sqrt(x - math.sqrt(y))
sqrt(124 - 6*sqrt(291))
>>> simplify(sqrt(x-2*sqrt(y//4)))
sqrt(124 - 6*sqrt(291))

>>> simplify(sqrt(a+b-2*sqrt(a*b))) # この形にしたい
sqrt(a + b - 2*sqrt(a*b)) # == sqrt(b) - sqrt(a), b > a
>>> eq1 = Eq(a+b, x)
>>> eq2 = Eq(a*b, y//4)
>>> res = solve([eq1, eq2]) # 連立方程式を解く
>>> res
[{a: 27, b: 97}, {a: 97, b: 27}]　# a = 27, b = 97 が解
>>> a2, b2 = list(res[0].values())
>>> a2+b2
124
>>> a2*b2
2619
>>> ans0 = sqrt(a2+b2-2*sqrt(a2*b2)) # math.sqrt(124 - 2*math.sqrt(2619))

>>> ans1 = sqrt(max(a2, b2)) - sqrt(min(a2, b2))　# math.sqrt(97) - math.sqrt(27)
>>> ans0 == ans1
False # 等しくはないが
>>> ans1.evalf(), ans0.evalf()
(4.65270537908947, 4.65270537908947) # 実質的に等しい
>>> (ans0 - ans1).evalf()
0.e-121

数学問題-2

Python では sympy を使う
WolframAlpha の「高等学校 数学」の例題を解いてみる
https://ja.wolframalpha.com/examples/mathematics/koukousugaku/

問題　f(x)=log(x)/x^2の最大値 を求める

>>> from sympy import *
>>> from sympy.plotting import plot
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

>>> func = log(x)/x**2 # 関数
>>> func
log(x)/x**2
>>> func2 = diff(func) # 導関数
>>> func2
-2*log(x)/x**3 + x**(-3)
>>> ans = solve(func2) # 導関数 = 0 を解く
>>> ans
[exp(1/2)]
>>> ans[0].evalf()
1.64872127070013
>>> func2.subs(x, ans[0])
0
>>> max = func.subs(x, ans[0]) # 最大値
>>> max
exp(-1)/2
>>> max.evalf()
0.183939720585721

図を描いて確かめておく

>>> plot(func, (x, 1.3, 2))