#!/usr/bin/env python
# coding: utf-8

# # Symbolic Expressions and Simplification

# In[ ]:


# Sage does not recognise two expressions as equal 
bool(arctan(1+abs(x)) == pi/2 - arctan(1/(1+abs(x))))


# In[ ]:


x = var('x')
plot( arctan(1+abs(x)) + arctan(1/(1+abs(x))) , (x, -5, 5))


# In[ ]:


a, x = var('a, x')
y = cos(x+a) * (x+1); y


# In[ ]:


y.subs(a=-x)


# In[ ]:


y.subs(x=pi/2, a=pi/3)


# In[ ]:


y.subs(x=0.5, a=2.3)


# In[ ]:


y(a=-x)


# In[ ]:


y(x=pi/2, a=pi/3)


# In[ ]:


y(x=0.5, a=2.3)


# In[ ]:


x, y, z = var('x, y, z')
q = x*y + y*z + z*x
bool(q(x=y, y=z, z=x) == q), bool(q(z=y)(y=x) == 3*x^2)


# In[ ]:


y, z = var('y, z'); f = x^3 + y^2 + z
f.substitute(x^3 == y^2, z==1)


# In[ ]:


f(x)=(2*x+1)^3 
f(-3)


# In[ ]:


f.expand()


# In[ ]:


y = var('y')
u = sin(x) + x*cos(y)
v = u.function(x, y)
v


# In[ ]:


w(x, y) = u
w


# In[ ]:


x, y = SR.var('x,y')


# In[ ]:


p = (x+y)*(x+1)^2
p2 = p.expand()
p2


# In[ ]:


p2.collect(x)


# In[ ]:


((x+y+sin(x))^2).expand().collect(sin(x))


# In[ ]:


(x^x/x).simplify()


# In[ ]:


f = (e^x-1) / (1+e^(x/2))
f.canonicalize_radical()


# In[ ]:


f = cos(x)^6 + sin(x)^6 + 3 * sin(x)^2 * cos(x)^2
f.simplify_trig()


# In[ ]:


f = cos(x)^6
f.reduce_trig().show()


# In[ ]:


f = sin(5 * x)
f.expand_trig().show()


# In[ ]:


n = var('n')
f = factorial(n+1)/factorial(n)
f.simplify_factorial()


# In[ ]:


f = sqrt(abs(x)^2)
f.canonicalize_radical()


# In[ ]:


f = log(x*y)
f.canonicalize_radical()


# In[ ]:


forget()
assume(x > 0)
bool(sqrt(x^2) == x)


# In[ ]:


forget(x > 0)
bool(sqrt(x^2) == x)


# In[ ]:


n = var('n')
assume(n, 'integer')
sin(n*pi)


# In[ ]:


a = var('a')
c = (a+1)^2 - (a^2+2*a+1)    # c = 0
eq = c * x == 0
eq


# In[ ]:


eq2 = eq / c     # divide by zero
eq2


# In[ ]:


solve(eq2, x)


# In[ ]:


solve(eq, x)     # Sage avoids this mistake.


# In[ ]:


expand(c)


# In[ ]:


c = cos(a)^2 + sin(a)^2 - 1    # c = 0
eq = c*x == 0
solve(eq, x)      # Sage does not avoid the pitfall.


# In[ ]:


c.simplify_trig()


# In[ ]:


c.is_zero()


# # Equations

# In[ ]:


z, phi = var('z, phi')
eq = z**2 - 2/cos(phi)*z + 5/cos(phi)**2 - 4 == 0
eq.show()


# In[ ]:


eq.lhs()


# In[ ]:


eq.rhs()


# In[ ]:


sols = solve(eq, z)
show(sols)


# In[ ]:


y = var('y')
sols = solve(y^7==y, y)
show(sols)


# In[ ]:


solve(x^2-1, x, solution_dict=True)     # dictionary


# In[ ]:


solve([x+y == 3, 2*x+2*y == 6], x, y)   # System of two equations


# In[ ]:


solve([cos(x)*sin(x) == 1/2, x+y == 0], x, y)


# In[ ]:


sols = solve(x^2+x-1 > 0, x)            # inequalities
show(sols)


# In[ ]:


x, y, z = var('x, y, z')
sols = solve([x^2 * y * z == 18, x * y^3 * z == 24, x * y * z^4 == 6], x, y, z)
for sol in sols: 
    print(sol)
len(sols)        # 17 solutions, among which 16 are approximate complex solutions


# In[ ]:


expr = sin(x) + sin(2 * x) + sin(3 * x)
solve(expr, x)


# In[ ]:


find_root(expr, 0.1, pi)     # numerical solution


# In[ ]:


f = expr.simplify_trig()
f


# In[ ]:


solve(f, x)


# In[ ]:


sols = (x^3+2*x+1).roots(x)    # with their multiplicity
show(sols)


# In[ ]:


(x^3+2*x+1).roots(x, ring=RR)    # RR means "Real Field with 53 bits of precision"


# In[ ]:


(x^3+2*x+1).roots(x, ring=CC)    # CC means "Complex Field with 53 bits of precision"


# In[ ]:


solve(x^(1/x)==(1/x)^x, x)     # no explicit solution


# In[ ]:


y = function('y')(x)
y


# In[ ]:


desolve(diff(y,x,x) + x*diff(y,x) + y == 0, y, [0,0,1]).show()


# # Analysis

# ## Sums

# In[ ]:


k, n = var('k, n')
sum(k, k, 1, n).factor().show()


# In[ ]:


n, k, y = var('n, k, y')
sum(binomial(n,k) * x^k * y^(n-k), k, 0, n).show()


# In[ ]:


k, n = var('k, n')
(
    sum(binomial(n,k), k, 0, n), 
    sum(k * binomial(n, k), k, 0, n), 
    sum((-1)^k*binomial(n,k), k, 0, n),
)


# In[ ]:


a, q, k, n = var('a, q, k, n')
sum(a*q^k, k, 0, n).show()


# In[ ]:


forget()
assume(abs(q) < 1)
sum(a*q^k, k, 0, oo).show()   # oo means infinity


# In[ ]:


forget()
assume(q > 1)
# sum(a*q^k, k, 0, infinity)    # error since the sum is divergent.


# ## Limits

# Let us compute the following limits:
# * $\displaystyle \lim_{x \to 8} \frac{\sqrt[3]{x} -2 }{\sqrt[3]{x+19} - 3} $
# * $\displaystyle \lim_{x \to \frac{\pi}{4}} \frac{\cos\left( \frac{\pi}{4} - x \right) - \tan x}{1 - \sin\left( \frac{\pi}{4} + x \right)} $

# In[ ]:


limit((x**(1/3) - 2) / ((x + 19)**(1/3) - 3), x = 8)


# In[ ]:


f(x) = (cos(pi/4-x)-tan(x))/(1-sin(pi/4 + x))
limit(f(x), x = pi/4)


# In[ ]:


limit(f(x), x = pi/4, dir='minus')


# In[ ]:


limit(f(x), x = pi/4, dir='plus')


# ## Sequences

# **Example.** (A sequence study) Let us consider the sequence 
# $$u_n = \frac{n^{100}}{100^n}.$$
# - Compute the first $10$ terms. 
# - How does the sequence vary?
# - What is the sequence limit? 
# - From which value of $n$ does $u_n \in (0, 10^{-8})$ hold?

# In[ ]:


u(n) = n^100 / 100^n
u(1.), u(2.), u(3.), u(4.), u(5.), u(6.), u(7.), u(8.), u(9.), u(10.)


# In[ ]:


plot(u(x), x, 1, 40)


# In[ ]:


v(x) = diff(u(x), x)
sol = solve(v(x) == 0, x)
sol


# In[ ]:


floor(sol[0].rhs())


# In[ ]:


limit(u(n), n=infinity)


# In[ ]:


n0 = find_root(u(n) - 1e-8 == 0, 22, 1000)
n0


# ## Power Series Expansions

# Let us determine the power series expansion of the following functions:
# * $(1+\arctan x)^{\frac{1}{x}}$ of order $3$, at $x_0 = 0$;
# * $\ln(2\sin x)$ of order $3$, at $x_0 = \frac{\pi}{6}$.

# In[ ]:


((1+arctan(x))^(1/x)).series(x==0, 3).show()


# In[ ]:


(ln(2*sin(x))).series(x==pi/6, 3).show()


# In[ ]:


(ln(2*sin(x))).series(x==pi/6, 3).truncate().show()


# $(x^3 + x)^{\frac{1}{2}} - (x^3 - x)^{\frac{1}{3}}$ behaves around $+\infty$.

# In[ ]:


taylor((x**3+x)**(1/3) - (x**3-x)**(1/3), x, infinity, 10)


# **Example.** (Machin’s formula) Prove the following formula:
# $$\frac{\pi}{4} = 4 \arctan \frac{1}{5} - \arctan \frac{1}{239}$$

# In[ ]:


tan(4*arctan(1/5)).simplify_trig()


# In[ ]:


tan(pi/4+arctan(1/239)).simplify_trig()


# In[ ]:


f = arctan(x).series(x, 10)
f.show()


# In[ ]:


(
    (16*f.subs(x==1/5) - 4*f.subs(x==1/239)).n(),
    pi.n(),
)


# ## Series

# **Example.** (Evaluation of the Riemann zeta function)

# In[ ]:


k = var('k')
expr = (
    sum(1/k^2, k, 1, infinity),
    sum(1/k^4, k, 1, infinity),
    sum(1/k^5, k, 1, infinity),
)
show(expr)


# **Example.** (A formula due to Ramanujan) Using the first $12$ terms of the following series, we give an approximation of and we compare it with the value given by Sage.
# $$\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum_{k=0}^{\infty} \frac{(4k)! \cdot (1103+26390 k)}{(k!)^4 \cdot 396^{4k}}$$

# In[ ]:


s = 2*sqrt(2)/9801 * (
    sum(
        ( factorial(4*k) * (1103+26390*k) ) / ((factorial(k))^4 * 396^(4*k)) 
        for k in (0..11)
       )
)
(1/s).n(digits=100)


# In[ ]:


(pi-1/s).n(digits=100).n()
# the partial sum of the first 12 terms already yields 95 significant digits of pi.


# **Example.** (Convergence of a series) Let us study the convergence of the series
# $$\sum_{n\geq 0 } \sin\left( \pi \sqrt{4n^2 + 1} \right)$$

# In[ ]:


n = var('n')
u = sin(pi*(sqrt(4*n^2+1)-2*n))
taylor(u, n, infinity, 3).show()


# We deduce $u_n \sim \frac{\pi}{4}$, so the series $\sum_{n\geq 0} u_n$ diverges.

# ## Derivatives

# In[ ]:


diff(sin(x^2), x)


# In[ ]:


reset()
f(x) = function('f')(x)
g(x) = function('g')(x)
diff(f(g(x)), x).show()


# In[ ]:


diff(ln(f(x)), x)


# ## Partial Derivatives

# In[ ]:


f(x,y) = x*y + sin(x^2) + e^(-x)
derivative(f, x)


# In[ ]:


derivative(f, y)    # or diff(f, y), or f.diff(y)


# Let us check that the following function is _harmonic_:
# $$f(x,y) = \frac{1}{2} \ln(x^2 + y^2) \quad \text{for all $(x,y) \neq (0,0)$.}$$
# Note that a function $f$ is said _harmonic_ when its Laplacian $\Delta f = \partial_1^2 f + \partial_2^2 f$.

# In[ ]:


x, y = var('x, y')
f = ln(x**2+y**2) / 2
delta = diff(f,x,2) + diff(f,y,2)
delta.simplify_rational()


# ## Integrals

# In[ ]:


sin(x).integral(x, 0, pi/2)


# In[ ]:


integrate(1/(1+x^2), x)


# In[ ]:


integrate(1/(1+x^2), x, -infinity, infinity)


# In[ ]:


integrate(exp(-x**2), x, 0, infinity)


# In[ ]:


# integrate(exp(-x), x, -infinity, infinity)   # error since integral is divergent.


# **Example.** Let us compute, for x 2 R, the integral
# $$\varphi(x) = \int_0^{\infty} \frac{x \cos u }{u^2+x^2} \ du.$$

# In[ ]:


forget()
u = var('u')
f = x * cos(u) / (u^2 + x^2)
# f.integrate(u, 0, infinity)      # error, Is x positive, negative or zero?


# In[ ]:


forget()
assume(x > 0)
A = f.integrate(u, 0, infinity, hold = 'True')
show(A == A.unhold())


# In[ ]:


forget()
assume(x < 0)
A = f.integrate(u, 0, infinity, hold = 'True')
show(A == A.unhold())


# We have: $\quad \varphi(x) = \frac{|x|}{x}  \frac{\pi}{2} e^{-|x|}$ for all $x \in \mathbb{R}^{*}.$

# In[ ]:


integral_numerical(sin(x)/x, 0, 1)    # ( the answer, an error estimate )


# In[ ]:


g = integrate(exp(-x**2), x, 0, infinity)
g, g.n()


# In[ ]:


approx = integral_numerical(exp(-x**2), 0, infinity)
approx


# In[ ]:


approx[0]-g.n()


# # Basic Linear Algebra

# ## Matrix Computations

# In[ ]:


A = matrix(QQ, [[1,2],[3,4]])
A


# In[ ]:


show(A)    # pretty_print(A)


# In[ ]:


B = A.inverse()
show(B)


# In[ ]:


A = matrix(QQ, [[2,4,3],[-4,-6,-3],[3,3,1]])
A.characteristic_polynomial()


# In[ ]:


A.eigenvalues()


# In[ ]:


A.minimal_polynomial().factor()


# In[ ]:


eigen = A.eigenvectors_right()
show(eigen)


# In[ ]:


J = A.jordan_form(transformation=True)  # A is triangularisable.
show(J)


# **Example.** Let us diagonalise the matrix 
# $$A = \begin{pmatrix}
# 1 & −1/2 \\
# −1/2 & −1 \\
# \end{pmatrix}
# .$$

# In[ ]:


A = matrix(QQ, [[1,-1/2],[-1/2,-1]])  
# A.jordan_form()          # error, A is NOT triangularisable in Rational Field.


# In[ ]:


A = matrix(QQ, [[1,-1/2],[-1/2,-1]])
A.minimal_polynomial()


# In[ ]:


R = QQ[sqrt(5)]
A = A.change_ring(R)
J = A.jordan_form(transformation=True, subdivide=False)
show(J)


# **Example.** Let us diagonalise the matrix 
# $$A = \begin{pmatrix}
# 2 & \sqrt{6} & \sqrt{2} \\
# \sqrt{6} & 3 & \sqrt{3} \\
# \sqrt{2} & \sqrt{3} & 1 \\
# \end{pmatrix}
# .$$

# In[ ]:


K.<sqrt2> = NumberField(x^2 - 2)
L.<sqrt3> = K.extension(x^2 - 3)
A = matrix(L, 
           [[2, sqrt2*sqrt3, sqrt2],
            [sqrt2*sqrt3, 3, sqrt3],
            [sqrt2, sqrt3, 1],
           ] 
          )
show(A)


# In[ ]:


J = A.jordan_form(transformation=True)
show(J)


# In[ ]: