Sage Days 114</h2>

Institute of Mathematical Sceinces, Chennai (India)</h3>

Ajit Kumar

Institute of Chemical Technology, Mumbai

What is Sagemath?

  • SageMath is a free and open source sophisticated computer algebra system (CAS).
  • It is built on to top of several free and open source mathemtical packages. 
          For a full list see: https://www.sagemath.org/links-components.html
  • It supports numerical computations, symbolic manipulations, graphing and complex progarmming.
  • Viable open source alternative to: Mathematica, MatLAB, Maple, Magma.
  • It was created by William Stein (https://www.wstein.org/) in 2005 as a platform which uses many free and open source mathematical software.
  • Native language of SageMath is PYTHON.
  • One can do python programming in Sagemath environment and also import Sage as a library in a Python script.
  • SageMath is a very nice mathematical tool which can help in learning, teaching and researching in mathematics.
  • SageMath has very nice interface with LaTeX.
  • It is available for all platforms: Linux, Windows, Mac.

Why to use SageMath?

  • It is FREE and with amazing capability.
  • Being Python based, it is easy to learn and use.
  • SageMath is a very nice Pedagogical tool.
  • Can be used as a very useful research tool.

Let us get started

Dealing with integers

In [1]:
726471264+36816
Out[1]:
726508080
In [2]:
sqrt(5.0)
Out[2]:
2.23606797749979
In [3]:
sqrt(5).n()
Out[3]:
2.23606797749979
In [4]:
show(sqrt(5))
5
In [5]:
a = 615711658618567
In [6]:
a.is_prime()
Out[6]:
False
In [7]:
a.factor()
Out[7]:
7 * 87958808374081
In [8]:
a.previous_prime()
Out[8]:
615711658618493
In [9]:
a.next_prime()
Out[9]:
615711658618597
In [ ]:

Over Finite Fields

In [10]:
K.<a> = GF(5^45);K
Out[10]:
Finite Field in a of size 5^45
In [11]:
a.multiplicative_order()
Out[11]:
229207334116161350281007828251
In [12]:
a.is_square()
Out[12]:
True
In [13]:
b = a.square_root()
b
Out[13]:
a^43 + 4*a^42 + 2*a^40 + 3*a^38 + a^37 + 2*a^36 + 4*a^35 + 3*a^32 + a^31 + 3*a^30 + 3*a^29 + 2*a^28 + 4*a^27 + 4*a^24 + 2*a^23 + 3*a^20 + 3*a^19 + 2*a^18 + 3*a^17 + 4*a^16 + a^15 + 3*a^13 + 2*a^12 + 4*a^11 + 3*a^10 + a^6 + 2*a^5 + 2*a^4 + 3*a^3
In [14]:
b^2
Out[14]:
a
In [15]:
a, b = 282, 384
gcd(a,b)
Out[15]:
6
In [16]:
d, x, y = xgcd(a,b)
In [17]:
print(d,x,y)

6 15 -11
In [18]:
x*a+y*b == d
Out[18]:
True
In [19]:
P.<x> = QQ[]
F = (x^2 + 2)*x^3 
G = (x^2+2)*(x-3)
g, u, v = QQ._xgcd_univariate_polynomial(F,G)
print(f'g={g}, u={u}, v={v}')
u*F + v*G==g # Verification of Bezout's Identity

g=x^2 + 2, u=1/27, v=-1/27*x^2 - 1/9*x - 1/3
Out[19]:
True
In [20]:
P1 = []   # list of pairs (n,n) when n is prime
P2 = []   # list of pairs (n,d) where d|n and n is not prime
for n in xsrange(1,500):
    if n.is_prime():
        P1.append((n,n))
    else:
        P2.extend([(n,d) for d in n.divisors()[1:]])
G1 = point2d(P1,color='red',pointsize=20)
G2 = point2d(P2,color='blue',pointsize=5)
(G1 + G2).show()

Exploring a function of on variable

Let us explore the function

x2+xe(−x2+2)−3x+sin⁡(x2)+1

In [21]:
f(x) = sin(x^2)+x*exp(2-x^2)+x^2-3*x+1
In [ ]:
In [22]:
show(f(x))
x2+xe(−x2+2)−3x+sin⁡(x2)+1
In [23]:
latex(f(x))
Out[23]:
x^{2} + x e^{\left(-x^{2} + 2\right)} - 3 \, x + \sin\left(x^{2}\right) + 1
In [24]:
f.plot(-2,4,color ='red',figsize=5)
Out[24]:
In [25]:
f.diff()
Out[25]:
x |--> -2*x^2*e^(-x^2 + 2) + 2*x*cos(x^2) + 2*x + e^(-x^2 + 2) - 3
In [26]:
f3 = f.diff(3)
In [27]:
f.plot(-1,1,color='red')+f3.plot(-1,1,color='blue')
Out[27]:
In [28]:
f.find_local_maximum(0,2)
Out[28]:
(3.033075870198519, 0.672856583892896)
In [29]:
f.find_root(1,2)
Out[29]:
1.6460113287910914
In [30]:
f.taylor(x,0,8)
Out[30]:
x |--> -1/6*x^7*e^2 - 1/6*x^6 + 1/2*x^5*e^2 - x^3*e^2 + 2*x^2 + x*(e^2 - 3) + 1

Solve Function

In [31]:
#help(solve)
In [32]:
reset()
In [33]:
var('y')
sols = solve([x^3==y,y^2==x], [x,y]); 
sols
Out[33]:
[[x == (0.3090169943749475 + 0.9510565162951535*I), y == (-0.8090169943749475 - 0.5877852522924731*I)], [x == (0.3090169943749475 - 0.9510565162951535*I), y == (-0.8090169943749475 + 0.5877852522924731*I)], [x == (-0.8090169943749475 + 0.5877852522924731*I), y == (0.3090169943749474 + 0.9510565162951536*I)], [x == (-0.8090169943749475 - 0.5877852522924731*I), y == (0.3090169943749475 - 0.9510565162951536*I)], [x == 1, y == 1], [x == 0, y == 0]]

Exploring function of more than one variables

In [34]:
var('x,y')
f(x,y)=(x^2-y^2)*exp(-x^2-y^2)
In [35]:
ct = contour_plot(f(x,y),(x,-2,2),(y,-2,2),contours=12,fill=False,
                  labels=True,cmap='hsv',label_colors='black')
gr = f.gradient()
vf = plot_vector_field(gr,(x,-2,2),(y,-2,2),color='red')
ct+vf
Out[35]:
In [36]:
latex(f(x,y))
Out[36]:
{\left(x^{2} - y^{2}\right)} e^{\left(-x^{2} - y^{2}\right)}
In [37]:
plot3d(f(x,y),(x,-2,2),(y,-2,2))
Out[37]:
In [38]:
f.diff(x)
Out[38]:
(x, y) |--> -2*(x^2 - y^2)*x*e^(-x^2 - y^2) + 2*x*e^(-x^2 - y^2)
In [39]:
f.diff(x,y)
Out[39]:
(x, y) |--> 4*(x^2 - y^2)*x*y*e^(-x^2 - y^2)
In [40]:
f.gradient()
Out[40]:
(x, y) |--> (-2*(x^2 - y^2)*x*e^(-x^2 - y^2) + 2*x*e^(-x^2 - y^2), -2*(x^2 - y^2)*y*e^(-x^2 - y^2) - 2*y*e^(-x^2 - y^2))
In [41]:
f.hessian()(x=1,y=-1)
Out[41]:
[-6*e^(-2)         0]
[        0  6*e^(-2)]
In [ ]:

Working with vectors and Matrices

In [42]:
v = vector(QQ,[-2,3])
u = vector(QQ,[1,-1])
In [43]:
reset()
In [44]:
A = matrix(QQ, [
    [1, 1, 2, -1 ,1, 1],
    [-6, 3, 7, -3, 0, 3],
    [-2, 0, 5, -1, 0, 1],
    [-10, 0, 10, -2, 1, 5],
    [2, 0, -2, 1, 3, -1],
    [-6, 2, 6, -3, 3, 6]])
In [45]:
A.jordan_form()
Out[45]:
[1|0 0 0|0 0]
[-+-----+---]
[0|3 1 0|0 0]
[0|0 3 1|0 0]
[0|0 0 3|0 0]
[-+-----+---]
[0|0 0 0|3 1]
[0|0 0 0|0 3]
In [ ]:
@interact 
def linear_transformation(A=matrix(QQ,[[1,1],[1/2,-1]]),
            u=matrix(QQ,[2,2]),v=matrix(QQ,[-1,2])): 
    u=vector(QQ,[u[0,0],u[0,1]])
    v=vector(QQ,[v[0,0],v[0,1]])
    r =1.1
    e=vector(QQ,[1/10,1/10])
    u1 = A*u
    v1 = A*v
    p1=plot(u,color='red')+plot(v,color='green')
    p1=p1+plot(u+v,color='blue')
    p1=p1+line((u,u+v),linestyle="--",color='black')
    p1=p1+line((v,u+v),linestyle="--",color='black')
    p2=plot(A*u,color='red')+plot(A*v,color='green')
    p2=p2+plot(A*(u+v),color='blue')
    p2=p2+line((A*u,A*(u+v)),linestyle="--",color='black')
    p2=p2+line((A*v,A*(u+v)),linestyle="--",color='black')
    t=text('$u$',r*u)+text('$v$',r*v)+text('$u+v$',r*(u+v))
    t=t+text('$Au$',r*A*u,color='red')+text('$Av$',r*A*v,color='red')+text('$A(u+v)$',r*A*(u+v),color='red')
    show(p1+p2+t,aspect_ratio=1)

Working with Groups

In [46]:
G =SymmetricGroup(6)
In [47]:
a = G.random_element()
In [48]:
a.order()
Out[48]:
6
In [49]:
G = GL(3,5)
S = SL(3,5)
In [50]:
G.order()
S.order()
Out[50]:
372000

Sylow Subgroup in GL(n, GF(5))

In [51]:
n,p=2,5
G = GL(n,GF(p))
print(G.order())
print(factor(G.order()))

480
2^5 * 3 * 5
In [52]:
O2 = [g for g in G if g.multiplicative_order()==2]
O4 = [g for g in G if g.multiplicative_order()==4]
(len(O2),len(O4))
Out[52]:
(31, 152)
In [53]:
for i in range(len(O2)):
    for j in range(len(O4)):
        A,B= O4[i],O4[j]
        H = MatrixGroup([matrix(GF(p),n,flatten(A.list())),
                              matrix(GF(p),n,flatten(B.list()))]); 
        if(H.order()==32):
            break
print(H.order())

32

Lattice of Subgroups

In [54]:
G = CyclicPermutationGroup(30)
subgroups = G.subgroups()
P = Poset((subgroups, lambda h,k: h.is_subgroup(k)))
relabeling = {h:h.order() for (i,h) in enumerate(P)}
Q = P.relabel(relabeling)
LZ30_2=Q.plot(figsize=5,element_size=500)
LZ30_2
Out[54]:

Image processing usng SVD from Python Numpy

In [55]:
from matplotlib.pyplot import imread
import pylab
import numpy as np
img=pylab.imread('Image2.png')
matrix_plot(img)
Out[55]:
In [56]:
R1 = img[:,:,0] # Red channel
G1 = img[:,:,1] # Green Channel
B1 = img[:,:,2] # Blue Channel
In [57]:
R1.shape,G1.shape,B1.shape
Out[57]:
((3456, 4608), (3456, 4608), (3456, 4608))
In [58]:
graphics_array([matrix_plot(R1),matrix_plot(G1),matrix_plot(B1)])
Out[58]:
In [59]:
img.shape
Out[59]:
(3456, 4608, 3)
In [60]:
img_transposed = np.transpose(img, (2, 0, 1))
img_transposed.shape
Out[60]:
(3, 3456, 4608)
In [61]:
U, s, Vt = np.linalg.svd(img_transposed)
In [62]:
Sigma = np.zeros(img_transposed.shape)
for j in range(3):
    np.fill_diagonal(Sigma[j, :, :], s[j, :])
In [65]:
Sigma.shape
Out[65]:
(3, 3456, 4608)
In [66]:
reconstructed = U @ Sigma @ Vt
In [67]:
k= 5
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
pylab.imshow(np.transpose(approx_img, (1, 2, 0)))

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[67]:
<matplotlib.image.AxesImage object at 0x6ffecc0bab10>
In [68]:
k= 10
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
pylab.imshow(np.transpose(approx_img, (1, 2, 0)))

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[68]:
<matplotlib.image.AxesImage object at 0x6ffec1a936d0>
In [69]:
k= 20
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
pylab.imshow(np.transpose(approx_img, (1, 2, 0)))

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[69]:
<matplotlib.image.AxesImage object at 0x6ffec1d4bed0>
In [70]:
k= 50
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
pylab.imshow(np.transpose(approx_img, (1, 2, 0)))

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[70]:
<matplotlib.image.AxesImage object at 0x6ffec9a7d750>
In [71]:
k= 100
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
pylab.imshow(np.transpose(approx_img, (1, 2, 0)))

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[71]:
<matplotlib.image.AxesImage object at 0x6ffec9d6a250>

References

  1. Sagemath documentations: https://doc.sagemath.org/
  2. Computational Mathematics with Sagemath by Paul Zimmerman
  3. Online book on Linear Algebra by Robert Beezer
  4. Series of YouTube video lectures (Twelve) on SageMath by Ajit Kumar
    ajitmathsoft.wordpress.com/sagemath
  5. NPTEL course on Computational Matheatics with Sagemath by Ajit Kumar https://onlinecourses.nptel.ac.in/noc22_ma24/preview
  6. Group Threory: A expedition with SageMath by Ajit Kumar and Vikas Bist.
In [ ]: