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

# # First contact with SageMath
# 
# This Jupyter notebook illustrates a few elementary features of [SageMath](https://www.sagemath.org/).
# 
# First we set up the display to have LaTeX outputs:

# In[1]:


get_ipython().run_line_magic('display', 'latex')


# ## Starting with numbers...
# 
# SageMath knows about $\pi$, $e$ and $i$ (well, it's mathematical software, isn't it ?):

# In[2]:


e^(i*pi) + 1


# *--- It's nice, but I thought SageMath was based on Python: shouldn't the above be written* `e**(i*pi) + 1`, 
# *given that* `^` *is the bitwise XOR operator in Python?*
# 
# Actually, the input cells are **preparsed** by SageMath before being sent to the Python interpreter. The action of the preparser is revealed by the function `preparse`:

# In[3]:


preparse("e^(i*pi) + 1")


# We see that indeed, the character `^` has been changed to `**`.
# 
# Another noticable change is `1` $\to$ `Integer(1)`. This means that the preparser is turning integers into SageMath integers, which belong to the class `Integer`:

# In[4]:


type(1)


# This type of integer is much more sophisticated than a mere Python integer (`int`). In particular, it knows to which mathematical set it belongs. The latter is returned by the function `parent`:

# In[5]:


parent(1)


# To have more information than just the symbol of the parent, one can use the `print` function:

# In[6]:


print(parent(1))


# Let us denote this object by the Python variable `Z`:

# In[7]:


Z = parent(1)
Z


# `Z` is endowed with many methods, which can be discovered via the **TAB key**:
# ```
#     Z.<TAB>
# ```
# Once a method has been selected, one can get some **documentation** about it via `?`:

# In[8]:


get_ipython().run_line_magic('pinfo', 'Z.cardinality')


# A double question mark leads directly to the **source code** (recall that SageMath is a free software!):

# In[9]:


get_ipython().run_line_magic('pinfo2', 'Z.cardinality')


# More documentation is available in the online [reference manual](https://doc.sagemath.org/html/en/reference/) as well as in various [tutorials and guides](https://doc.sagemath.org/html/en/).

# In[10]:


Z.cardinality()


# Other methods:

# In[11]:


Z.is_ring()


# In[12]:


Z.is_field()


# In[13]:


Z.zero(), Z.one()  # the ring zero and unit elements


# Expressing $1\in \mathbb{Z}$ in SageMath:

# In[14]:


1 in Z


# 1 is actually the unit element of the ring $\mathbb{Z}$:

# In[15]:


1 == Z.one()


# Note that there exits a predefined variable, `ZZ`, for $\mathbb{Z}$:

# In[16]:


Z is ZZ


# SageMath integers can be arbitrarily large (up to the limit of the computer memory):

# In[17]:


m = 2^127 - 1
m


# In[18]:


m.is_prime()


# Note that one has still access to the Python integers, via `int`:

# In[19]:


int(1)


# In[20]:


type(_)  # the underscore stands for the latest output


# They are automatically converted (coerced) into SageMath integers if necessary, for instance if they are
# added to some SageMath integer:

# In[21]:


int(1) + 2


# In[22]:


type(_)


# Of course, SageMath also knows about **rational numbers**:

# In[23]:


2/3


# Python would have returned 0.6666666666666666. This is not the case here because

# In[24]:


preparse('2/3')


# and the division of the (Sage) integer 2 by the (Sage) integer 3 is the (Sage) rational number $2/3$: 

# In[25]:


parent(2/3)


# In[26]:


print(_)


# As the above examples illustrate, SageMath is based on a **Parent/Element scheme**: `parent(a)` is the algebraic or topological/differential structure to which `a` belongs. 
# 
# SageMath has also the concept of **mathematical categories**:

# In[27]:


print(category(Z))


# ### Real numbers
# 
# SageMath can compute numerical values with an arbitrary number of digits:

# In[28]:


n(pi, digits=1000)


# `n` is the shortcut alias for `numerical_approx`.

# Another interesting computation regards the *Hermite-Ramanujan constant*:

# In[29]:


a = exp(pi*sqrt(163))
a


# Actually, this number is very close to an integer, as announced by [Charles Hermite in 1859](https://gallica.bnf.fr/ark:/12148/bpt6k6209489m/f60) (probably without using SageMath...):

# In[30]:


n(a, digits=50)


# That's clear if we turn off scientific notation:

# In[31]:


n(a, digits=50).str(no_sci=2)


# ## Symbolic computations
# 
# Beside numerical computations, SageMath can perform symbolic ones. For instance, it can compute a **derivative**:

# In[32]:


f = diff(sin(x^2), x)
f


# Due to the command `%display latex` in the first cell of this notebook, SageMath displays all results, such as the one above, in LaTeX format. To get them in console mode, use the function `print`:

# In[33]:


print(f)


# The explicit LaTeX code can be obtained:

# In[34]:


print(latex(f))


# A pdf file with the LaTeX typeset formula is generated by the function `view`:

# In[35]:


view(f)


# `x` is the only generic symbolic variable predefined in a SageMath session. All other symbolic variables must be declared explicitly, via the function `var`:

# In[36]:


y = var('y')

alp = var('alp', latex_name=r'\alpha')


# In[37]:


f = alp*y + x
f


# In[38]:


diff(f, alp)


# SageMath can also compute a **primitive**:

# In[39]:


integrate(x^5/(x^3 - 2*x +1), x, hold=True) 


# In[40]:


integrate(x^5/(x^3 - 2*x + 1), x) 


# and of course definite integrals: it suffices to provide the two boundaries after the argument `x`. A first attempt with 0 and 1 fails:

# In[41]:


# integrate(x^5/(x^3 - 2*x + 1), x, 0, 1)


# **Don't panic:** error messages are usually long because they display the whole call stack. The important information lies at the end:
# ```
# ValueError: Integral is divergent
# ```

# Choosing the boundaries to be 2 and 3 yields a convergent integral:

# In[42]:


integrate(x^5/(x^3 - 2*x + 1), x, 2, 3)


# In[43]:


n(_)  # numerical approximation of the above result


# SageMath proposes various symbolic engines to evaluate a primitive. The default one is [Maxima](https://maxima.sourceforge.io/):

# In[44]:


integrate(x^5/(x^3 - 2*x + 1), x, algorithm='maxima')  # same result as above


# but [SymPy](https://www.sympy.org) is also available:

# In[45]:


integrate(x^5/(x^3 - 2*x + 1), x, algorithm='sympy')


# as well as [Giac](https://www-fourier.ujf-grenoble.fr/~parisse/giac.html): 

# In[46]:


integrate(x^5/(x^3-2*x+1), x, algorithm='giac')


# One can even use [WolframAlpha](https://www.wolframalpha.com/calculators/integral-calculator/) via some internet connection:

# In[47]:


integrate(x^5/(x^3-2*x+1), x, algorithm='mathematica_free')


# Actually SageMath has some interface to Mathematica:

# In[48]:


cos(x^2)._mathematica_init_()


# as well as to SymPy, Maxima and Giac:

# In[49]:


cos(x^2)._sympy_()


# In[50]:


print(cos(x^2)._maxima_())


# In[51]:


cos(x^2)._giac_()


# Let us check the primitive computation:

# In[52]:


f = integrate(x^5/(x^3-2*x+1), x) 
f


# In[53]:


diff(f, x)


# In[54]:


_.simplify_full()


# Indefinite integrals can also be computed:

# In[55]:


integrate(exp(-x^2), x, -oo, +oo)


# ### Other examples of symbolic computations: 
# 
# - Taylor expansions:

# In[56]:


exp(x).taylor(x, 0, 8)


# - limits:

# In[57]:


lim(sin(x)/x, x=0)


# - series: Riemann's zeta function $\zeta(s)$ for $s=2$ and $s=3$ (Apéry's constant):

# In[58]:


n = var('n')  # declaring n as a symbolic variable
sum(1/n^2, n, 1, +oo, hold=True) == sum(1/n^2, n, 1, +oo)


# In[59]:


sum(1/n^3, n, 1, +oo)


# In[60]:


numerical_approx(_)


# Of course, as many standard functions, Riemann's zeta function is already implemented in SageMath:

# In[61]:


zeta(2)


# - Solving an equation:

# In[62]:


solve(x^2 == x + 1, x)


# - Solving a differential equation:

# In[63]:


y = function('y')
eq = diff(y(x), x) - y(x) == x*y(x)^4
eq


# In[64]:


desolve(eq, y(x))


# ## The power of Python

# Finally, to illustrate the advantage of being built atop of Python, here is a loop for displaying Pascal's triangle with only two instruction lines:

# In[65]:


for n in range(10): 
    print([binomial(n, p) for p in range(n+1)])


# ## Next step:
# 
# ### Plots and functions in SageMath

# In[66]:


plot(sin(x^2), (x, 0, 4), axes_labels=['$x$', '$y$'])


# See the dedicated [notebook](https://nbviewer.org/github/egourgoulhon/SageMathTour/blob/master/Notebooks/plot_tour_2D.ipynb).