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 bools 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$.

Avatar
Aritanan Gruber
Assistant Professor

“See, if y’all haven’t the same feeling for this, I really don’t give a damn. If you ain’t feeling it, then dammit this ain’t for you!“
(desconheço a autoria; agradeço a indicação)