MAC0122  Princípios de Desenvolvimento de Algoritmos

Máximo divisor comum

"Analysis of algorithms is the name the author likes
to use to describe investigations such as this. The general idea is
to take a particular algorithm and to determine its quantitative behaviour; occasionally
we also study wheher or not an algorithm is optimal in some sense.

—Donald E. Knuth, The Art of Computer Programming: Fundamental Algorithms.

  Está página é baseada na seção 2 do capítulo 31 do CLRS e capítulo 1 do

D.E. Knuth,
The art of computer programming, vol 1: Fundamental algorithms,
Addison-Wesley, Reading, Mass., 1968, QA758 K74a

 

  Suponha que m e n são números inteiros. Dizemos que n divide m se m % n == 0. Note que,

m % n == 0   se e somente se   m == dn para algum inteiro d.
As vezes escrevemos n\m como abreviação de "n divide m.

O máximo divisor comum de dois números inteiros m e n, denotado por mdc(m,n), é o maior número inteiro que divide ambos. Assim,

mdc(m,n)  =  max{ d  |  d divide m e k divide n }.
É claro que se m == n == 0, então o máximo não faz sentido. Logo, estaremos sempre supondo que pelo menos um dos números dados é diferente de zero. Nossa atual preocupação é com o seguinte problema
PROBLEMA:  Dados número inteiros m e n determinar mdc(m,n).

  Por simplicidade, na discussão a seguir estamos supondo que m >= 0 e n >= 0. Isto não é nenhum pecado, já que mdc(m,n) == mdc(abs(m),abs(n))

Vamos examinar dois algoritmos para o problema: um óbvio e lento e outro sofisticado e bem mais rápido.

 

Algoritmo café-com-leite

  Vamos começar com um algoritmo simples e óbvio: 

 /*
  * Ao receber dois números inteiros não-negativos m e n, não
  * ambos nulos, está função devolve o máximo divisor de m e n.
  */  
  long
  cafe_com_leite(long m, long n)
  {
    long d;

    if (n == 0) return m;
    if (m == 0) return n;
  
    d =  min(m,n);

    while (/*1*/ m % d != 0 || n % d != 0)
      { /*2*/
        d--;
      }

    /*3*/
    return d;
  }

  Aa descrição do algoritmo utilizou a linguagem C, mas nem sempre será assim. Muitas vezes optaremos por um descrição mais "alto nível" que enfatise os pontos realmente cruciais e não nos percamos com detalhes menos relevantes.

 

Correção

  Passamos agora a verificar a correção do algoritmo.
[Correção da função == a função funciona == a função faz o que promete.]

  A correção de algoritmos iterativos é comumente baseadas em demonstrações da validade de invariantes. Estes invariantes são afirmações ou relações envolvendo os objetos mantidos pelo algoritmo. Eis os nossos invariantes para o algoritmo. Em /*1*/ vale que

(i0) 1 <= d <= min(m,n), e
(i1) m % t != 0 ou n % t != 0 para todo número inteiro t > d,
e em /*2*/ vale que
(i2) em /*2*/ vale que m % d != 0 ou n % d != 0.
Verificaremos que essas afirmações, no momento candidatos a invariantes, são de fato válidas.
[A propósito, eu não sei se é a invariante ou o invariante.]

  Antes de mais nada, note que se de fato (i1) vale então o algoritmo devolve o prometido. É evidente que em /*3*/, antes do algoritmo devover d, vale que

m % d == 0 e n % d == 0.
Como (i1) vale em /*1*/, então (1) também vale em /*3*/. Donde, nenhum número inteiro maior que o número d devolvido divide m e n. Portanto, se o algoritmo devolve algum valor, este é de fato o mdc(m,n). O "se" é um certo excesso de zelo, já que, a rigor, ainda não verificamos que o algoritmo sempre pára quando alimentado com números inteiros satisfazendo as especificações.

  Invariantes são assim mesmo. A validade de alguns principais torna a correção do algoritmo (muitas vezes) evidente. Os invariantes secundários servem para confirmar a validade dos principais.
[Hmmmm, talvez os adjetivos principal e secundário sejam inapropriados.]

  Passemos a verificar a validade dos invariantes. A condição do while faz a validade (i2) em /*2*/ óbvia. Este foi fácil!

  Inicialmente d vale min(m,n) que é pelos menos 1, já que m e n são não nulos. Assim, quando a execução do programa alcança pela primeira vez o while tem-se que todo número inteiro maior que d == min(m,n) não divide m ou não divide n; em particular não divide min(m,n). Logo, (i0) e (i1) valem na primeira vez que a execução do algoritmo alcança o ponto /*1*/.

  Considere agora uma iteração do while em que o comando "d--;" é executado. Suponha que nesta iteração (i0) e (i1) valem em /*1*/, antes da execução do comando. Logo, (i0) e (i1) também valem em /*2*/, imediatamente antes da execução do comando d--;. Em /*2*/ também vale (i2). Portando, em /*2*/ valem (i0) e (i1) com d' == d-1 no papel de d. Como a iteração decrementa o valor de d, então (i0) e (i1) valem na primeira vez que a execução do algoritmo alcança o ponto /*1*/ novamente.

  Resumindo, acabamos de demonstrar que

se (i0) e (i1) valem em /*1*/ antes de uma execução do comando "d--", (i0) e (i1) valem na próxima vez que a execução alcaça o /*1*/
Como (i0) e (i1) valem na primeira vez que a execução passa por /*1*/, então (i0) e (i1) valem sempre que a execução do algoritmo está no ponto /*1*/. Portanto, (i0) e (i1) são invariantes genuinos.

  Invariantes, além de serem uma ferramente útil para demonstra a correção de algoritmos iterativo, eles nos ajudam a compreender o funcionamento do algoritmo. De certa forma, eles "espelham" a maneira que entendemos o algoritmo.

 

Desempenho

  Quantas iterações do while essa função faz? Equivalentemente, quantas vezes o comando "d--" é executado? A resposta é

min(m,n)-1   no pior caso.

Esse número de iterações ocorre quando mdc(m,n) = 1, ou seja, m e n são relativamente primos.

Por exemplo, para a chamada mdc(317811,514229) a função executará 317811 - 1 iterações e para a chamada mdc(2147483647,2147483646) a função faz 2147483646 - 1 iterações, pois mdc(317811,514229) == 1 e mdc(2147483647,2147483646) == 1.

Exercício

  1. Como diz um professor meu, só de farra, use o programa mdc.c para determinar mdc(2147483647,2147483646). Quanto tempo o programa gastou para fazer o serviço?

    Olha o tempo que demorou no computador que eu uso: 

        meu_prompt> time mdc 2147483647 2147483646 
    mdc de 2147483647 e 2147483646 e' 1.

    real    1m8.292s
    user    1m4.062s
    sys     0m1.101s

  Neste caso, costuma-se dizer que a quantidade de tempo consumida pelo algoritmo, no pior caso, é proporcional a min(m,n), ou ainda, que o tempo gasto pelo algoritmo é da ordem de min(m,n). A abreviação de "ordem de blá" é O(blá).
[Em breve, tornaremos o significado de O(blá) mais precisa.]

  Portanto, o número de operações feitas pela função cresce (no ior caso) como uma função linear de min(m,n). Isto significa que se o valor de min(m,n) dobra então espera-se que o tempo gasto pela função também dobre.

  É possível resolver o problema de uma maneira mais eficiente? A resposta é SIM.

 

 

Algoritmo de Euclides

  O máximo divisor comum entre dois número pode ser determinado através de um algoritmo de 2300 anos (cerca de 300 A.C.), o algoritmo de Euclides. Para calcular o mdc(m,n) para 0 <= n < m, o algoritmo de Euclides usa a seguinte recorrência:

mdc(m,0)   ==  m;
mdc(m,n)   ==  mdc(n, m % n),  para n > 0.

Assim, por exemplo,

mdc(12,18) == mdc(18,12) == mdc(12,6) == mdc(6,0) == 6.

 

Correção

A correção da recorrência proposta por Euclides é baseada no seguinte exercício.

Exercício

  1. Seja m um inteiro não-negativo e n um inteiro positivo. Para cada inteiro positivo d tem-se que

    (m % d == 0) e (n % d == 0)
      se e somente se  
    (n % d == 0) e ((m % n) % d == 0).

    Em outras palavras, o pares m,n e m,m%n têm o mesmo conjunto de divisores.

  Eis a implementação recursiva em C do algoritmo de Euclides para determinar mdc(m,n). 

 /*
  * Ao receber dois números inteiros não-negativos m e n, não
  * ambos nulos, está função devolve o máximo divisor de m e n.
  */  
  long
  euclides_r(long m, long n)
  {
    if (n==0) return m;
    return euclides_r (n, m % n);
  }

  A seqüência abaixo ilustra as chamadas do algoritmo de Euclides para determinar o mdc de 317811 e 514229: 

mdc(317811,514229)
  mdc(514229,317811)
    mdc(317811,196418)
      mdc(196418,121393)
        mdc(121393,75025)
          mdc(75025,46368)
            mdc(46368,28657)
              mdc(28657,17711)
                mdc(17711,10946)
                  mdc(10946,6765)
                    mdc(6765,4181)
                      mdc(4181,2584)
                        mdc(2584,1597)
                          mdc(1597,987)
                            mdc(987,610)
                              mdc(610,377)
                                mdc(377,233)
                                  mdc(233,144)
                                    mdc(144,89)
                                      mdc(89,55)
                                        mdc(55,34)
                                          mdc(34,21)
                                            mdc(21,13)
                                              mdc(13,8)
                                                mdc(8,5)
                                                  mdc(5,3)
                                                    mdc(3,2)
                                                      mdc(2,1)
                                                        mdc(1,0)
mdc de 317811 e 514229 e' 1.

Desempenho

  Quanto tempo o algoritmo leva para parar? Este tempo é certamente proporcional ao número de chamadas recursivas. Suponha que a função euclides_r faz k chamadas recursivas e que no início da 1a. chamada ao algoritmo tem-se que 0 < n <= m. Sejam

(m,n) == (m0,n0), (m1,n1),  (m2,n2),  ... ,  (mk,nk) == (mdc(m,n),0),

os valores dos parâmetros no início de cada uma das chamadas da função. Por exemplo, para m==514229 e n==317811 tem-se

(m0,n0) == (514229,317811),
(m1,n1) == (317811,196418),
(m2,n2) == (196418,121393),
        ... ,
(m27,n27) == (1,0).

Estimaremos o valor de k em função de n == min(m,n). Note que mi+1 == ni e ni+1 == mi % ni para i=1,2,...,k. Note ainda que para inteiros a e b, 0 < b <= a vale que

a % b   <   a/2   (verifique!).

Desta forma tem-se que

n2  ==  m1 % n1  ==  n0 % n1  <  n0/2  ==  n/2  ==  n/21
n4  ==  m3 % n3  ==  n2 % n3  <  n2/2  <  n/4   ==  n/22
n6  ==  m5 % n5  ==  n4 % n5  <  n4/2  <  n/8  ==  n/23
n8  ==  m7 % n7  ==  n6 % n7  <  n6/2  <  n/16  ==  n/24
n10  ==  m9 % n9  ==  n8 % n9  <  n8/2  <  n/32  ==  n/25
...
...
...

  Percebe-se que depois de cada 2 chamadas recursivas o valor do segundo parâmetro é reduzido a menos da sua metade. Seja t o número inteiro tal que

2t <= n < 2t+1

Da primeira desigualdade temos que

t <= lg(n),

onde lg(n) denota o logaritmo de n na base 2, ou seja lg(2) == log2(n). Da desigualde estrita, concluimos que

k <= 2t+1

Logo, o número k de chamadas recursivas é não superior a

2 (t+1) - 1 = 2t + 1 <= 2lg(n)+1.

[Quanto n > m tem-se que k <= 2lg(n)+2.]

  Para o exemplo acima, onde m==514229 e n==317811, temos que

2 lg(n) + 1 == 2 lg(317811) + 1 < 2 x 18.3 + 1 < 37.56

e o número de chamadas recursivas feitas por euclides_r(514229,317811) foram 27.

 

  Resumindo, a quantidade de tempo consumida pelo algoritmo de Euclides é, no pior caso, proporcional a lg(n). Este desempenho é significativamente melhor que o desempenho do algoritmo Café-com-leite, já que a função f(n) = lg(n) cresce muito mais lentamente que a função g(n) = n.

n   (int) lg(n)
4 2
5 2
6 2
10 3
64 6
100 6
128 7
1000 9
1024 10
1000000 19
1000000000 29

Exercícios

  1. Demonstre que para inteiros m, k e n mdc(m,n) == mdc(m + kn, n).

    2.  

  2. É verdade que para todo m >= n > 0 a função mdc faz um número de iterações (do comando while) maior ou igual ao feito (pelo comando do ... while) da função euclides_i abaixo. 
       long
   euclides_i(long m, long n)
   {
     long r;

     do
       {
         r = m % n;
         m = n;
         n = r;
       }
     while (r != 0);

     return m;
   }

    3.  

  3. Só de farra, compare, na prática, o consumo de tempo das funções mdc.c, euclides_r e euclides_i.

     

  4. Qual função consome uma quantidade maior de memória, euclides_r ou euclides_i? Faça uma estimativa da quantidade de memória consumida por cada uma dessas funções.

 

 

Números de Fibonacci e o algoritmo de Euclides

  Os 10 primeiros números da seqüência de Fibonacci são:

0 1 2 3 4 5 6 7 8 9 10
0 1 1 2 3 5 8 13 21 34 55

  Os números de Fibonacci são definidos (ou podem ser calculados) pela seguinte função recursiva: 

 /*
  * A função receber um número inteiros não-negativos  n, e
  * devolve o n-ésimo número de Fibonacci.
  */  
   long
   fibonacci(long n)
   {
     if (n == 0) return 0;
     if (n == 1) return 1;
     return  fibonacci(n-1) + fibonacci(n-2);
   }

Os números de Fibonacci estão intimamente relacionados com o algoritmo de Euclides, como mostra o exercício abaixo.

Exercício

  1. Se m > n >= 0 e se a chamada euclides_r(m,n) faz k >=1 chamadas recursivas, então m >= fibonacci(k+2) e n >= fibonacci(k+1).
    [Sugestao: demonstre por indução em k.]

  Uma conseqüência imediata deste exercício é

Para todo inteiro t >= 1, se m > n >= 0 e n < fibonacci(t+1), então a chamada euclides_r(m,n) faz menos de t chamadas recursivas.

Por exemplo, fibonacci(28) == 317811 e fibonacci(29) = 514229 e a chamada euclides_r(514229,317811) faz 27 chamadas recursivas. Na verdade, números de Fibonacci consecutivos são os que dão mais trabalho para o algoritmo de Euclides (veja exercício 8)..

  Vejamos agora como utilizando os fatos acima e algo que um passarinho me contou podemos obter uma estimativa melhor para o número de chamadas recursivas de euclides_r(m,n).

  Suponha que m > n >= 4 e que um passarinho me contou que para todo número inteiro não-negativo vale que

phit-2 <= fibonacci(t) <= phit-1,

para t >= 2, onde phi == (1+sqrt(5))/2 que é um número entre 1.618 e 1.619. Por exemplo,

phi26 < 1.61926 < 275689 < fibonacci(28)==317811 < 438954 < 1.61827 < phi27

[Hmmm, parece que o passarinho não mentiu.]

  Seja t o número inteiro tal que [(neste ponto precisamos da hipótese n >= 4]

phit+2 <= n < phit+3 < fibonacci(t+1)

então o número k de chamadas recursivas de euclides_r(m,n) é no máximo

t - 1 < (logphi(n)+2) - 1 == logphi(n) + 1.

  Segundo está nova estimativa temos que euclides_r(514229,317811) faz no máximo

logphi(317811) + 1 > log1.619(317811) + 1 > 26.28 + 1 = 27.45.

e o número de chamadas é 27.
[Uauuuu. Isto foi muito perto! Talvez tenha alguma conta errada.]

Exercícios

  1. Demonstre que a chamada euclides_r(fibonacci(k+1),fibonacci(k)) faz k-1 chamadas recursivas.
    [Sugestao: demonstre por indução em k.]

    2.  

  2. Demonstre que números de Fibonacci consecutivos são relativamente primos.

    3.  

  3. Eis a versão iterativa da função fibonacci. 
      /*
  * A função receber um número inteiros não-negativos  n, e
  * devolve o n-ésimo número de Fibonacci.
  */      
  long
  fibonacci(long n)
  {
    long i;
    long fib_ant, fib_atual, fib_prox;

    if (n == 0) return 0;

    fib_ant = 0; fib_atual = 1; 
    for (i = 1; i < n; i++)
      { 
        fib_prox  = fib_ataul + fib_ant;
        fib_ant   = fib_atual;
        fib_atual = fib_prox;
      }

    return  fib_atual;
  }
    Qual função você usaria? A recursiva ou a iterativa? Por quê? A propósito, note que, no início da i-ésima iteração vale que
    (i0) fib_ant é o (i-1)-ésimo número de fibonacci,
    (i1) fib_atual é o i-ésimo número de finonacci.
    Quando a função pára tem-se que i == n e portanto, de (i1), conclui-se que a função devolve o n-ésimo número de Fibonacci.

     

  4. Matematicamente falando, quanto é 7 % -3? Qual é o valor que o computador devolve? No computador vale que m == (m/n) * n + (m%n)?

     

  5. Demonstre que 1 + phi == phi2, onde phi == (1+sqrt(5))/2.

     

  6. Demonstre a afirmação do passarinho de que

    phit-2 <= fibonacci(t) <= phit-1,

    para t >= 2.

 

 

Algoritmo de Euclides estendido

  O algoritmode Euclides pode ser facilmente modificado para nos fornecer mais que o mdc(m,n). Podemos estender o algoritmo para que ele compute números inteiros a e b tais que

am + bn == mdc(m,n).
Se n == 0, então podemos tomar a == 1 e b == 0. Caso contrário, aplique o método recursivamente com n e m%n nos papéis de m e n, respectivamente, computando assim inteiros e e f tais que
en + f(m%n) == mdc(n,m%n).
Como m%n == m-(m/n)m  e  mdc(n,m%n) == mdc(m,n), esta equação nos diz que
en + f(m-(m/n)m) == mdc(m,n).
O lado esquerdo da equação pode ser reescrito para mostrar a dependência de m e n
fm + (e-f(m/n))n == mdc(m,n).
Portanto a == f  e  b = (e - f(m/n)) são os inteiros desejados. Por exemplo, para 
  mdc(3267,2893)
    mdc(2893,374)
      mdc(374,275)
        mdc(275,99)
          mdc(99,77)
            mdc(77,22)
              mdc(22,11)
                mdc(11,0)
  mdc de 3267 e 2893 e' 11
e temos que 
  11 ==    1 x 11   +    0 x 0   ==   0 x 22  +  1 x 11
     ==    1 x 77   +   -3 x 22  ==  -3 x 99  +  4 x 77
     ==    4 x 275  +  -11 x 99  == -11 x 374 + 15 x 275
     =  = 15 x 2893 + -116 x 374
     == -116 x 3267 +  131 x 2893.

  Eis uma versão em C da descrição acima 

 /*
  * Ao receber dois números inteiros não-negativos m e n, não
  * ambos nulos, está função devolve o máximo divisor de m e n
  * e nos parâmetros *a e *b devolve números inteiros
  * tais que
  *           (*a)*m + (*b)*n == mdc(m,n).
  */  
  long
  euclides_r_x(long m, long n, long *a, long *b)
  {
     long d;
    
     if (n==0) 
       { 
         *a = 1;
         *b = 0; 
         return m;
       }
     else 
       {
         long e, f;

         d = euclides_r_x (n, m % n, &e, &f);

        /* Neste ponto vale que
              d == e*n + f*(m%n)
                == e*n + f*(m - (m/n)*n)
                == f*m + (e - f*(m/n))*n   */

         *a = f;
         *b = e - f*(m/n);

        return d;
      }
  
    }

  Mas qual a razão de tanta confusão? A principal razão (pelo menos neste ponto da disciplina) é que de uma certa maneira os números a e b são uma prova ou certificado da correção da resposta. Suponha que nosso programa afirme que o máximo divisor comum de dois números gigantescos m e n é d e alám disso nos fornece números inteiros a e b tais que

am + bn == d.
Suponha ainda ques estamos meio deconfiados da correção da resposta. Será que pode existir um inteiro positivo d' maior que d e que divida m e n? A resposta é não! Se existe d' < d que divida m e n, então d' divide am + bn == d, e portanto d' <= d.

  Para verificar que a resposta do nosso programa está correta, basta verificarmos que

am + bn == d,
e que d divide m e n. Esta verificação custa 2 multiplicações, 2 divisões, 1 adição e 3 comparações.

  Algoritmos que devolvem certificados da sua correção sao chamados self-certifying ou simpesmente certifying. Na página pessoal de Kurt Mehlhorn há um link para uma palestra sobre Certifying Algorithms.

 

 

 

Last modified: Thu Mar 11 10:47:28 BRT 2004