#!/usr/bin/env python
# coding: utf-8
# # Table of Contents
#
#
#
#
# 式変形(equation manipulation)
#
#
#
# file:/~/python/doing_math_with_python/equation_manipulation.ipynb
#
# cc by Shigeto R. Nishitani 2017
#
# # 数式処理
#
# 数式の変形は，手で直すほうが圧倒的に早くきれいになる場合が多い．しかし，テイラー展開や，複雑な積分公式，三角関数とexp関数の変換などの手間がかかるところを，数式処理ソフトは間違いなく変形してくれる．ここで示すコマンドを全て覚える必要は全くない．というか忘れるもの．ここでは，できるだけコンパクトにまとめて，悩んだときに参照できるようにする．初めての人は，ざっと眺めた後，鉄則からじっくりフォローせよ．
#
# # 式変形に関連したコマンド
# ## 簡単化(simplify)
#
# In[1]:
from sympy import *
x,y=symbols('x y')
simplify(3*x+4*x+2*y)
# In[5]:
exp1=3*sin(x)**3-sin(x)*cos(x)**2
simplify(exp1)
# ## 展開(expand)
#
# In[6]:
expand((x+y)**2)
# ## 因数分解(factor)
#
# In[7]:
factor(4*x**2-6*x*y+2*y**2)
# ## 約分・通分(normal)
#
# In[11]:
simplify((x+y)/(x**2-3*x*y-4*y**2))
# In[19]:
together(1/x+1/y)
# In[9]:
get_ipython().run_line_magic('pinfo', 'sqrtdenest')
# ## 項を変数でまとめる(collect)
#
# In[18]:
a,b,c = symbols('a b c')
collect(4*a*x**2-3*y**2/x+6*b*x*y+3*c*y+2*y**2,y)
# ## 分母，分子(denom,numer)
#
# 分数を英語でfractionというが，これでとり出せる．出力は，numerator(分子)，denominator(分母)の順．
# In[21]:
a,b,c, x = symbols('a b c x')
frac = (1+x)**2/(a*(x-1)+b*(x+3))
n,d = fraction(frac)
pprint(n)
pprint(d)
# ## 係数(coeff)
#
# In[25]:
eq = a+x+b*x+c*x**2
pprint(eq)
eq.coeff(x)
# In[34]:
eq.as_coefficients_dict()
# In[49]:
eq0 = Poly(eq, x)
coeffs = eq0.coeffs()
print(coeffs[::-1])
# ## 分数をそのまま使いたい
#
# 方程式などで，係数として分数を使うとうまく扱ってくれません．
# そんなときは，Rationalで有理数と指定してしまいます．
# In[7]:
factor(1/4*x**2-y**2)
# In[8]:
factor(Rational(1,4)*x**2-y**2)
# ## 課題
# ### 因数分解
# 1. $x^3-8$
# 1. $abc−abx−acx+ax^2+bcx−bx^2−cx^2+x^3$
#
# ### 部分分数
# 1. $ \frac{x+1}{(x-1)(x^2+1)^2}$
# 1. $ \frac{3}{1-x^4} = \frac{a}{x^2+1}+\frac{b}{x+1}+\frac{c}{x-1}$が常に成立する$a, b, c$を定めよ．
#
# ### 分母の有理化
# $ \frac{8}{3-\sqrt{5}}-\frac{2}{2+\sqrt{5}}$を簡単化せよ．
#
# ### 複合問題
# $ x^2+2kx+5-k=0$が重根をもつように$k$を定めよ．
#
# ### 解答例
# In[8]:
factor(x**3-8)
# In[14]:
c = symbols('c')
factor(a*b*c - a*b*x - a*c*x + a*x**2 + b*c*x - b*x**2 - c*x**2 + x**3)
# In[17]:
apart((x+1)/((x-1)*(x**2+1)**2))
# In[21]:
apart(3/(1-x**4))
# In[5]:
from sympy import *
radsimp(8/(3-sqrt(5)))
radsimp(2/(2+sqrt(5)))
simplify(8/(3-sqrt(5))-2/(2+sqrt(5)))
# In[12]:
k,x = symbols('k x')
eq=x**2 +2*k*x+5-k
sol = solve(eq,x)
solve(sol[0]-sol[1],k)
# # 式の変形の基本(BottomLine)
#
#
# どうしても解かなければならない課題を前にコマンドリファレンスのあちこちを参照しながら解いていくのが数式処理を修得する最速法である．とびかかる前にちょっとした共通
# のコツがある．それをここでは示す．数式処理ソフトでの数式処理とは，数式処理ソフトが『自動的にやって』くれるのではなく，実際に紙と鉛筆で解いていく手順を数式処理ソ
# フトに『やらせる』ことであることを肝に銘じよ．
#
# ## 鉄則
#
# 数式処理ソフトの習得にあたって初心者がつまづく共通の過ちを回避する鉄則がある．
#
#
# - 鉄則０：restart をかける
#
# 続けて入力すると前の入力が生きている．違う問題へ移るときや，もう一度入力をし直すときには，Kernel->Restartを押して初期状態からはじめる．入力した順番が狂っている場合もある．頭から順にshift-returnをやり直す．
# さらに，
#
# ``` python
# # kernel->Restart
# from sympy import *
# init_session()
# ```
#
# として，sympyのimport，init_sessionによる綺麗な表示を強制しておく．
#
# - 鉄則１：出力してみる
#
# テキストではページ数の関係で出力を抑止しているが，初心者が問題を解いていく段階ではデータやグラフをできるだけ多く出力する．途中の結果を明示するためにpprint文を入れる．
#
# - 鉄則２：関数に値を代入してみる
#
# 数値が返ってくるべき時に変数があればどこかで入力をミスっている．plotで以下のようなエラーが出た場合にチェック．
#
# ``` python
# ValueError: The same variable should be used in
# all univariate expressions being plotted.
# ```
#
#
# - 鉄則３：内側から順に入力する
#
# 長い入力やfor-loopを頭から打ち込んではいけない！！ 内側から順に何をしているか解読・確認しながら打ち込む．括弧が合わなかったり，読み飛ばしていたりというエラーが回避できる．
#
#
#
# ## 数式変形の実践（大学入試センター試験の解法を通して）
#
# 今まで出てきたコマンドを使えば，典型的なセンター試験の問題を解くのも容易である．以下の例題を参照して課題を解いてみよ．使うコマンドは，unapply, solve, diff, expand(展開), factor(因数分解)とsubs(一時的代入)である．
#
#
# ### 数式処理問題の典型例:2次関数の頂点
#
# $a,b$を定数とし, $a \neq 0$とする．2次関数
#
# $$
# y = a\,x^2-b\,x-a+b\,\cdots(1)
# $$
# のグラフが点 $(-2, 6)$ を通るとする．このとき
#
# $$
# b = -a+\fbox{ ア }
# $$
# であり，グラフの頂点の座標を$a$を用いて表すと
#
# $$
# \left(\frac{-a+\fbox{ イ }}{\fbox{ ウ }\,a}, -\frac{(\fbox{ エ }\,a- \fbox{ オ })^2}{\fbox{ カ }\,a}\right)
# $$
# である (2008 年度大学入試センター試験数学 I 第2問より抜粋)．
#
#
#
# ### 解答例
#
# In[1]:
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b
bb = solve(f.subs({x:-2})-6,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
y0=factor(f.subs({x:x0,b:bb}))
pprint(y0)
# 解答の一例は以上のとおりであるが，これはだいぶ整理されたあとの記述になる．
# 普通の答案ではここまで綺麗に初めから求まるものではない．
#
# そこで，綺麗な解答を導く一つ一つのステップをはしょることなくここで示しておく．
# 今後巡り会うであろう解答例やコードサンプルにおいても，
# ここで示した導出のコンセプトを参考に数式処理コードを理解していってもらいたい．
#
# まずはRestartをかける．これもcodeに示しておくと忘れない．
# それとsympyのimportとinit_sessionによる綺麗なprint outを癖付けておく．
# In[1]:
# kernel->Restart
from sympy import *
init_session()
# 与えられた2次関数を$f(x)$で定義する
# In[2]:
f = a*x**2-b*x-a+b
# すると早速エラー．sympyでは変数を明示的に宣言しないといけない．
# In[3]:
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
# 関数fが点(−2,6)を通る式を立てる．まずはx=-2を代入(subs)
# In[4]:
f.subs({x:-2})
# こいつがy=-6をみたすのだから，yに前の式を代入して，左辺にまとめて見る．
# In[5]:
y = f.subs({x:-2})
y-6
# これをeq1と置いて，bについて解く．
# In[6]:
eq1 = y-6
solve(eq1,b)
# これがbの表式．配列の一番目の要素となっているので，それを取り出しておく．
# In[7]:
bb = solve(eq1,b)[0]
pprint(bb)
# 頂点の座標(x0,y0)で傾きが0になることを用いて解いていく．鉄則3を思い出して，中から見ていく．
# 微分diffを思い出して，
# In[8]:
diff(f,x)
# これをxについて解く．
# In[9]:
solve(diff(f,x),x)
# In[10]:
solve(diff(f,x),x).subs({b:bb})
# oops. このエラーはlistに対してはsubsは使えないとのこと．
# 先に中身を取り出しておいて，代入．
# In[11]:
solve(diff(f,x),x)[0].subs({b:bb})
# これが頂点x0．これをfに代入してやる．ついでにb=bbも．．．
# In[12]:
f.subs({x:x0,b:bb})
# せやった．x0に代入してない．
# In[13]:
x0=solve(diff(f,x),x)[0].subs({b:bb})
f.subs({x:x0,b:bb})
# こいつを簡単化simplifyすればOK.
# In[14]:
simplify(f.subs({x:x0,b:bb}))
# あかん．ここは明示的に因数分解factorを指定する．
# In[15]:
factor(f.subs({x:x0,b:bb}))
# でや？！ これをまとめて表示すると，
# In[16]:
# kernel->Restart
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
y = f.subs({x:-2})
eq1 = y-6
bb = solve(eq1,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
factor(f.subs({x:x0,b:bb}))
# あれれ．．．なんか途中が抜けとるで？？
#
# 途中の重要な結果は変数に置き換えて明示pprintしておくのが鉄則１．
# それを意識して修正すると，
# In[2]:
# kernel->Restart
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
bb = solve(f.subs({x:-2})-6,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
y0=factor(f.subs({x:x0,b:bb}))
pprint(y0)
# ## 課題
# ### 課題1
# （例題1に続いて)
#
# さらに，2次関数(1)のグラフの頂点の$y$座標が-2であるとする．このとき，$a$は
# \begin{equation*}
# \fbox{ キ }\,a^2-\fbox{ クケ }\,a+\fbox{ コ } = 0
# \end{equation*}
# を満たす．これより，$a$の値は
# \begin{equation*}
# a = \fbox{ サ }, \frac{\fbox{ シ }}{\fbox{ ス }}
# \end{equation*}
# である．
# 以下，$a = \frac{\fbox{ シ }}{\fbox{ ス }}$であるとする．
#
# このとき，2次関数(1)のグラフの頂点のx座標は$\fbox{ セ }$であり，(1)のグラフと$x$軸の2交点の$x$座標は$\fbox{ ソ },\fbox{ タ }$である．
#
# また，関数(1)は$0 \leqq x \leqq 9$において
#
# $x = \fbox{ チ }$のとき，最小値$\fbox{ ツテ }$をとり，
#
# $x = \fbox{ ト }$のとき，最大値$\frac{\fbox{ ナニ }}{\fbox{ ヌ }}$をとる．
#
#
# (2008 年度大学入試センター試験数学 I より抜粋).
# #### 解答例
# In[4]:
# 実践の解答例に続けて．．．
eq2 = expand((y0 +2)*(-4*a)) # -4aをかけて，展開
pprint(eq2) # キクケコ
s_a=solve(eq2,a) # 方程式を解く
pprint(s_a) # サシス
aa = min(s_a) #二つの解の小さい方が題意，シス
pprint(aa)
pprint(x0.subs({a:aa})) # セ
eq3 = f.subs({b:bb}).subs({a:aa}) # fに得られたa,bを代入
pprint(solve(eq3,x)) # そいつを解いてx軸との交点を，ソタ
# In[5]:
get_ipython().run_line_magic('matplotlib', 'inline')
plot(eq3, (x,0,9)) # 得られたeq3をplot
# In[6]:
print(eq3.subs({x:4})) # min at x=4
print(eq3.subs({x:9})) # 目視でmax at x=9
# ### 課題2
# $P = x(x+3)(2x-3)$とする．
# また，$a$を定数とする．
# $x = a+1$のときの
# $P$の値は
# \begin{equation*}
# 2a^3+\fbox{ ア }a^2+\fbox{ イ }a-\fbox{ ウ }
# \end{equation*}
# である．
#
# $x=a+1$のときの$P$の値と，$x=a$のときの$P$の値が等しいとする．このとき，$a$は
# \begin{equation*}
# 3a^2+\fbox{ エ }a-\fbox{ オ } = 0
# \end{equation*}
# を満たす．したがって
# \begin{equation*}
# a = \frac{\fbox{ カキ }\pm \sqrt{\fbox{ クケ }}}{\fbox{ コ }}
# \end{equation*}
# である.
#
# (2008 年度大学入試センター試験数学 I 第4問より抜粋)
# #### 解答例
# In[21]:
from sympy import *
a,b,x = symbols('a b x')
P = x*(x+3)*(2*x-3) # $P$を$x$の関数として定義
pprint(expand(P.subs({x:a+1}))) # $$P(a)$を形式的に出してみる．
eq1 = expand(P.subs({x:a+1}) - P.subs({x:a})) #2式を差し引く
pprint(eq1/2) # 出題にそろえるため2で割っている．
solve(eq1,a) # eq1=0をaについて解く(solve)
# # 追記(ちょっとしたコツ)
# ## 式フォローのデフォルト
#
# 数式処理ソフトで実際に数式をいじる状況というのは，ほとんどの場合が既知の数式変形のフォローだろう．例えば，論文で「(1)式から(2)式への変形は自明である」とかいう
# 文章で済ましている変形が本当にあっているのかを確かめたい時．一番単純なやり方は自明と言われた前後の式が一致していることを確かめるだけで十分である．
# 最も単純な確認法は以下の通り，変形の前後の式を手入力してその差をexpandした結果が0か否かでする．
#
# In[19]:
from sympy import *
x = symbols('x')
ex1=(x-3)**4
ex2 = x**4-12*x**3+54*x**2-108*x+81
expand(ex1-ex2)
# 0ならば式の変形は保証されているので，その導出が間違いでなく誤植などもないことが確認できる．ただ，これだけでは変形の哲学や技法が身に付くわけではない．あくまでも
# 苦し紛れのデフォルトであることは心に留めておくように．
# 論理値(True or False)として確かめたいときには，==0で出てくる．
# ## assumeの具体例：無限積分
# 以下に示す積分を実行せよ．
#
# $$
# \int _{-\infty }^{\infty }x \exp(-\beta c\,{x}^{2}) \left( 1+\beta g\,{x}^{3} \right) {dx}
# $$
# 最新版のMapleでは改良が施されていて，このような複雑な積分も一発で求まる．
#
# ```maple
# > int(f1(x),x=-infinity..infinity);
# ```
# $$
# \left\{\, \begin {array}{cc} { \frac{3}{4}\,\frac {g \sqrt{\pi }}{\beta\,{c}^{2} \sqrt{c\beta}}}&csgn \left( c\beta \right) =1\\ \infty&otherwise\end {array} \right.
# $$
# ここでは，$c \beta$が正の場合(csgn(beta c)=1)とそれ以外の場合(otherwise)に分けて答えを返している．しかしこのような意図したきれいな結果をいつも数式処理ソフトが返してくれるとは限らない．これだけしか知らないと，なにかうまくいかないときにお手上げになってしまう．
# In[1]:
from sympy import *
init_printing()
x,beta,c,g = symbols('x beta c g')
f1 = x*exp(-beta*c*x**2)*(1+beta*g*x**3)
pprint(f1)
pprint(integrate(f1,(x,-oo,oo)))
# この辺りの記述をやり直したけど，
# assumeのためにassumpitons moduleを使ってやったんだけど，
# どうもうまくいかない．
# ``` python
# assumes = Q.positive(beta), Q.positive(c)
# with assuming(*assumes):
# print ask(beta*c, positive)
# ```
# とかは通るけど，integrateの結果は変わらない．
# 結局，symbols('x',positive = True)とかで初めに指定するとすっと通る．
#
#
# In[ ]: