MAP-i sage tour

234 days ago by FMC-MAPi

A tour at SAGEMATH

Get SAGE!

 

Sage is available at http://www.sagemath.org

Use its online notebook for trying-out!

http://demo.sagenb.org/

 

Disclaimer:

This notes are an introduction to SageMath, and most of details were ommited. These will be carefully filled in during classes.

Getting started

Assignment & Arithmetic
x=2; x 
        
2
2
x^2==4 
        
True
True
x^3-2!=x^2+2 
        
False
False
2<=3 
        
True
True
3>2 
        
True
True
a=x**4%3; a 
        
1
1
y=pi/2; y 
        
1/2*pi
1/2*pi
n(y) 
        
1.57079632679490
1.57079632679490
MyPi=numerical_approx(pi,prec=150) 
       
MyPi 
        
3.1415926535897932384626433832795028841971694
3.1415926535897932384626433832795028841971694
cos(MyPi) 
        
-1.0000000000000000000000000000000000000000000
-1.0000000000000000000000000000000000000000000
cos(pi) 
        
-1
-1
type(cos(pi)) 
        
<type 'sage.symbolic.expression.Expression'>
<type 'sage.symbolic.expression.Expression'>
type(cos(MyPi)) 
        
<type 'sage.rings.real_mpfr.RealNumber'>
<type 'sage.rings.real_mpfr.RealNumber'>

Changing types on-the-fly...
z=3 
       
type(z) 
        
<type 'sage.rings.integer.Integer'>
<type 'sage.rings.integer.Integer'>
z=z+2/3; z 
        
11/3
11/3
type(z) 
        
<type 'sage.rings.rational.Rational'>
<type 'sage.rings.rational.Rational'>

Changing basis...
num1=10; num2=7 [num1.str(2),num2.str(2)] 
        
['1010', '111']
['1010', '111']
num1.str(2)+num2.str(2) 
        
'1010111'
'1010111'

HEEEEEEEEEELP!!!!!!!!!!

Just add "?" at the end of the command.

Use autocompletion pressing the TAB key.
IntegerRing? 
       

File: /home/pedro/sage-4.7/devel/sage/sage/rings/integer_ring.pyx

Type: <type ‘builtin_function_or_method’>

Definition: IntegerRing()

Docstring:

Return the integer ring

EXAMPLE:

sage: IntegerRing()
Integer Ring
sage: ZZ==IntegerRing()
True

File: /home/pedro/sage-4.7/devel/sage/sage/rings/integer_ring.pyx

Type: <type ‘builtin_function_or_method’>

Definition: IntegerRing()

Docstring:

Return the integer ring

EXAMPLE:

sage: IntegerRing()
Integer Ring
sage: ZZ==IntegerRing()
True
factor? 
       

File: /home/pedro/sage-4.7/local/lib/python2.6/site-packages/sage/rings/arith.py

Type: <type ‘function’>

Definition: factor(n, proof=None, int_=False, algorithm=’pari’, verbose=0, **kwds)

Docstring:

Returns the factorization of n. The result depends on the type of n.

If n is an integer, returns the factorization as an object of type Factorization.

If n is not an integer, n.factor(proof=proof, **kwds) gets called. See n.factor?? for more documentation in this case.

Warning

This means that applying factor to an integer result of a symbolic computation will not factor the integer, because it is considered as an element of a larger symbolic ring.

EXAMPLE:

sage: f(n)=n^2
sage: is_prime(f(3))
False
sage: factor(f(3))
9

INPUT:

  • n - an nonzero integer
  • proof - bool or None (default: None)
  • int_ - bool (default: False) whether to return answers as Python ints
  • algorithm - string
    • 'pari' - (default) use the PARI c library
    • 'kash' - use KASH computer algebra system (requires the optional kash package be installed)
    • 'magma' - use Magma (requires magma be installed)
  • verbose - integer (default: 0); PARI’s debug variable is set to this; e.g., set to 4 or 8 to see lots of output during factorization.

OUTPUT:

  • factorization of n

The qsieve and ecm commands give access to highly optimized implementations of algorithms for doing certain integer factorization problems. These implementations are not used by the generic factor command, which currently just calls PARI (note that PARI also implements sieve and ecm algorithms, but they aren’t as optimized). Thus you might consider using them instead for certain numbers.

The factorization returned is an element of the class Factorization; see Factorization?? for more details, and examples below for usage. A Factorization contains both the unit factor (+1 or -1) and a sorted list of (prime, exponent) pairs.

The factorization displays in pretty-print format but it is easy to obtain access to the (prime,exponent) pairs and the unit, to recover the number from its factorization, and even to multiply two factorizations. See examples below.

EXAMPLES:

sage: factor(500)
2^2 * 5^3
sage: factor(-20)
-1 * 2^2 * 5
sage: f=factor(-20)
sage: list(f)
[(2, 2), (5, 1)]
sage: f.unit()
-1
sage: f.value()
-20
sage: factor( -next_prime(10^2) * next_prime(10^7) )
-1 * 101 * 10000019

sage: factor(-500, algorithm='kash')      # optional - kash
-1 * 2^2 * 5^3

sage: factor(-500, algorithm='magma')     # optional - magma
-1 * 2^2 * 5^3

sage: factor(0)
Traceback (click to the left of this block for traceback)
...

File: /home/pedro/sage-4.7/local/lib/python2.6/site-packages/sage/rings/arith.py

Type: <type ‘function’>

Definition: factor(n, proof=None, int_=False, algorithm=’pari’, verbose=0, **kwds)

Docstring:

Returns the factorization of n. The result depends on the type of n.

If n is an integer, returns the factorization as an object of type Factorization.

If n is not an integer, n.factor(proof=proof, **kwds) gets called. See n.factor?? for more documentation in this case.

Warning

This means that applying factor to an integer result of a symbolic computation will not factor the integer, because it is considered as an element of a larger symbolic ring.

EXAMPLE:

sage: f(n)=n^2
sage: is_prime(f(3))
False
sage: factor(f(3))
9

INPUT:

  • n - an nonzero integer
  • proof - bool or None (default: None)
  • int_ - bool (default: False) whether to return answers as Python ints
  • algorithm - string
    • 'pari' - (default) use the PARI c library
    • 'kash' - use KASH computer algebra system (requires the optional kash package be installed)
    • 'magma' - use Magma (requires magma be installed)
  • verbose - integer (default: 0); PARI’s debug variable is set to this; e.g., set to 4 or 8 to see lots of output during factorization.

OUTPUT:

  • factorization of n

The qsieve and ecm commands give access to highly optimized implementations of algorithms for doing certain integer factorization problems. These implementations are not used by the generic factor command, which currently just calls PARI (note that PARI also implements sieve and ecm algorithms, but they aren’t as optimized). Thus you might consider using them instead for certain numbers.

The factorization returned is an element of the class Factorization; see Factorization?? for more details, and examples below for usage. A Factorization contains both the unit factor (+1 or -1) and a sorted list of (prime, exponent) pairs.

The factorization displays in pretty-print format but it is easy to obtain access to the (prime,exponent) pairs and the unit, to recover the number from its factorization, and even to multiply two factorizations. See examples below.

EXAMPLES:

sage: factor(500)
2^2 * 5^3
sage: factor(-20)
-1 * 2^2 * 5
sage: f=factor(-20)
sage: list(f)
[(2, 2), (5, 1)]
sage: f.unit()
-1
sage: f.value()
-20
sage: factor( -next_prime(10^2) * next_prime(10^7) )
-1 * 101 * 10000019

sage: factor(-500, algorithm='kash')      # optional - kash
-1 * 2^2 * 5^3

sage: factor(-500, algorithm='magma')     # optional - magma
-1 * 2^2 * 5^3

sage: factor(0)
Traceback (most recent call last):
...
ArithmeticError: Prime factorization of 0 not defined.
sage: factor(1)
1
sage: factor(-1)
-1
sage: factor(2^(2^7)+1)
59649589127497217 * 5704689200685129054721

Sage calls PARI’s factor, which has proof False by default. Sage has a global proof flag, set to True by default (see sage.structure.proof.proof, or proof.[tab]). To override the default, call this function with proof=False.

sage: factor(3^89-1, proof=False)
2 * 179 * 1611479891519807 * 5042939439565996049162197

sage: factor(2^197 + 1)  # long time
3 * 197002597249 * 1348959352853811313 * 251951573867253012259144010843

Any object which has a factor method can be factored like this:

sage: K.<i> = QuadraticField(-1)
sage: factor(122+454*i)
(-1) * (-2*i - 3) * (-i + 4) * (i + 1)^3 * (-i + 2)^3

To access the data in a factorization:

sage: f = factor(420); f
2^2 * 3 * 5 * 7
sage: [x for x in f]
[(2, 2), (3, 1), (5, 1), (7, 1)]
sage: [p for p,e in f]
[2, 3, 5, 7]
sage: [e for p,e in f]
[2, 1, 1, 1]
sage: [p^e for p,e in f]
[4, 3, 5, 7]

Create functions

Define your own functions.
def divtrivial(k): primo=true for p in prime_range(1,floor(sqrt(k))+1): if k%p==0: primo=false print p, " divide ", k return primo 
       
 
       
def divtrivial2(k): primo=true p=1 while primo and p<=floor(sqrt(k)): p=next_prime(p) print "testando com p=",p if k%p==0: primo=false print p," divide ",k return primo 
       
 
       

Solving equations
x, y = var('x, y') solve([2*x+y==6, x-y==4], x, y) 
        
[[x == (10/3), y == (-2/3)]]
[[x == (10/3), y == (-2/3)]]
x, b, c = var('x b c') solve([x^2 + b*x + c == 0],x) 
        
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]
[x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]

Drawing and painting
c = circle((0,0), 1, rgbcolor=(0.5,1,0.5)) c.show() 
       

Scaling matters!
c.show(aspect_ratio=1) 
       

The command show(c, aspect_ratio=1) accomplishes the same thing, or you can save the picture using c.save('filename.png',aspect_ratio=1).
plot(prime_pi, 2, 10000) 
       
plot((x/log(x)),2,10000) 
       
A=plot(prime_pi, 1, 100000, rgbcolor='blue') B=plot((x/log(x)),2,100000, rgbcolor='green') show(A+B) 
       
P=plot(Li, 2, 10000, rgbcolor='red') Q=plot(prime_pi, 2, 10000, rgbcolor='blue') R=plot(sqrt(x)*(log(x)),2,10000, rgbcolor='green') show(P+Q+R,xmin=0, figsize=[10,5]) 
       
x = var('x') parametric_plot((cos(x),sin(x)^3),(x,0,2*pi),rgbcolor=hue(0.6)) 
       
x = var('x') p1 = parametric_plot((cos(x),sin(x)),(x,0,2*pi),rgbcolor=hue(0.2)) p2 = parametric_plot((cos(x),sin(x)^2),(x,0,2*pi),rgbcolor=hue(0.4)) p3 = parametric_plot((cos(x),sin(x)^3),(x,0,2*pi),rgbcolor=hue(0.6)) show(p1+p2+p3, axes=false) 
       

3D!
x, y = var('x,y') plot3d(x^2 + y^2, (x,-2,2), (y,-2,2)) 
       
u, v = var('u,v') fx = (3+sin(v)+cos(u))*cos(2*v) fy = (3+sin(v)+cos(u))*sin(2*v) fz = sin(u)+2*cos(v) parametric_plot3d([fx, fy, fz], (u, 0, 2*pi), (v, 0, 2*pi), frame=False, color="red") 
       
x, y, z = var('x,y,z') f(x, y, z) = 4*x^2 * (x^2 + y^2 + z^2 + z) + y^2 * (y^2 + z^2 - 1) implicit_plot3d(f, (x, -0.5, 0.5), (y, -1, 1), (z, -1, 1)) 
       

Interact
import random def ftree(rows, v, i, F): if len(v) > 0: # add a row to g at the ith level. rows.append(v) w = [] for i in range(len(v)): k, _, _ = v[i] if k is None or is_prime(k): w.append((None,None,None)) else: d = random.choice(divisors(k)[1:-1]) w.append((d,k,i)) e = k//d if e == 1: w.append((None,None)) else: w.append((e,k,i)) if len(w) > len(v): ftree(rows, w, i+1, F) def draw_ftree(rows,font): g = Graphics() for i in range(len(rows)): cur = rows[i] for j in range(len(cur)): e, f, k = cur[j] if not e is None: if is_prime(e): c = (1,0,0) else: c = (0,0,.4) g += text(str(e), (j*2-len(cur),-i), fontsize=font, rgbcolor=c) if not k is None and not f is None: g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)], alpha=0.5) return g @interact def factor_tree(n=100, font=(10, (8..20)), redraw=['Redraw']): n = Integer(n) rows = [] v = [(n,None,0)] ftree(rows, v, 0, factor(n)) show(draw_ftree(rows, font), axes=False) 
       

Click to the left again to hide and once more to show the dynamic interactive window
@interact def diffie_hellman(bits=slider(8, 513, 4, 8, 'Numero de bits', False), button=selector(["Novo exemplo"],label='',buttons=True)): maxp = 2 ^ bits p = random_prime(maxp) k = GF(p) if bits > 100: g = k(2) else: g = k.multiplicative_generator() a = ZZ.random_element(10, maxp) b = ZZ.random_element(10, maxp) print """ <html> <style> .gamodp, .gbmodp { color:#000; padding:5px } .gamodp { background:#846FD8 } .gbmodp { background:#FFFC73 } .dhsame { color:#000; font-weight:bold } </style> <h2 style="color:#000;font-family:Arial, Helvetica, sans-serif">Troca de chaves de Diffie-Hellman %s-Bits</h2> <ol style="color:#000;font-family:Arial, Helvetica, sans-serif"> <li>Alice e Bob escolhem um numero primo p = %s e uma base g = %s.</li> <li>Alice escolhe, secretamente, um inteiro a = %s, e envia a Bob (<span class="gamodp">g<sup>a</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gamodp">%s</span>.</li> <li>Bob escolhe, secretamente, um inteiro b=%s, e envia a Alice (<span class="gbmodp">g<sup>b</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gbmodp">%s</span>.</li> <li>Alice calcula (<span class="gbmodp">g<sup>b</sup> mod p</span>)<sup>a</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> <li>Bob calcula (<span class="gamodp">g<sup>a</sup> mod p</span>)<sup>b</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> </ol></html> """ % (bits, p, g, a, g, a, p, (g^a), b, g, b, p, (g^b), (g^b), a, p, (g^ b)^a, g^a, b, p, (g^a)^b) 
       

Click to the left again to hide and once more to show the dynamic interactive window
@interact def RSA(bits=slider(8, 513, 4, 8, 'Numero de bits', False), button=selector(["Novo exemplo"],label='',buttons=True)): maxp = 2 ^ ((bits/2).floor()) maxq = 2 ^ ((bits/2).floor()) p = random_prime(maxp) q = random_prime(maxq) while q==p: q = random_prime(maxq) n=p*q m=(p-1)*(q-1) mens=randint(2,n) a=randint(2,m) while gcd(a,m)!= 1: a=randint(2,m) b=Mod(1/a,m) cifrado=Mod(mens,n)^a decifrado=Mod(cifrado,n)^b print """ <html> <h2 style="color:#000;font-family:Arial, Helvetica, sans-serif">RSA com %s Bits</h2> <ol style="color:#000;font-family:Arial, Helvetica, sans-serif"> <li> Alice escolhe dois primos p = %s e q = %s, e calcula <br/> n = p*q = %s. <br/> Calcula ainda <br/> m = (p-1)*(q-1) = %s e escolhe a = %s primo relativo com <br/> m = %s.<br/> Calcula b = a<sup>-1</sup> mod m = %s.<br/> Publica n = %s e a = %s (chave publica de Alice)</li> <li> Bob pretende enviar a mensagem mens = %s a Alice. Com a chave publica de Alice, calcula<br/> cifrado = mens<sup>a</sup> mod n = %s.</li> <li> Alice recebe o texto cifrado = %s de Bob. Usando o inverso b = %s = a<sup>-1</sup> mod m, calcula <br/> decifrado = cifrado<sup>b</sup> mod n = %s. </ol></html> """ % (bits, p, q, n, m, a, m, b, n, a, mens, cifrado, cifrado, b, decifrado) 
       

Click to the left again to hide and once more to show the dynamic interactive window
@interact def _(N=(100,(2..2000))): html("<font color='red'>$\pi(x)$</font> e <font color='blue'>$x/(\log(x)-1)$</font> para $x < %s$"%N) show(plot(prime_pi, 0, N, rgbcolor='red') + plot(x/(log(x)-1), 5, N, rgbcolor='blue')) # obtido de http://wiki.sagemath.org/interact/number_theory#FactoringanInteger 
       

Click to the left again to hide and once more to show the dynamic interactive window

Rings, fields and vector spaces
Z7=IntegerModRing(7) 
       
a=Z7(3) 
       
1/a 
        
5
5
Z7.is_field() 
        
True
True
Z8=IntegerModRing(8) 
       
Z8.is_field() 
        
False
False
F8=GF(2^3,'a'); F8 
        
Finite Field in a of size 2^3
Finite Field in a of size 2^3
[i for i in F8] 
        
[0, a, a^2, a + 1, a^2 + a, a^2 + a + 1, a^2 + 1, 1]
[0, a, a^2, a + 1, a^2 + a, a^2 + a + 1, a^2 + 1, 1]
F8.characteristic() 
        
2
2

Constructing the Galois field.
p=2; n=3 
       
R.<x>=PolynomialRing(GF(p)) 
       
pol=x^n+x+1 
       
type(pol) 
        
<type 'sage.rings.polynomial.polynomial_gf2x.Polynomial_GF2X'>
<type 'sage.rings.polynomial.polynomial_gf2x.Polynomial_GF2X'>
pol in R 
        
True
True
pol.is_irreducible() 
        
True
True
pol.factor() 
        
x^3 + x + 1
x^3 + x + 1
S.<x>=R.quotient(pol) S 
        
Univariate Quotient Polynomial Ring in x over Finite Field of size 2
with modulus x^3 + x + 1
Univariate Quotient Polynomial Ring in x over Finite Field of size 2 with modulus x^3 + x + 1
f=x+1 
       
type(f) 
        
<class
'sage.rings.polynomial.polynomial_quotient_ring_element.PolynomialQuotie\
ntRingElement'>
<class 'sage.rings.polynomial.polynomial_quotient_ring_element.PolynomialQuotientRingElement'>
f in S 
        
True
True
type(pol) 
        
<type 'sage.rings.polynomial.polynomial_gf2x.Polynomial_GF2X'>
<type 'sage.rings.polynomial.polynomial_gf2x.Polynomial_GF2X'>
1/f 
        
x^2 + x
x^2 + x
S.order() 
        
8
8
[(i,f^i) for i in [1..7]] 
        
[(1, x + 1), (2, x^2 + 1), (3, x^2), (4, x^2 + x + 1), (5, x), (6, x^2 +
x), (7, 1)]
[(1, x + 1), (2, x^2 + 1), (3, x^2), (4, x^2 + x + 1), (5, x), (6, x^2 + x), (7, 1)]
g=S(x) 
       
1/g 
        
x^2 + 1
x^2 + 1
f*g 
        
x^2 + x
x^2 + x
f^3+g^(-2) 
        
x + 1
x + 1
[(i,g^i) for i in [1..7]] 
        
[(1, x), (2, x^2), (3, x + 1), (4, x^2 + x), (5, x^2 + x + 1), (6, x^2 +
1), (7, 1)]
[(1, x), (2, x^2), (3, x + 1), (4, x^2 + x), (5, x^2 + x + 1), (6, x^2 + 1), (7, 1)]

We have considered the euclidean domain \mathbb{Z}_2 [x] and its irreducible element pol(x)=x^3+x+1, which in turn is equivalent to the ideal (pol(x)) being maximal.

The quocient \mathbb{F}=\mathbb{Z}_2[x] /(pol(x)) is a field and \mathbb{F}^* is a cyclic group with order 2^3-1=7.
K.<x>=FiniteField(2^7,name='x', modulus='conway') 
       
K.modulus() 
        
x^7 + x + 1
x^7 + x + 1
def chaves(p,n): if is_prime(p): R.<x>=FiniteField(p^n, name='x', modulus='conway') exp=IntegerModRing(p^n-1).random_element() return [charpoly(x),R,exp] 
       
minhachave=chaves(2,10) minhachave 
        
[x^10 + x^6 + x^5 + x^3 + x^2 + x + 1, Finite Field in x of size 2^10,
827]
[x^10 + x^6 + x^5 + x^3 + x^2 + x + 1, Finite Field in x of size 2^10, 827]
corpo.<x>=minhachave[1] a=minhachave[2] 
       
mess=corpo(x^3+x+1) mess 
        
x^3 + x + 1
x^3 + x + 1

Linear algebra
M=Matrix(ZZ,[[1, 3, 2, 6],[3, 5, 2, 1],[6,0,1,-1]]) M 
        
[ 1  3  2  6]
[ 3  5  2  1]
[ 6  0  1 -1]
[ 1  3  2  6]
[ 3  5  2  1]
[ 6  0  1 -1]
M.echelon_form() 
        
[ 1  1  7 37]
[ 0  2  9 48]
[ 0  0 14 79]
[ 1  1  7 37]
[ 0  2  9 48]
[ 0  0 14 79]
M.row_space() 
        
Free module of degree 4 and rank 3 over Integer Ring
Echelon basis matrix:
[ 1  1  7 37]
[ 0  2  9 48]
[ 0  0 14 79]
Free module of degree 4 and rank 3 over Integer Ring
Echelon basis matrix:
[ 1  1  7 37]
[ 0  2  9 48]
[ 0  0 14 79]
Z7=IntegerModRing(11) 
       
M=matrix(Z7,[[1, 3, 2, 6],[3, 5, 2, 1],[6,0,1,-1]]) M 
        
[ 1  3  2  6]
[ 3  5  2  1]
[ 6  0  1 10]
[ 1  3  2  6]
[ 3  5  2  1]
[ 6  0  1 10]
transpose(M).nullity() 
        
1
1
M.rref() 
        
[ 1  0  0  4]
[ 0  1  0 10]
[ 0  0  1  8]
[ 1  0  0  4]
[ 0  1  0 10]
[ 0  0  1  8]
K.<x>=FiniteField(2^7,name='x', modulus='conway') 
       
K 
        
Finite Field in x of size 2^7
Finite Field in x of size 2^7
K.modulus() 
        
x^7 + x + 1
x^7 + x + 1
MK=MatrixSpace(K,3) MK 
        
Full MatrixSpace of 3 by 3 dense matrices over Finite Field in x of size
2^7
Full MatrixSpace of 3 by 3 dense matrices over Finite Field in x of size 2^7
A=MK.random_element() A 
        
[  x^5 + x^4 + x^3 + x^2 + x + 1                     x^3 + x + 1 x^6 +
x^5 + x^4 + x^3 + x^2 + x]
[                        x^4 + x             x^5 + x^4 + x^3 + 1        
x^6 + x^5 + 1]
[        x^6 + x^3 + x^2 + x + 1   x^6 + x^5 + x^3 + x^2 + x + 1   x^6 +
x^5 + x^4 + x^2 + x + 1]
[  x^5 + x^4 + x^3 + x^2 + x + 1                     x^3 + x + 1 x^6 + x^5 + x^4 + x^3 + x^2 + x]
[                        x^4 + x             x^5 + x^4 + x^3 + 1                   x^6 + x^5 + 1]
[        x^6 + x^3 + x^2 + x + 1   x^6 + x^5 + x^3 + x^2 + x + 1   x^6 + x^5 + x^4 + x^2 + x + 1]
A.det() 
        
x^6 + x^4 + 1
x^6 + x^4 + 1
A[2]=A[0]+A[1] A 
        
[  x^5 + x^4 + x^3 + x^2 + x + 1                     x^3 + x + 1 x^6 +
x^5 + x^4 + x^3 + x^2 + x]
[                        x^4 + x             x^5 + x^4 + x^3 + 1        
x^6 + x^5 + 1]
[            x^5 + x^3 + x^2 + 1                   x^5 + x^4 + x        
x^4 + x^3 + x^2 + x + 1]
[  x^5 + x^4 + x^3 + x^2 + x + 1                     x^3 + x + 1 x^6 + x^5 + x^4 + x^3 + x^2 + x]
[                        x^4 + x             x^5 + x^4 + x^3 + 1                   x^6 + x^5 + 1]
[            x^5 + x^3 + x^2 + 1                   x^5 + x^4 + x         x^4 + x^3 + x^2 + x + 1]
A.det() 
        
0
0
pretty_print(A.echelon_form()) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr}
1 & 0 & x^{6} + x^{4} + x^{3} + x^{2} \\
0 & 1 & x^{5} + x^{4} + x^{2} + x + 1 \\
0 & 0 & 0
\end{array}\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr}
1 & 0 & x^{6} + x^{4} + x^{3} + x^{2} \\
0 & 1 & x^{5} + x^{4} + x^{2} + x + 1 \\
0 & 0 & 0
\end{array}\right)
latex(A.echelon_form()) 
        
\left(\begin{array}{rrr}
1 & 0 & x^{6} + x^{4} + x^{3} + x^{2} \\
0 & 1 & x^{5} + x^{4} + x^{2} + x + 1 \\
0 & 0 & 0
\end{array}\right)
\left(\begin{array}{rrr}
1 & 0 & x^{6} + x^{4} + x^{3} + x^{2} \\
0 & 1 & x^{5} + x^{4} + x^{2} + x + 1 \\
0 & 0 & 0
\end{array}\right)
 
       

Numbers!

Euclid's algorithm allows to compute the gcd of the integers a,b, n.s.z., and also to solve the linear equation ax+by=(a,b), with x,y\in \mathbb{Z}.
a=356; b=281; g,x,y=xgcd(a,b); g,x,y 
        
(1, 15, -19)
(1, 15, -19)
a*x+b*y 
        
1
1

Set P as the set of prime numbers..
P=Primes(); P 
        
Set of all prime numbers: 2, 3, 5, 7, ...
Set of all prime numbers: 2, 3, 5, 7, ...
[n for n in range(2,40) if not is_prime(n)] 
        
[4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30,
32, 33, 34, 35, 36, 38, 39]
[4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 38, 39]
[n for n in range(2,40) if n in P] 
        
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]

Create your own function that show the non primes \le n.
def listaP(m): lista=[n for n in range(2,m) if is_prime(n)] return lista 
       

Twin primes...
gemeos=[[n,n+2] for n in range(2,200) if n in P and n+2 in P]; gemeos 
        
[[3, 5], [5, 7], [11, 13], [17, 19], [29, 31], [41, 43], [59, 61], [71,
73], [101, 103], [107, 109], [137, 139], [149, 151], [179, 181], [191,
193], [197, 199]]
[[3, 5], [5, 7], [11, 13], [17, 19], [29, 31], [41, 43], [59, 61], [71, 73], [101, 103], [107, 109], [137, 139], [149, 151], [179, 181], [191, 193], [197, 199]]
def gemeos(m): P=Primes() lista=[[n,n+2] for n in range(2,2000) if n in P and n+2 in P] return lista, len(lista) 
       
gemeos(100) 
        
([[3, 5], [5, 7], [11, 13], [17, 19], [29, 31], [41, 43], [59, 61], [71,
73], [101, 103], [107, 109], [137, 139], [149, 151], [179, 181], [191,
193], [197, 199], [227, 229], [239, 241], [269, 271], [281, 283], [311,
313], [347, 349], [419, 421], [431, 433], [461, 463], [521, 523], [569,
571], [599, 601], [617, 619], [641, 643], [659, 661], [809, 811], [821,
823], [827, 829], [857, 859], [881, 883], [1019, 1021], [1031, 1033],
[1049, 1051], [1061, 1063], [1091, 1093], [1151, 1153], [1229, 1231],
[1277, 1279], [1289, 1291], [1301, 1303], [1319, 1321], [1427, 1429],
[1451, 1453], [1481, 1483], [1487, 1489], [1607, 1609], [1619, 1621],
[1667, 1669], [1697, 1699], [1721, 1723], [1787, 1789], [1871, 1873],
[1877, 1879], [1931, 1933], [1949, 1951], [1997, 1999]], 61)
([[3, 5], [5, 7], [11, 13], [17, 19], [29, 31], [41, 43], [59, 61], [71, 73], [101, 103], [107, 109], [137, 139], [149, 151], [179, 181], [191, 193], [197, 199], [227, 229], [239, 241], [269, 271], [281, 283], [311, 313], [347, 349], [419, 421], [431, 433], [461, 463], [521, 523], [569, 571], [599, 601], [617, 619], [641, 643], [659, 661], [809, 811], [821, 823], [827, 829], [857, 859], [881, 883], [1019, 1021], [1031, 1033], [1049, 1051], [1061, 1063], [1091, 1093], [1151, 1153], [1229, 1231], [1277, 1279], [1289, 1291], [1301, 1303], [1319, 1321], [1427, 1429], [1451, 1453], [1481, 1483], [1487, 1489], [1607, 1609], [1619, 1621], [1667, 1669], [1697, 1699], [1721, 1723], [1787, 1789], [1871, 1873], [1877, 1879], [1931, 1933], [1949, 1951], [1997, 1999]], 61)

Interact...
@interact def gemeos(m=slider(5, 10000, 2, 3, 'm', False), button=selector(["1000","10000"],label='aa',buttons=False)): P=Primes() print "m =",m #, button print([[n,n+2] for n in range(2,m) if n in P and n+2 in P]) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Fermat's little theorem: n prime implies a^{n-1}\equiv 1 \text{ mod } n, if (a,n)=1., This gives an elementary way for checking if a number is not prime.
n=12^31+1; n 
        
2848515765597237675947403497177089
2848515765597237675947403497177089
2^(n-1)%n 
        
Traceback (click to the left of this block for traceback)
...
RuntimeError: exponent must be at most 9223372036854775807
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_4.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Ml4obi0xKSVu"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpnxvD2p/___code___.py", line 3, in <module>
    exec compile(u'_sage_const_2 **(n-_sage_const_1 )%n
  File "", line 1, in <module>
    
  File "integer.pyx", line 1839, in sage.rings.integer.Integer.__pow__ (sage/rings/integer.c:13120)
RuntimeError: exponent must be at most 9223372036854775807

This gives an overflow error. We may avoid this by using a fast modular exponential algorithm.

Or by taking, in sagemath, the right algebraic struture.
R=IntegerModRing(n); R 
        
Ring of integers modulo 2848515765597237675947403497177089
Ring of integers modulo 2848515765597237675947403497177089
R(2)^(n-1) 
        
825326789816013506777461377270418
825326789816013506777461377270418
2.powermod(n-1,n) 
        
825326789816013506777461377270418
825326789816013506777461377270418
factor(n) 
        
13 * 1263499 * 173420475484059478781201647
13 * 1263499 * 173420475484059478781201647
R(2)^(n-1)==R(1) 
        
False
False

We may also compute the multicative inverse of an element, in case it exists.
1/R(322) 
        
2114271018564409330905060359705976
2114271018564409330905060359705976

Revisiting Erastosthenes
lista=[n for n in range(0,1000)]; lista 
        
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55,
56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73,
74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107,
108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121,
122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135,
136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163,
164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177,
178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191,
192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205,
206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219,
220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,
234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247,
248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261,
262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275,
276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289,
290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303,
304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317,
318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331,
332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345,
346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359,
360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373,
374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387,
388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401,
402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415,
416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429,
430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443,
444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457,
458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471,
472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485,
486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499,
500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513,
514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527,
528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541,
542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555,
556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569,
570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583,
584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597,
598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611,
612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625,
626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639,
640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653,
654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667,
668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681,
682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695,
696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709,
710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737,
738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751,
752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765,
766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779,
780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793,
794, 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807,
808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821,
822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835,
836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849,
850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863,
864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877,
878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891,
892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905,
906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919,
920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933,
934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947,
948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961,
962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975,
976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989,
990, 991, 992, 993, 994, 995, 996, 997, 998, 999]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999]
lista[0] 
        
0
0
floor(sqrt(1000)) 
        
31
31
for p in prime_range(sqrt(1001)): for k in range(2,floor(1001/p)): lista[k*p]=0 
       
lista 
        
[0, 1, 2, 3, 0, 5, 0, 7, 0, 0, 0, 11, 0, 13, 0, 0, 0, 17, 0, 19, 0, 0,
0, 23, 0, 0, 0, 0, 0, 29, 0, 31, 0, 0, 0, 0, 0, 37, 0, 0, 0, 41, 0, 43,
0, 0, 0, 47, 0, 0, 0, 0, 0, 53, 0, 0, 0, 0, 0, 59, 0, 61, 0, 0, 0, 0, 0,
67, 0, 0, 0, 71, 0, 73, 0, 0, 0, 0, 0, 79, 0, 0, 0, 83, 0, 0, 0, 0, 0,
89, 0, 0, 0, 0, 0, 0, 0, 97, 0, 0, 0, 101, 0, 103, 0, 0, 0, 107, 0, 109,
0, 0, 0, 113, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 127, 0, 0, 0, 131,
0, 0, 0, 0, 0, 137, 0, 139, 0, 0, 0, 0, 0, 0, 0, 0, 0, 149, 0, 151, 0,
0, 0, 0, 0, 157, 0, 0, 0, 0, 0, 163, 0, 0, 0, 167, 0, 0, 0, 0, 0, 173,
0, 0, 0, 0, 0, 179, 0, 181, 0, 0, 0, 0, 0, 0, 0, 0, 0, 191, 0, 193, 0,
0, 0, 197, 0, 199, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 211, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 223, 0, 0, 0, 227, 0, 229, 0, 0, 0, 233, 0, 0, 0, 0,
0, 239, 0, 241, 0, 0, 0, 0, 0, 0, 0, 0, 0, 251, 0, 0, 0, 0, 0, 257, 0,
0, 0, 0, 0, 263, 0, 0, 0, 0, 0, 269, 0, 271, 0, 0, 0, 0, 0, 277, 0, 0,
0, 281, 0, 283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 293, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 307, 0, 0, 0, 311, 0, 313, 0, 0, 0, 317, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 331, 0, 0, 0, 0, 0, 337, 0, 0, 0, 0, 0, 0, 0, 0,
0, 347, 0, 349, 0, 0, 0, 353, 0, 0, 0, 0, 0, 359, 0, 0, 0, 0, 0, 0, 0,
367, 0, 0, 0, 0, 0, 373, 0, 0, 0, 0, 0, 379, 0, 0, 0, 383, 0, 0, 0, 0,
0, 389, 0, 0, 0, 0, 0, 0, 0, 397, 0, 0, 0, 401, 0, 0, 0, 0, 0, 0, 0,
409, 0, 0, 0, 0, 0, 0, 0, 0, 0, 419, 0, 421, 0, 0, 0, 0, 0, 0, 0, 0, 0,
431, 0, 433, 0, 0, 0, 0, 0, 439, 0, 0, 0, 443, 0, 0, 0, 0, 0, 449, 0, 0,
0, 0, 0, 0, 0, 457, 0, 0, 0, 461, 0, 463, 0, 0, 0, 467, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 479, 0, 0, 0, 0, 0, 0, 0, 487, 0, 0, 0, 491, 0, 0, 0,
0, 0, 0, 0, 499, 0, 0, 0, 503, 0, 0, 0, 0, 0, 509, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 521, 0, 523, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 541, 0, 0, 0, 0, 0, 547, 0, 0, 0, 0, 0, 0, 0, 0, 0, 557, 0, 0, 0, 0,
0, 563, 0, 0, 0, 0, 0, 569, 0, 571, 0, 0, 0, 0, 0, 577, 0, 0, 0, 0, 0,
0, 0, 0, 0, 587, 0, 0, 0, 0, 0, 593, 0, 0, 0, 0, 0, 599, 0, 601, 0, 0,
0, 0, 0, 607, 0, 0, 0, 0, 0, 613, 0, 0, 0, 617, 0, 619, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 631, 0, 0, 0, 0, 0, 0, 0, 0, 0, 641, 0, 643, 0, 0, 0,
647, 0, 0, 0, 0, 0, 653, 0, 0, 0, 0, 0, 659, 0, 661, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 673, 0, 0, 0, 677, 0, 0, 0, 0, 0, 683, 0, 0, 0, 0, 0, 0,
0, 691, 0, 0, 0, 0, 0, 0, 0, 0, 0, 701, 0, 0, 0, 0, 0, 0, 0, 709, 0, 0,
0, 0, 0, 0, 0, 0, 0, 719, 0, 0, 0, 0, 0, 0, 0, 727, 0, 0, 0, 0, 0, 733,
0, 0, 0, 0, 0, 739, 0, 0, 0, 743, 0, 0, 0, 0, 0, 0, 0, 751, 0, 0, 0, 0,
0, 757, 0, 0, 0, 761, 0, 0, 0, 0, 0, 0, 0, 769, 0, 0, 0, 773, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 787, 0, 0, 0, 0, 0, 0, 0, 0, 0, 797, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 809, 0, 811, 0, 0, 0, 0, 0, 0, 0, 0, 0, 821,
0, 823, 0, 0, 0, 827, 0, 829, 0, 0, 0, 0, 0, 0, 0, 0, 0, 839, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 853, 0, 0, 0, 857, 0, 859, 0, 0, 0, 863,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 877, 0, 0, 0, 881, 0, 883, 0, 0,
0, 887, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 907, 0,
0, 0, 911, 0, 0, 0, 0, 0, 0, 0, 919, 0, 0, 0, 0, 0, 0, 0, 0, 0, 929, 0,
0, 0, 0, 0, 0, 0, 937, 0, 0, 0, 941, 0, 0, 0, 0, 0, 947, 0, 0, 0, 0, 0,
953, 0, 0, 0, 0, 0, 0, 0, 961, 0, 0, 0, 0, 0, 967, 0, 0, 0, 971, 0, 0,
0, 0, 0, 977, 0, 0, 0, 0, 0, 983, 0, 0, 0, 0, 0, 989, 0, 991, 0, 0, 0,
0, 0, 997, 0, 999]
[0, 1, 2, 3, 0, 5, 0, 7, 0, 0, 0, 11, 0, 13, 0, 0, 0, 17, 0, 19, 0, 0, 0, 23, 0, 0, 0, 0, 0, 29, 0, 31, 0, 0, 0, 0, 0, 37, 0, 0, 0, 41, 0, 43, 0, 0, 0, 47, 0, 0, 0, 0, 0, 53, 0, 0, 0, 0, 0, 59, 0, 61, 0, 0, 0, 0, 0, 67, 0, 0, 0, 71, 0, 73, 0, 0, 0, 0, 0, 79, 0, 0, 0, 83, 0, 0, 0, 0, 0, 89, 0, 0, 0, 0, 0, 0, 0, 97, 0, 0, 0, 101, 0, 103, 0, 0, 0, 107, 0, 109, 0, 0, 0, 113, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 127, 0, 0, 0, 131, 0, 0, 0, 0, 0, 137, 0, 139, 0, 0, 0, 0, 0, 0, 0, 0, 0, 149, 0, 151, 0, 0, 0, 0, 0, 157, 0, 0, 0, 0, 0, 163, 0, 0, 0, 167, 0, 0, 0, 0, 0, 173, 0, 0, 0, 0, 0, 179, 0, 181, 0, 0, 0, 0, 0, 0, 0, 0, 0, 191, 0, 193, 0, 0, 0, 197, 0, 199, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 211, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 223, 0, 0, 0, 227, 0, 229, 0, 0, 0, 233, 0, 0, 0, 0, 0, 239, 0, 241, 0, 0, 0, 0, 0, 0, 0, 0, 0, 251, 0, 0, 0, 0, 0, 257, 0, 0, 0, 0, 0, 263, 0, 0, 0, 0, 0, 269, 0, 271, 0, 0, 0, 0, 0, 277, 0, 0, 0, 281, 0, 283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 293, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 307, 0, 0, 0, 311, 0, 313, 0, 0, 0, 317, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 331, 0, 0, 0, 0, 0, 337, 0, 0, 0, 0, 0, 0, 0, 0, 0, 347, 0, 349, 0, 0, 0, 353, 0, 0, 0, 0, 0, 359, 0, 0, 0, 0, 0, 0, 0, 367, 0, 0, 0, 0, 0, 373, 0, 0, 0, 0, 0, 379, 0, 0, 0, 383, 0, 0, 0, 0, 0, 389, 0, 0, 0, 0, 0, 0, 0, 397, 0, 0, 0, 401, 0, 0, 0, 0, 0, 0, 0, 409, 0, 0, 0, 0, 0, 0, 0, 0, 0, 419, 0, 421, 0, 0, 0, 0, 0, 0, 0, 0, 0, 431, 0, 433, 0, 0, 0, 0, 0, 439, 0, 0, 0, 443, 0, 0, 0, 0, 0, 449, 0, 0, 0, 0, 0, 0, 0, 457, 0, 0, 0, 461, 0, 463, 0, 0, 0, 467, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 479, 0, 0, 0, 0, 0, 0, 0, 487, 0, 0, 0, 491, 0, 0, 0, 0, 0, 0, 0, 499, 0, 0, 0, 503, 0, 0, 0, 0, 0, 509, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 521, 0, 523, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 541, 0, 0, 0, 0, 0, 547, 0, 0, 0, 0, 0, 0, 0, 0, 0, 557, 0, 0, 0, 0, 0, 563, 0, 0, 0, 0, 0, 569, 0, 571, 0, 0, 0, 0, 0, 577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 587, 0, 0, 0, 0, 0, 593, 0, 0, 0, 0, 0, 599, 0, 601, 0, 0, 0, 0, 0, 607, 0, 0, 0, 0, 0, 613, 0, 0, 0, 617, 0, 619, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 631, 0, 0, 0, 0, 0, 0, 0, 0, 0, 641, 0, 643, 0, 0, 0, 647, 0, 0, 0, 0, 0, 653, 0, 0, 0, 0, 0, 659, 0, 661, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 673, 0, 0, 0, 677, 0, 0, 0, 0, 0, 683, 0, 0, 0, 0, 0, 0, 0, 691, 0, 0, 0, 0, 0, 0, 0, 0, 0, 701, 0, 0, 0, 0, 0, 0, 0, 709, 0, 0, 0, 0, 0, 0, 0, 0, 0, 719, 0, 0, 0, 0, 0, 0, 0, 727, 0, 0, 0, 0, 0, 733, 0, 0, 0, 0, 0, 739, 0, 0, 0, 743, 0, 0, 0, 0, 0, 0, 0, 751, 0, 0, 0, 0, 0, 757, 0, 0, 0, 761, 0, 0, 0, 0, 0, 0, 0, 769, 0, 0, 0, 773, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 787, 0, 0, 0, 0, 0, 0, 0, 0, 0, 797, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 809, 0, 811, 0, 0, 0, 0, 0, 0, 0, 0, 0, 821, 0, 823, 0, 0, 0, 827, 0, 829, 0, 0, 0, 0, 0, 0, 0, 0, 0, 839, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 853, 0, 0, 0, 857, 0, 859, 0, 0, 0, 863, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 877, 0, 0, 0, 881, 0, 883, 0, 0, 0, 887, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 907, 0, 0, 0, 911, 0, 0, 0, 0, 0, 0, 0, 919, 0, 0, 0, 0, 0, 0, 0, 0, 0, 929, 0, 0, 0, 0, 0, 0, 0, 937, 0, 0, 0, 941, 0, 0, 0, 0, 0, 947, 0, 0, 0, 0, 0, 953, 0, 0, 0, 0, 0, 0, 0, 961, 0, 0, 0, 0, 0, 967, 0, 0, 0, 971, 0, 0, 0, 0, 0, 977, 0, 0, 0, 0, 0, 983, 0, 0, 0, 0, 0, 989, 0, 991, 0, 0, 0, 0, 0, 997, 0, 999]
def eratostenes(m): lista=[n for n in range(0,m+1)] lista[1]=0 for p in prime_range(sqrt(m)+1): for k in range(2, floor(m/p)+1): lista[k*p]=0 return lista 
       

Classic cripto...

Sagemath has some classic criptosystems...
S = AlphabeticStrings(); S 
        
Free alphabetic string monoid on A-Z
Free alphabetic string monoid on A-Z
E = VigenereCryptosystem(S,5) E 
        
Vigenere cryptosystem on Free alphabetic string monoid on A-Z of period
5
Vigenere cryptosystem on Free alphabetic string monoid on A-Z of period 5
K = E.random_key(); K 
        
VBONV
VBONV
L = E.inverse_key(K); L 
        
FZMNF
FZMNF
M = S("THISISASECRETMESSAGE") 
       
MyCipher=E(K) MyCipher 
        
Cipher on Free alphabetic string monoid on A-Z
Cipher on Free alphabetic string monoid on A-Z
hide=MyCipher(M) hide 
        
OIWFDNBGRXMFHZZNTOTZ
OIWFDNBGRXMFHZZNTOTZ
MyDecipher=E(L) 
       
MyDecipher(hide) 
        
THISISASECRETMESSAGE
THISISASECRETMESSAGE
HillCipher= HillCryptosystem(S,5); HillCipher 
        
Hill cryptosystem on Free alphabetic string monoid on A-Z of block
length 5
Hill cryptosystem on Free alphabetic string monoid on A-Z of block length 5
MyKey=HillCipher.random_key() MyKey 
        
[ 1 24 14  9 22]
[24 11  5  2  9]
[ 9 20 10 21  1]
[15  5  4  0 10]
[14  8 25  0 19]
[ 1 24 14  9 22]
[24 11  5  2  9]
[ 9 20 10 21  1]
[15  5  4  0 10]
[14  8 25  0 19]
MyKey.determinant() 
        
21
21
MyKey.inverse() 
        
[24 19  4  9 13]
[ 6 25  0 18 21]
[25  2 25 20 13]
[19 17 20 15 13]
[ 3 18 23  6 22]
[24 19  4  9 13]
[ 6 25  0 18 21]
[25  2 25 20 13]
[19 17 20 15 13]
[ 3 18 23  6 22]
M 
        
THISISASECRETMESSAGE
THISISASECRETMESSAGE
MyCipher=HillCipher(MyKey) MyCipher 
        
Hill cipher on Free alphabetic string monoid on A-Z of block length 5
Hill cipher on Free alphabetic string monoid on A-Z of block length 5
hide=MyCipher(M) hide 
        
RPDPPIWEUYAOYOBYQYQS
RPDPPIWEUYAOYOBYQYQS
MyInverseKey=HillCipher.inverse_key(MyKey) MyInverseKey 
        
[24 19  4  9 13]
[ 6 25  0 18 21]
[25  2 25 20 13]
[19 17 20 15 13]
[ 3 18 23  6 22]
[24 19  4  9 13]
[ 6 25  0 18 21]
[25  2 25 20 13]
[19 17 20 15 13]
[ 3 18 23  6 22]
MyKey*MyInverseKey 
        
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
Decipher=HillCipher(MyInverseKey) Decipher(hide) 
        
THISISASECRETMESSAGE
THISISASECRETMESSAGE

Elliptic curves

 

E=EllipticCurve([-5, 4])
P=E.plot()
P.show()

The "standard elliptic curve" has the form

y^2= x^3 + ax + b,

for some fixed a and b. (This is called a Weierstrass equation).

In characteristic 2 or 3 the form is slightly more complicated.

 
E=EllipticCurve([-5, 4]) P=E.plot() P.show() 
       
E 
        
Elliptic Curve defined by y^2 = x^3 - 5*x + 4 over Rational Field
Elliptic Curve defined by y^2 = x^3 - 5*x + 4 over Rational Field
E.discriminant() 
        
1088
1088
E=EllipticCurve([-1, 1]) P=E.plot() show(P,figsize=[8,4]) 
       
E.discriminant() 
        
-368
-368
E=EllipticCurve([-5, 4]); E E.plot() 
       
P=E([0,-2]) Q=E([1,0]) R=P+Q 
       
G1=plot(E) GP=plot(P, pointsize=50, rgbcolor='red') GQ=plot(Q, pointsize=50, rgbcolor='red') GmenosR=plot(-R, pointsize=50, rgbcolor='purple') GR=plot(R, pointsize=50, rgbcolor='green') G2R=plot(2*R,pointsize=50, rgbcolor='black') show(G1+GP+GQ+GR+GmenosR+G2R, figsize=[10,6], gridlines=True) 
       
E=EllipticCurve(GF(37),[1,0]) E.plot() 
       
E 
        
Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 37
Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 37

Diffie-Helmann with EC

Alice and Bob choose an e.c. E modulus p.
p=random_prime(100000000) R=IntegerModRing(p) 
       
p,R 
        
(85938953, Ring of integers modulo 85938953)
(85938953, Ring of integers modulo 85938953)
a=R.random_element() b=R.random_element() while 4*a^3+27*b^2==0: a=R.random_element() b=R.random_element() 
       

We now have an "authorized" e.c.
E=EllipticCurve(R,[a,b]); E 
        
Elliptic Curve defined by y^2 = x^3 + 60201846*x + 12864088 over Ring of
integers modulo 85938953
Elliptic Curve defined by y^2 = x^3 + 60201846*x + 12864088 over Ring of integers modulo 85938953

Alice and Bob pick a random point P in E.
P=E.random_point(); P 
        
(12808456 : 23729794 : 1)
(12808456 : 23729794 : 1)

Alice picks a and Bob picks b. 

Alice sends aP to Bob, and Bob sends bP to Alice. Now both can compute abP.
a=randint(2,order(E)) b=randint(2,order(E)) a, b 
        
(58323290, 55794135)
(58323290, 55794135)
A=a*P; B=b*P; A, B 
        
((16143726 : 52623811 : 1), (4152859 : 44091311 : 1))
((16143726 : 52623811 : 1), (4152859 : 44091311 : 1))
b*A 
        
(75728692 : 3973049 : 1)
(75728692 : 3973049 : 1)
a*B 
        
(75728692 : 3973049 : 1)
(75728692 : 3973049 : 1)