4750-A tour at SAGEMATH.sagews

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

x^2==4

True

x^3-2!=x^2+2

False

2<=3

True

3>2

True

a=x**4%3;
a

1

y=pi/2;
y

1/2*pi

n(y)

1.57079632679490

MyPi=numerical_approx(pi,prec=150)

MyPi

3.1415926535897932384626433832795028841971694

cos(MyPi)

-1.0000000000000000000000000000000000000000000

cos(pi)

-1

type(cos(pi))

<type 'sage.symbolic.expression.Expression'>

type(cos(MyPi))

<type 'sage.rings.real_mpfr.RealNumber'>

Changing types on-the-fly...

z=3

type(z)

<type 'sage.rings.integer.Integer'>

z=z+2/3; z

11/3

type(z)

<type 'sage.rings.rational.Rational'>

Changing basis...

num1=10; num2=7
[num1.str(2),num2.str(2)]

['1010', '111']

num1.str(2)+num2.str(2)

'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

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 (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, 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)]

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)

@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)

@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)

@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

Rings, fields and vector spaces

Z7=IntegerModRing(7)

a=Z7(3)

1/a

5

Z7.is_field()

True

Z8=IntegerModRing(8)

Z8.is_field()

False

F8=GF(2^3,'a');
F8

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]

F8.characteristic()

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'>

pol in R

True

pol.is_irreducible()

True

pol.factor()

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

f=x+1

type(f)

<class 'sage.rings.polynomial.polynomial_quotient_ring_element.PolynomialQuotientRingElement'>

f in S

True

type(pol)

<type 'sage.rings.polynomial.polynomial_gf2x.Polynomial_GF2X'>

1/f

x^2 + x

S.order()

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)]

g=S(x)

1/g

x^2 + 1

f*g

x^2 + x

f^3+g^(-2)

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)]

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

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]

corpo.<x>=minhachave[1]
a=minhachave[2]

mess=corpo(x^3+x+1)
mess

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]

M.echelon_form()

[ 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]

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]

transpose(M).nullity()

1

M.rref()

[ 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

K.modulus()

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

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]

A.det()

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]

A.det()

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)

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)

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)

a*x+b*y

1

Set $P$ as the set of prime numbers..

P=Primes(); P

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]

[n for n in range(2,40) if n in P]

[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]]

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)

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])

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

2^(n-1)%n

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"))' + '\n', '', 'single') 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' + '\n', '', 'single') 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

R(2)^(n-1)

825326789816013506777461377270418

2.powermod(n-1,n)

825326789816013506777461377270418

factor(n)

13 * 1263499 * 173420475484059478781201647

R(2)^(n-1)==R(1)

False

We may also compute the multicative inverse of an element, in case it exists.

1/R(322)

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]

lista[0]

0

floor(sqrt(1000))

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]

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

E = VigenereCryptosystem(S,5)
E

Vigenere cryptosystem on Free alphabetic string monoid on A-Z of period 5

K = E.random_key(); K

VBONV

L = E.inverse_key(K); L

FZMNF

M = S("THISISASECRETMESSAGE")

MyCipher=E(K)
MyCipher

Cipher on Free alphabetic string monoid on A-Z

hide=MyCipher(M)
hide

OIWFDNBGRXMFHZZNTOTZ

MyDecipher=E(L)

MyDecipher(hide)

THISISASECRETMESSAGE

HillCipher= HillCryptosystem(S,5);
HillCipher

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]

MyKey.determinant()

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]

M

THISISASECRETMESSAGE

MyCipher=HillCipher(MyKey)
MyCipher

Hill cipher on Free alphabetic string monoid on A-Z of block length 5

hide=MyCipher(M)
hide

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]

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]

Decipher=HillCipher(MyInverseKey)
Decipher(hide)

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

E.discriminant()

1088

E=EllipticCurve([-1, 1])
P=E.plot()
show(P,figsize=[8,4])

E.discriminant()

-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

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)

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

Alice and Bob pick a random point $P$ in $E$.

P=E.random_point(); P

(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)

A=a*P;
B=b*P;
A, B

((16143726 : 52623811 : 1), (4152859 : 44091311 : 1))

b*A

(75728692 : 3973049 : 1)

a*B

(75728692 : 3973049 : 1)