BCM0505-15 - Processamento da Informação
Aula 04 – 28/02
Fatoração de Inteiros
Problema 4.1: Dado um inteiro $n\geq 2$, determine primos $p_1,p_2,\ldots,p_k$, não necessariamente distintos, tais que $$n = p_1 p_2 \cdots p_k.$$
Sabemos, pelo teorema fundamental da aritmética, que todo inteiro positivo admite uma única fatoração em $\mathbb{Z}$. A ideia ingênua consiste em dividir $n$ por $2$ o maior número de vezes possível, pegar o resultado e fazer o mesmo por $3$, por $4$, e assim sucessivamente.
Claramente, $4$ não é um primo. No entanto, para não termos de determinar uma coleção de primos à priori — faremos isto no Problema 4.2 — podemos tentar dividir $n$ por $4$. Seguindo a ordem acima, $n$ certamente não será divisível por $4$ (ou qualquer outro não primo) pois já foi dividido por $2$ (ou por primos menores). O código abaixo implementa esta ideia.
# recebe : um inteiro n >= 2
# devolve: uma lista de primos [p1,p2,...,pk] tal que n == p1*p2*...*pk
def factoring (n: int) -> [int]:
factors = []
k = 2
while n >= 2:
if n % k == 0:
factors.append (k)
n //= k
else:
k += 1
return factors
Nota: O problema de fatoração de inteiros em primos é altamente relevante para alguns esquemas de segurança criptográfica. Por exemplo, um algoritmo eficiente de fatoração seria capaz de quebrar o protocolo RSA. Até o momento, não se sabe de um algoritmo clássico eficiente para tal tarefa — nosso algoritmo acima certamente não é um, pois requer tempo exponencial no comprimento da representação binária de $n$. Peter Shor introduziu em 1994 um algoritmo quântico que requer tempo apenas polinomial no comprimento da representação! No entanto, a construção de um computador quântico capaz de executar o algoritmo de Shor ainda está distante de ocorrer — e pode nem ser possível.
Crivo de Erathostenes.
Problema 4.2: Dado um inteiro $n\geq 0$, determine todos os primos menores ou iguais a $n$.
Podemos resolver o problema utilizando o teste de primalidade (is_prime()
) construído na aula passada, mas vamos fazer algo melhor.
Vamos utilizar uma lista $q[0\,..n]$ de bool
s para indicar se o inteiro $j$, associado à entrada $j$ de $q$, é primo ou não. Observe que isto não faz sentido para os índices $0$ e $1$ e vamos ignorar essas posições — é mais econômico fazer isso do que subtrair $2$ unidades do índice a cada acesso a $q$.
Inicialmente, todas as entradas de $q$ são True
. À medida que descobrimos que um valor $j$ não é primo, faremos q[j] = False
. No final, somente primos conterão True
em suas correspondentes entradas.
Resta especificarmos o principal: e como descobrimos que um valor $j$ não é primo sem utilizar is_prime()
?
Da forma mais simples possível. Pegamos o valor $2$ e marcamos todos os seus múltiplos como não primos. Depois, pegamos o valor $3$ e marcamos todos os seus múltiplos como não primos. O valor $4$ já foi marcado como não primo, logo avançamos para o próximo valor, o $5$. Fazemos o mesmo que foi feito para o $2$ e $3$, marcando todos os múltiplos de $5$ como não primos. Continuamos este processo até atingirmos $\lfloor\sqrt{n}\rfloor+1$, pois todos os não primos entre ele e $n$ já terão sido marcados.
O código abaixo implementa essa ideia.
# recebe : um inteiro n >= 0.
# devolve: uma lista com todos os primos <= n.
def erathostenes_sieve (n: int) -> [int]:
q = [True for i in range(n+1)]
k = 2
while k*k < n+1:
for i in range (k*k,n+1,k): q[i] = False
k += 1
while k*k <= n and not q[k]: k += 1
p = []
for i in range (2,n+1):
if q[i]: p.append (i)
return p
Nota: Podemos pegar a lista de primos gerada por erathostenes_sieve()
e revisitar o Problema 4.1.
# recebe : um inteiro n >= 2 e uma lista Q de primos menores ou iguais a \sqrt{n}
# devolve: uma lista de primos [p1,p2,...,pk] tal que n == p1*p2*...*pk
def factoring_with_primes (n: int, Q: [int]) -> [int]:
factors = []
k = 0
while n >= 2:
if n % Q[k] == 0:
factors.append (Q[k])
n //= Q[k]
else:
k += 1
return factors
O ganho de performance, no entanto, é irrisório frente ao RSA.
Derivação Numérica
Problema 4.3: Dados uma função diferenciável $f$ em um intervalo $I=(a,b)\subset\mathbb{R}$ e um ponto $x\in\ I$, calcule uma aproximação para a derivada de $f$ em $x$.
Vamos determinar $p\in\mathbb{R}$ tal que $$f'(x) = \frac{d}{dx}f(x) \approx px.$$ Isto é, $px$ é aproximadamente tangente a $f$ em $x$.
Usamos o método quociente de diferença simétrica: $$p = \frac{f(x-h) + f(x+h)}{2h},$$ para um valor positivo pequeno de $h$.
# recebe : uma função f diferenciável em I=[a,b] e um ponto x em I
# devolve: p e q tais que p*x+q é uma aproximação à derivada de f em x
def differentiate (f, x: float) -> float:
h = 0.001
return (f(x-h) + f(x+h))/(2*h)
Integração Numérica
Problema 4.4: Dados um intervalo $I=[a,b]\subset\mathbb{R}$, um número $n\geq 1$ de divisões para $I$ e uma função $f$ integrável em $I$, calcule uma aproximação para $$F(a,b)=\int_a^b f(x)dx,$$ a integral definida de Riemann de $f$ em $I$.
Observe que a função $f$ a ser integrada é um parâmetro do problema. Python permite passarmos funções como parâmetros da mesma forma que passamos escalares ou listas. Neste caso, qualquer função que receba e devolva um float
, como
def function (x: float) -> float:
...
return y
pode ser utilizada (desde que seja integrável em $I$).
Vamos apresentar três algoritmos de integração numérica para (4.4).
O primeiro, chamado método do ponto médio, consiste em:
- particionar $I$ em $[a = x_0 < x_1 < … < x_n = b]$ tal que $|x_{i}-x_{i+1}|=(b-a)/n$,
- para cada $i\in\{0,1,\ldots,n-1\}$, calcular $$M(x_{i},x_{i+1}) = (x_{i+1}-x_i)\cdot f\left(\frac{x_i+x_{i+1}}{2}\right),$$
- somar os termos $M(x_{i},x_{i+1})$, de forma que $$\begin{align*} F(a,b) &\approx M(x_0,x_1) + M(x_1,x_2) + \cdots + M(x_{n-1},x_n) \\ &= \frac{b-a}{n} \left( f\left(\frac{x_0+x_1}{2}\right) + f\left(\frac{x_1+x_2}{2}\right) + \cdots + f\left(\frac{x_{n-1}+x_n}{2}\right) \right). \end{align*}$$
# Integração numérica: método do ponto médio.
#
# Recebe: uma função real f integrável no intervalo I=[a,b],
# as extremidades a e b de I,
# n >= 1 tal que [a = x_0 < x_1 < ... < x_n = b].
#
# Para cada intervalo [x,y], M(x,y) = (y-x)*f((x+y)/2).
#
# Devolve:
# integral_a^b f(x)dx = M(x_0,x_1) + M(x_1,x_2) + ... + M(x_{n-1},x_n)
# = (b-a)/n*(f((x_0+x_1)/2) + f((x_1+x_2)/2) + ... + f((x_{n-1}+x_n)/2)).
#
def int_midpoint (f, a: float, b: float, n: int) -> float:
h, s = (b-a)/n, 0.0
for i in range(n):
s += f((2*(a+i*h)+h)/2)
return h*s
O segundo, chamado método dos trapézios, consiste em:
- particionar $I$ em $[a = x_0 < x_1 < … < x_n = b]$ tal que $|x_{i}-x_{i+1}|=(b-a)/n$,
- para cada $i\in\{0,1,\ldots,n-1\}$, calcular $$T(x_{i},x_{i+1}) = \frac{x_{i+1}-x_i}{2}\cdot \Big(f(x_i)+f(x_{i+1})\Big),$$
- somar os termos $T(x_{i},x_{i+1})$, de forma que $$\begin{align*} F(a,b) &\approx T(x_0,x_1) + T(x_1,x_2) + \cdots + T(x_{n-1},x_n) \\ &= \frac{b-a}{n} \left(\frac{1}{2}f(x_0) + f(x_1) + \cdots + f(x_{n-1}) + \frac{1}{2}f(x_n)\right). \end{align*}$$
# Integração numérica: método dos trapézios.
#
# Recebe: uma função real f integrável no intervalo I=[a,b],
# as extremidades a e b de I,
# n >= 1 tal que [a = x_0 < x_1 < ... < x_n = b].
#
# Para cada intervalo [x,y], T(x,y) = (y-x)/2*(f(x)+f(y)).
#
# Devolve:
# integral_a^b f(x)dx = T(x_0,x_1) + T(x_1,x_2) + ... + T(x_{n-1},x_n)
# = ((b-a)/n)*(1/2*f(x_1) + f(x_2) + ... + f(x_{n-1}) + 1/2*f(x_n)).
#
def int_trapezoid (f, a: float, b: float, n: int) -> float:
h, s = (b-a)/n, 1/2*(f(a) + f(b))
for i in range(1,n):
s += f(a+i*h)
return h*s
O terceiro, chamado método de Simpson, é uma combinação dos dois anteriores e consiste em:
- particionar $I$ em $[a = x_0 < x_1 < … < x_n = b]$ tal que $|x_{i}-x_{i+1}|=(b-a)/n$,
- para cada $i\in\{0,1,\ldots,n-1\}$, calcular $$\begin{align*} S(x_{i},x_{i+1}) &= \frac{2}{3}M(x_{i},x_{i+1}) + \frac{1}{3}T(x_{i},x_{i+1}) \\ &= \frac{x_{i+1}-x_i}{6}\cdot \left(f(x_i) + f\left(\frac{x_i+x_{i+1}}{2}\right) + f(x_{i+1})\right), \end{align*}$$
- somar os termos $T(x_{i},x_{i+1})$, de forma que $$\begin{align*} F(a,b) &\approx S(x_0,x_1) + S(x_1,x_2) + \cdots + S(x_{n-1},x_n) \\ &= \frac{b-a}{n} \left(\frac{1}{3}f(x_0) + \frac{4}{3}f(x_1) + \frac{2}{3}f(x_2) + \frac{4}{3}f(x_3) + \cdots + \frac{2}{3}f(x_{n-2}) + \frac{4}{3}f(x_{n-1}) + \frac{1}{3}f(x_n)\right). \end{align*}$$
# Integração numérica: método de Simpson.
#
# Recebe: uma função real f integrável no intervalo I=[a,b],
# as extremidades a e b de I,
# n >= 1 tal que [a = x_0 < x_1 < ... < x_n = b].
#
# Para cada intervalo [x,y], S(x,y) = 2/3*M(x,y)+1/3*T(x,y)
# = (y-x)/6*(f(x) + f((x+y)/2) + f(y)).
#
# Devolve:
# integral_a^b f(x)dx = S(x_0,x_1) + S(x_1,x_2) + ... + S(x_{n-1},x_n)
# = ((b-a)/n)*(1/3*f(x_0) + 4/3*f(x_1) + 2/3*f(x_2) + 4/3*f(x_3) + ...
# + 2/3*f(x_{n-2}) + 4/3*f(x_{n-1}) + 1/3*f(x_n)).
#
def int_simpson (f, a: float, b: float, n: int) -> float:
return 2/3*int_midpoint (f,a,b,n) + 1/3*int_trapezoid (f,a,b,n)
É instrutivo codificarmos o método diretamente, sem depender das funções anteriores.
def int_simpson (f, a: float, b: float, n: int) -> float:
h, s = (b-a)/n, 1/3*(f(a) + f(b))
for i in range(1,n):
c = 2*(i%2) + 2
s += c/3*f(a+i*h)
return h*s
Nota: O método de Simpson é exato para funções quadráticas quando $n$ é par e tal que $|x_{i-1}-x_i|$ é “suficientemente” pequeno.
Alguns exemplos de utilização.
(I) O código abaixo,
import math
def flin (x): return x
def fsqr (x): return x*x
a, b = -1, 3
for n in [2, 20, 200, 2000]:
print ("flin: ", int_midpoint (flin, a, b, n), int_trapezoid (flin, a, b, n), int_simpson (flin, a, b, n))
print ("fsqr: ", int_midpoint (fsqr, a, b, n), int_trapezoid (fsqr, a, b, n), int_simpson (fsqr, a, b, n))
print()
print ("sqrt: ", int_midpoint (math.sqrt, 0, 4, 20), int_trapezoid (math.sqrt, 0, 4, 20), int_simpson (math.sqrt, 0, 4, 20))
produz a saída:
flin: 4.0 4.0 4.0
fsqr: 8.0 12.0 9.333333333333332
flin: 4.0 4.000000000000001 4.000000000000001
fsqr: 9.320000000000004 9.360000000000001 9.333333333333334
flin: 4.000000000000001 4.0 4.0
fsqr: 9.3332 9.3336 9.333333333333336
flin: 3.9999999999999996 4.000000000000001 4.000000000000002
fsqr: 9.333331999999999 9.333335999999997 9.33333333333333
sqrt: 5.338362719363362 5.315572731413136 5.326072063999334
A título de curiosidade, vamos calcular os valores exatos (são de cálculo simples neste caso). Temos que $$\int_{-1}^3 x\,dx = \frac{x^2}{2}\bigg|_{-1}^3 = \frac{9}{2} - \frac{1}{2} = 4,$$ que $$\int_{-1}^3 x^2\,dx = \frac{x^3}{3}\bigg|_{-1}^3 = 9 + \frac{1}{3} = \frac{28}{3} = 9.\overline{3},$$ e que $$\int_{0}^4 x^{1/2}\,dx = \frac{2}{3}x^{3/2}\bigg|_{0}^4 = \frac{2}{3}8 = \frac{16}{3} = 5.\overline{3}.$$
Observe que alguns erros de precisão (ponto flutuante) podem ocorrer e acumular durante as contas.
(II) A distribuição de probabilidade Gaussiana (ou Normal) $\Phi_{0,1}$, com média $0$ e desvio padrão $1$, possui $$\phi(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$$ como função densidade de probabilidade. Isto é, $$\Phi_{0,1}(y) = \int_{-\infty}^y \phi(x)dx.$$ Se $X\sim\Phi_{0,1}$ é uma variável aleatória (distribuída de acordo com tal Gaussiana), temos para $a\leq b$ que $$\mathbb{P}\mathsf{r}[a\leq X\leq b] = \Phi_{0,1}(b) - \Phi_{0,1}(a) = \int_a^b \phi(x)dx.$$
import math
def fgaussian (x):
return math.exp(-x*x/2) / math.sqrt(2*math.pi)
for j in [10, 100, 1000, 10000, 100000]:
print (int_simpson (fgaussian, -3, 3, j))
Os valores impressos são:
0.9971953090849658
0.9973001924543536
0.9973002039355907
0.9973002039367428
0.9973002039367469
Pequenas Funções Anônimas (lambda) *
Esta seção é opcional (e não foi feita em sala).
As funções que utilizamos como exemplo nos métodos de derivação e integração podem ser consideradas “simples.” Nestes casos, e supondo que estas funções não precisam ser utilizadas em outra ocasião — para outro cálculo por exemplo — suas definições explícitas podem ser chatas e inconvenientes.
Python possui alguns mecanismos de programação funcional, dentre eles, a possibilidade de definirmos pequenas funções anonimamente. Chamadas de expressões lambda, estas funções são definidas como
lambda [lista_de_parâmetros]: expressão
e equivalem a
def <lambda>(parâmetros):
return expressão
O código
def foo (x):
return x**4 * math.log (x + math.sqrt (x*x+1))
print (int_simpson (foo, -2, 4, 200))
pode ser re-escrito como
print (int_simpson (lambda x: x**4 * math.log (x + math.sqrt (x*x+1)), -2, 4, 200))
Nota: Como lambda
define uma expressão, é possível utilizá-la de diversas formas, como
foo = lambda x: x**4 * math.log (x + math.sqrt (x*x+1))
print (foo(3))
ou
def bar (y: int) -> lambda:
return lambda x: x**4 * math.log (x + math.sqrt (x*x+y))
print (bar(1)(3) == foo(3), bar(2)(3) == foo(3)) # imprime True e False
A função bar()
cria uma função anônima diferente para cada valor de $y$.