Introdução ao SageMath¶
Sage¶
O que é o Sage?¶
- Sage é um sistema de computação simbólica
- Oferece uma interface em Python para Álgebra, Grafos, Análise Numérica, e muito mais
- No Backend, usa sistemas como LAPACK, ATLAS, PARI/GP e M4RI
- Distribuído sob a licença GPLv3
Formalmente ( : $$ \text{Sage} > \left( \text{Magma} \cup \text{Maple} \cup \text{Mathematica} \cup \text{MATLAB} \right) $$
Há duas semanas dei uma monitoria básica sobre Sage, que pode ser acessada aqui: https://github.com/thalespaiva/monitoria-cripto/blob/master/intro-sage-ep1.ipynb
Mas o que é computação simbólica?
Sistemas de computação simbólica permitem trabalhar com expressões matemáticas:
x, y = var('x, y')
f(x, y) = x**2 + y**2
print(f)
f.derivative()
plot3d(f, (x, -1, 1), (y, -1, 1))
Importante. Quando queremos respostas numéricas, computação simbólica mais atrapalha do que ajuda.
93/13 + (130 - 93)/13
log(1201, 2)
Para obtermos a aproximação numérica de uma expressão, usamos a função n que toma como parâmetros uma expressão e, opcionalmente, o número de dígitos da aproximação.
n(93/13)
n(log(1201, 2), digits=30)
numerical_approx(log(1201, 2))
Para que usar Sage?¶
Vamos comparar código em C e Sage para fazer multiplicações de matrizes binárias.
Em C, a biblioteca mais rápida para fazer multiplicação binária eficiente é chamada M4RI.
Considere a implementação abaixo:
mzd_t* multiply_two_random_matrices(rci_t n) {
mzd_t* A = mzd_init(n, n);
mzd_t* B = mzd_init(n, n);
mzd_randomize(A);
mzd_randomize(B);
mzd_t* C = mzd_mul(NULL, A, B, NULL);
mzd_free(A);
mzd_free(B);
return C;
}
Agora o código em Sage que realiza exatemente a mesma computação:
def multiply_two_random_matrices(n):
A = random_matrix(GF(2), n, n)
B = random_matrix(GF(2), n, n)
return A*B
C = multiply_two_random_matrices(8)
print(C)
As duas implementações rodam em tempo similar pois Sage chama a rotina da M4RI para fazer a multiplicação!
Além disso, temos a vantagem de ter a mesma interface para vários sistemas bem diferentes entre si. Por exemplo: podemos usar $\mathbf{C}$ diretamente como matriz de adjacência para construir um grafo
G = sage.graphs.graph.Graph(C, format='adjacency_matrix')
plot(G)
Podemos obter informações importantes sobre o grafo gerado. Por exemplo, se o grafo é conexo ou não.
G.is_connected()
Por onde você começaria a implementar em C isso que acabamos de fazer tão facilmente?
Como usar o Sage?¶
O SageMath pode ser usado de 3 formas:
- Usando o Jupyter Notebook (essa é a forma sendo usada neste documento)
- Com um interpretador interativo (como Python):
- Por meio de scripts (como Python)
Para o EP, vocês deverão entregar um script .py.
Um dos melhores editores para isso é o Atom (https://atom.io/), porém qualquer editor de texto pode ser usado.
Para que seu script possa usar todas as funções do sage, ele deve conter from sage.all import * no topo.
script.py
from sage.all import *
def func1():
print('Esta é a função 1')
def func2():
print('Esta é a função 2')
if __name__ == "__main__":
func1() # Chama a função 1
func2() # Chama a função 2
Para rodar o arquivo script.py, use:
sage script.py
Para a depuração do seu script, é muito útil rodar o script dentro do interpretador interativo do Sage. Da mesma forma que fazemos isso no iPython, podemos rodar %run script.py dentro do Sage.
/home/user/$ sage
┌────────────────────────────────────────────────────────────────────┐
SageMath version 9.0, Release Date: 2020-01-01
Using Python 3.8.2. Type "help()" for help.
└────────────────────────────────────────────────────────────────────┘
sage: %run script.pyAritmética modular¶
Nesta seção, trataremos de $Z_n$, que é o anel de inteiros módulo $n$. Primeiro relembremos a definição de Anel.
Definição. Um anel $R$ é um conjunto munido de duas operações $+$ (adição) e $.$ (multiplicação), que, para todos $a, b, c \in R$, satisfazem:
- Associatividade de $+$ $$ (a + b) + c = a + ( b + c) $$
- Comutatividade de $+$ $$ a + b = b + a $$
Existência do identidade de $+$. Existe um elemento $0 \in R$ tal que $$ a + 0 = a $$
Existência da inversa aditiva. Existe um elemento $(-a) \in R$ tal que $$ a + (-a) = 0 $$
Associatividade de $.$ $$ (a b) c = a ( b c) $$
Existência da identidade de $.$
Existe um elemento $1 \in R$ tal que $$ a1 = 1a = a $$
- Distributividade de $.$ sobre $+$: $$ (a + b)c = ac + bc \\ c(a + b) = ca + cb $$
Importante. Note que num anel não é necessário que $.$ seja comutativa, e nem que exista a inversa multiplicativa de todos os elementos.
Em Sage, é muito simples trabalhar com $Z_n$, o anel de inteiros $\bmod n$.
n = 13
Zn = Integers(n)
print(Zn)
Zn.addition_table(names='elements')
Zn.multiplication_table(names='elements')
O Sage cuida para que as operações sejam feitas dentro deste anel.
10 + 9 # Operação fora de Zn
(10 + 9) % 13
10**20 # Operação fora de Zn
10**20 % 13
Zn(10) + Zn(9) # Operação em Zn
Zn(10)**20 # Operação em Zn
Note que, quando $\gcd(a, n) = 1$, existe inversa de $a \bmod n$, que pode é computada diretamente:
Zn(4)^-1
Zn(10)*Zn(4)
Além disso, para casos em que precisamos apenas da exponenciação ou inversão modular, sem que seja necessário fazer outras operações em $Z_n$, podemos usar as funções power_mod ou inverse_mod. Isso evita a sobrecarga de criar o $Z_n$, como fizemos anteriormente.
power_mod(10, 20, 13)
inverse_mod(4, 13)
Compare as respotas acima com as calculadas sobre elementos de $Z_n$ diretamente.
Funções importantes como o Máximo Divisor Comum, Algoritmo de Euclides estendido, e Teorema Chinês do Resto são implementados em Sage.
O MDC (Greatest Common Divisor - GCD) pode ser computado usando a função gcd
gcd(133, 48)
gcd(12, 18)
O algoritmo de Euclides estendido é feito pela função xgcd
res, u, v = xgcd(133, 48)
res
u*133 + v*48
Considere os seguintes resíduos para os quais queremos aplicar o Teorema Chinês do Resto (TCR). $$ \begin{cases} x \equiv 17 \mod 23 \\ x \equiv 91 \mod 111 \\ x \equiv 40 \mod 57 \\ \end{cases} $$
O Teorema Chinês do Resto (Chinese Remainder Theorem - CRT) pode ser usado para encontrar uma solução para as congruências usando a função crt
x = crt([17, 91, 40], [23, 111, 57])
x
x % 23, x % 111, x % 57
Matrizes¶
Como vimos no começo da aula, Sage é capaz de trabalhar com álgebra de matrizes.
Criando matrizes e vetores¶
Para criar uma matriz usando uma lista de listas, usamos a função matrix
A = matrix(
[
[1, 3, 2],
[2, 5, 1],
[7, 8, 3],
]
)
print(A)
Como não foi explicitado a que anel pertencem os elementos de A, o Sage supõe que são inteiros. Portanto, ao calcular A*A, a resposta é dada em inteiros.
A**10
Em criptografia, é muito comum trabalharmos com matrizes de inteiros $\bmod n$. Para isso, podemos fazer a inicialização da seguinte forma.
Zn = Integers(13) # Apenas para lembrar do que é Zn
An = matrix(Zn,
[
[1, 3, 2],
[2, 5, 1],
[7, 8, 3],
]
)
print(An)
Agora, o produto An*An e soma An + An são calculado em Zn
An**10
An + An
Muitas vezes, em criptografia, precisamos gerar matrizes aleatórias. Isto pode ser feito usando a função random_matrix, que recebe como parâmetro um Anel de base, o número de linhas e o número de colunas.
random_matrix(Zn, 3, 3)
Operações sobre matrizes¶
O método inverse calcula inversa de uma matriz quadrada, quando esta é não singular.
An_inv = An.inverse()
print(An_inv)
An_inv*An
An**(-1)
O determinante pode ser calculado usando o método determinant
A.determinant()
O método echelon_form devolve a matriz escalonada. Não confunda com o método echelonize, que escalona a matriz, alterando-a.
B = random_matrix(Zn, 5, 10)
print('Matriz B')
print(B)
print('Matriz B escalonada')
print(B.echelon_form())
Posto e nulidade:
B.rank()
B.nullity()
Resolvendo sistemas lineares¶
Para resolver um sistema da forma $\mathbf{D} \mathbf{x} = \mathbf{b}$, usamos o método solve_right.
D = random_matrix(Zn, 4, 4)
b = random_vector(Zn, 4)
print('D = ')
print(D)
print('b = ', b)
x = D.solve_right(b)
print(x)
D*x == b
Corpos¶
Definição. Um Corpo K é um conjunto munido de duas operações $+$ e $.$ tais que:
- $(K, +, .)$ é um anel
A operação $.$ é comutativa. Ou seja, para todos $a, b \in K$: $$ a b = b a, $$
Para todo elemento $a$ não nulo de $K$, existe o seu inverso multiplicativo $a^{-1} \in K$ tal que: $$ aa^{-1} = a^{-1}a = 1 $$
Observação. Corpo é uma abstração do conjunto dos Racionais, onde todos os elementos diferentes de 0 têm inversa.
Em criptografia, tipicamente estamos interessados em corpos finitos, pois os elementos devem ser representados num número finito e prefixado de bits.
Teorema. Se $p$ é número primo, então $Z_p$ é um corpo finito.
Demontração. Sabemos que $(Z_p, +, .)$ é um anel, então só precisamos mostrar que valem as propriedades 2 e 3. Sabemos também que a multiplicação modular é comutativa, o que prova que a propriedade 2 vale. Além disso, como $p$ é primo, então $\gcd(a, p) = 1$ para todo $a \in Z_p - \{0\}$. Portanto existe a inversa de $a \mod p$, e vale a propriedade 3.
Em sage, para determinar se um determinado anel é também um corpo, podemos usar o método is_field
Integers(13).is_field()
Integers(10).is_field()
Em geral, ao programar, devemos ser o mais explícitos possível. Portanto, ao declarar um conjunto que sabemos que é um corpo, é preferível usar o construtor GF (de Galois Field)
K = GF(13)
Este construtor quebra quando passamos um tamanho incompatível com um corpo finito
try:
GF(10)
except Exception as e:
print(e)
Há corpos de tamanho não-primo. Em particular o tamanho de um corpo é sempre $q = p^k$ para algum primo $p$ e algum inteiro positivo $k \geq 1$.
Porém, pelo modo como é feita a multiplicação em corpos finitos de tamanhos não primos, seus elementos são melhor representados por polinômios e não por inteiros.
GF(7**2, 't').random_element()
Integers(7**2).is_field()
GF(7**2).random_element()
GF(7).random_element()
Para entender como chegar nisso, primeiro vejamos como construir anéis polinômiais.
Anéis Polinomiais¶
Um anel polinomial $K[x]$ nada mais é que o conjunto de polinômios em $x$ com coeficientes em $K$. Tipicamente, estamos interessados em $K$ sendo um corpo, mas a construção vale para $K$ sendo apenas um anel. Para ver que $K[x]$ é um anel, basta usar o fato de que $K$ é um corpo (e portanto um anel) e usar as regras de adição e multiplicações de polinômios.
Para criar um anel polinomial R sobre a variavel x e com coeficientes em GF(13), em Sage, fazemos:
R.<x> = GF(13)[]
R
R2 = GF(13).polynomial_ring('x2')
R2
v = R.random_element(4) # 10 é o grau do elemento aleatório
print(v)
Note que $GF(n)[x]$ é infinito.
R.cardinality()
Para contruir um anel polinomial finito, podemos fazer o análogo do que fizemos para construir $Z_n$ a partir de $Z$: pegar os resíduos quando dividimos por um elemento $p_1 \in R$. O anel construído é denotado por $R/p_1$
p1 = x^2 + x
R_mod_p1 = R.quotient(p1)
Note que o número de resíduos quando dividimos por um polinômio de grau 2 é dado por $$ ( \text{possibilidades para coeficiente de grau 0}) \times (\text{possibilidades para coeficiente de grau 1}) = n^2 $$
R_mod_p1.cardinality()
Porém note que o anel quociente construído não é um corpo:
R_mod_p1.is_field()
O motivo por que $R/p_1$ não é um corpo é análogo ao motivo de $Z_{10}$ não ser um corpo: o módulo não é primo, e há elementos $a$ em $R/p_1$ tais que $\gcd(a, p_1) \neq 1$
gcd(x + 1, x^2 + x)
Definição. Um polinômio $m \in GF(q)[x]$ é primo (ou irredutível) quando ele não pode ser fatorado em dois polinômios não-constantes de $GF(q)[x]$.
No caso anterior, $p_1= x^2 + x$ claramente pode ser fatorado em $x(x + 1)$.
Observação. Para mostrar que um polinômio não é primo, basta mostrar uma fatoração. Agora, para mostrar que um polinômio é primo não é tão fácil...
Por sorte, o Sage vem com um método que devolve um polinômio irredutível num anel polinomial, dado um grau.
m = R.irreducible_element(2)
print(m)
Podemos agora fazer a mesma construção anterior só que com $m$ irredutível no lugar de $p_1$.
R_mod_m = R.quotient(m)
Neste caso, $R/m$ será um corpo finito de $13^2 = 169$ elementos.
R_mod_m.is_field()
R_mod_m.cardinality()
Note que $R/m$ e $R/p_1$ contam com exatamente os mesmos polinômios! A diferença ocorre na multiplicação de elementos, que, quando feita $\bmod m$, mantém a propriedade da existência do inverso.
Importante. Se tomarmos outro polinômio irredutível $m_2$ do mesmo grau de $m$, os corpos finitos $R/m$ e $R/m_2$ são isomórficos. Há um teorema importante que diz que todos os corpos finitos de mesmos tamanhos são isomórficos.
Para ver isso na prática, considere.
P.<t1> = GF(2)[]
m1 = P.irreducible_element(3, algorithm='random')
m2 = P.irreducible_element(3, algorithm='random')
print("Queremos m1 != m2")
print(m1)
print(m2)
P.quotient(m1, 'i').multiplication_table()
P.quotient(m2, 'j').multiplication_table()
Observação. É possível achar uma bijeção entre os polinômios dos dois corpos que transforma uma tabela na outra. (Não que seja uma coisa muito fácil de fazer agora, mas acho que a ideia do isomorfismo fica bem ilustrada pelas tabelas.)
O Corpo finito $GF(2^m)$¶
Vimos que todos os corpos finitos são isomórficos a $GF(p^k)$ para algum primo $p$ e algum inteiro $k \geq 1$. Vimos também que para usar qualquer um desses corpos finitos em Sage, usamos o construtor GF, não sendo necessário criar anéis polinomiais e selecionar um polinômio irredutível.
Em criptografia é comum usarmos corpos finitos de característica 2, isto é, $p = 2$, pois algumas operações podem ser simplificadas usando operações de bits.
Uma aplicação interessante é que podemos definir uma operação de multiplicação sobre bytes.
Considere por exemplo o $GF(2^8)$ na representação usada pelo AES.
t = GF(2).polynomial_ring().gen()
F8 = GF(2^8, 't', modulus=t^8 + t^4 + t^3 + t + 1)
b1 = F8.random_element()
b2 = F8.random_element()
print(b1)
print(b2)
Para obter a representação de b1 e b2 como sequências de bits, podemos fazer:
print(list(b1.polynomial()))
print(list(b2.polynomial()))
Note como a soma de b1 e b2 é simplesmente o XOR bit a bit
print(list((b1 + b2).polynomial()))
E agora temos uma definição para o produto de dois bytes b1*b2
b1*b2
print(list((b1*b2).polynomial()))
Para pegar a representação em bytes de um elemento de F8, basta usarmos o método integer_representation
n1 = b1.integer_representation()
print(n1)
Para o processo inverso, não há uma função pronta. Podemos usar a função bin para implementar a conversão de inteiros para elementos de F8
def int_to_F8_poly(n):
t = F8.gen()
if not (0 <= n < 2**8) or not (0 <= n < 2**8):
raise ValueError('Deve-se obedecer à condição 0 <= int1, int2 < 256')
coeffs = reversed(bin(n)[2:])
return sum(GF(2)(v)*(t**i) for i, v in enumerate(coeffs))
int_to_F8_poly(23)
Teste de sanidade:
int_to_F8_poly(b1.integer_representation()) == b1
Para explicitar as operações que queremos implementar, podemos usar as seguintes funções bytes_sum e bytes_product:
def bytes_sum(n1, n2):
'''
n1: Inteiro >= 0 e < 256
n2: Inteiro >= 0 e < 256
'''
b1 = int_to_F8_poly(n1)
b2 = int_to_F8_poly(n2)
return (b1 + b2).integer_representation()
def bytes_product(n1, n2):
'''
n1: Inteiro >= 0 e < 256
n2: Inteiro >= 0 e < 256
'''
b1 = int_to_F8_poly(n1)
b2 = int_to_F8_poly(n2)
return (b1 * b2).integer_representation()
hex(bytes_product(int('53', 16), int('ca', 16))) # Deve dar 0x01
hex(bytes_product(int('45', 16), int('02', 16))) # Deve dar 0x8a
hex(bytes_product(int('c5', 16), int('02', 16))) # Deve dar 0x91
Importante. Este modo de implementar a soma e o produto não é o mais eficiente. Em geral, as operações são implementadas usando operações de bits, que são muito mais rápidas.
Curvas Elípticas¶
O que são curvas Elípticas?¶
Definição. Seja K. Uma curva elíptica é uma curva que pode ser descrita na forma: $$ y^2 + a_1 xy + a_3 y = x^3 + a_2 x^2 + a_4 x + a_6, $$ onde $a_1, a_2, a_3, a_4, a_6 \in K$ tais que $\Delta(a_1, a_2, a_3, a_4, a_6) \neq 0$.
Esta forma de definir a curva se chama Forma longa de Weierstrass.
Observação. Omiti intencionalmente a forma como o discriminante $\Delta$ é calculado para evitar deixar o texto mais pesado. Ele serve para garantir que a curva não é degenerada.
Importante. Esta é a forma mais geral para uma curva elíptica, porém ela pode ser simplificada de acordo a característica do corpo $K$.
Observação. Curvas elípticas não têm nada a ver com elipses.
E1 = EllipticCurve([0, -2, 0, 2, 10])
E1
plot(E1, xmin=-10, xmax=20)
Curvas elípticas definem um grupo¶
Primeiro, vejamos a definição de grupo. Agora que já vimos Anéis e Corpos, Grupos são bem tranquilos.
Definição. Um Grupo é um conjunto $G$ munido de uma operação $+$, tal que, para todos $a, b \in G$ vale que:
- $a + b \in G$
- $(a +b) + c = a + (b + c)$
- Existe $0 \in G$ tal que $a + 0 = 0 + a = a$
- Existe $-a \in G$ tal que $a + (-a) = 0$
Quando a operação $+$ é comutativa, o grupo é chamado de abeliano.
Queremos mostrar que uma curva elíptica pode definir um grupo abeliano. Para isso, precisamos definir uma operação de adição e mostrar que ela satisfaz as propriedades de grupos.
P = E1.an_element()*20
Q = E1.an_element()*30
Soma de dois pontos diferentes¶
(plot(E1, xmin=-10, xmax=20) +
plot(P, color='red', size=30) +
plot(Q, color='brown', size=30) +
plot(P + Q, color='green', size=30) +
line2d([P.xy(), (-P -Q).xy()], color='black') +
line2d([(-P -Q).xy(), (P+Q).xy()], color='black')
)
Soma de dois pontos iguais¶
(plot(E1, xmin=-10, xmax=10) +
plot(P, color='red', size=30) +
plot(P + P, color='green', size=30) +
line2d([P.xy(), (-P -P).xy()], color='black') +
line2d([(-P -P).xy(), (P + P).xy()], color='black')
)
Oposto de um ponto¶
(plot(E1, xmin=-10, xmax=10) +
plot(Q, color='red', size=30) +
plot(-Q, color='green', size=30) +
line2d([Q.xy(), (-Q).xy()], color='black')
)
O elemento neutro¶
Para manter a operação de soma consistente, é necessário adicionar um novo elemento, chamado ponto no infinito que é o ponto de encontro entre todas as retas verticais. (pausa dramática)
(plot(E1, xmin=-10, xmax=10) +
plot(Q, color='red', size=30) +
plot(-Q, color='green', size=30) +
line2d([Q.xy(), (-Q).xy()], color='black') +
line2d([(Q*15).xy(), (-Q*15).xy()], color='black')
)
Observação mais importante¶
Para calcular a soma de dois pontos, só são usadas operações de multiplicações e divisões de coordenadas, que são elementos do corpo $K$.
Logo, se $P$ tem coordenadas racionais, então $2P, 3P, 4P, \ldots$ todos terão coordenadas racionais.
Em particular, se $F$ é um corpo qualquer, e $P$ tem coordenadas em $F$, então $2P, 3P, 4P, \ldots$ todos terão coordenadas em $F$. Para criptografia, tipicamente usamos como $F$ um corpo finito.
Curvas Elípticas sobre corpos finitos¶
plot(E1, -30, 10) + plot(E1.change_ring(GF(17)))
poly = E1.defining_polynomial()
poly
poly.substitute(x=3, z=1)
poly.substitute(x=9, z=1)
Logaritmo discreto em curvas elípticas¶
A dificuldade aqui vêm de que múlitplos de um ponto $P$ parecem aleatórios.
E1_307 = E1.change_ring(GF(307))
g = E1_307.an_element()
Ps = []
for i in range(200):
Ps.append(scatter_plot([g*j for j in range(0, i)]))
a = animate(Ps)
a
Definição. Dados $P$ e $Q = kP$ pontos de uma curva elíptica construída sobre o corpo finito $F$. O problema do logaritmo discreto (em Curvas Elípticas) pede para encontrar $k$ dados $P$ e $Q$.
Curvas Elípticas em $GF(2^m)$¶
Podemos usar definir uma curva usando $F = GF(2^m)$. Como $F$ tem característica 2, as equações que definem as curvas elípticas podem ser simplificadas para a forma:
$$ y^2 + cy = x^3 + ax + b $$E2 = EllipticCurve(GF(2^10), [0, 0, 1, 2, 10])
E2
len(E2.rational_points())
P = E2.an_element()
Px, Py = P.xy()
Px.integer_representation(), Py.integer_representation()
points = []
for p in E2:
if p.is_zero():
continue
px, py = p.xy()
points.append([px.integer_representation(), py.integer_representation()])
scatter_plot(points)
