Dia desses o Lameiro me passou uma lista de livros para desenvolvedores, e eu achei a lista bastante curiosa. Se eu fosse criar uma similar, eu não colocaria nenhum dos livros citados naquela lista! Ao invés disso, a minha lista com os top 10 livros de computação seria a abaixo:
1 - Gödel, Escher, Bach (Hofstadter): Eu começaria com o GEB, por dois motivos. O primeiro é que ele é um ótimo reality check: se você não gostar do GEB, então mude de área, porque computação não é a sua praia :) O segundo motivo é que esse livro tem uma excelente introdução à lógica (proposicional e de predicados), que é a ferramenta básica onde você constrói a ciência da computação.
2 - Concrete Mathematics (Knuth): Você não precisa saber matemática para programar, mas se você quiser ser um bom programador, então matemática é essencial. O Concrete tem todo o básico que você precisa pra fazer análises de complexidade computacional, e tudo escrito de maneira extremamente bem-humorada.
3 - Algorithms in C++ (Sedgewick): Se você já sabe lógica e matemática, então agora pode partir pro estudo de algoritmos. O Sedgewick tem todos os algoritmos básicos, e é uma leitura bem leve: se você está começando com algoritmos agora, esse precisa ser o seu primeiro livro. Ele não vai muito fundo em nenhum tópico, mas isso é compensado pela extrema didática nos tópicos. O livro ainda tem código exemplo pra todos os algoritmos, e várias edições, uma pra cada linguagem (eu sei que tem, pelo menos, C, C++, Java e Pascal).
4 - Introduction to Algorithms (Cormen): Os algoritmos que você aprendeu no Sedgewick, você vai estudar em detalhes no Cormen. Esse livro é extremamente formal, e talvez por isso é o livro-texto usado nos cursos de computaçao do MIT. Ele também cobre algoritmos mais avançados, que o Sedgewick apenas cita (por exemplo, Fibonacci Heap)
5 - The Art of Computer Programming (Knuth): O tAoCP está para o Cormen assim como o Cormen está pro Segdewick, aqui você vai dissecar os algoritmos até o último bit deles. Além de ser uma coleção excelente, a encadernação é muito bonita (mas não se engane, ele fica ótimo na prateleira, mas melhor ainda na sua cabeça).
6 - Effective C++ (Meyers): Agora que você sabe os algoritmos, precisa de uma linguagem para programá-los. Pra falar a verdade, a escolha de linguagem nem importa tanto assim, mas se você escolher uma, aprenda-a tão bem quanto possível. Eu escolhi C++, e esse livro do Meyers é o que diferencia as crianças dos adultos (especialmente para aquelas horas quando você cria uma classe sem destrutor virtual e não sabe por que a memória está vazando).
7 - Effective STL (Meyers): Já teve uma época em que eu não gostava de C++, mas isso é porque eu sou velho o suficiente pra ter mexido em C++ antes que os compiladores tivessem templates. Com templates a linguagem fica muito mais atraente, e esse é o livro que vai te ensinar a dominar a STL.
8 - The Practice of Programming (Kernighan, Pike): Se você leu tudo até agora, então você já é um programador muito bom na teoria. Na prática, entretanto, tem um monte de skills que ainda faltam. Nesse livro você aprende sobre as coisas que usualmente não se aprende na escola: debugging, otimização, unit testing, documentação.
9 - Programming Pearls (Bentley): Com os livros lidos até agora, você já deve ser um excelente programador. O passo final é passar de programador para um true hacker, e esse é um passo que não requer só conhecimento, você também precisa de manha, criatividade e insight. Eu não sei se dá pra ensinar essas coisas, mas esse livro certamente é o que chega mais próximo disso.
10 - The Mythical Man-Month (Brooks): Depois de ler todos os livros acima você estará próximo do nirvana, mas pra atingir o zen da programação de verdade, é preciso lembrar que projetos não precisam só de computadores, precisam de pessoas também. O Mythical Man-Month é um livro de gerência de projetos de software escrito em 1975, mas é surpreendente como ele continua atual. A tecnologia avança, mas as pessoas continuam as mesmas :)
É claro que pra manter uma lista com só dez itens, muita coisa boa fica de fora. Mas a lista acima tem um mérito: foi com esses livros que eu aprendi computação de verdade (vale lembrar que eu sou autodidata, minha graduação foi em engenharia elétrica, e eu quase não tive computação em aulas). Se funcionou pra mim, pode ser que funcione pra você também :)
segunda-feira, 23 de junho de 2008
quarta-feira, 18 de junho de 2008
Um cientista em minha vida
Eu li lá no blog do Kentaro que o meme da semana era "Um Cientista em minha vida", onde deveríamos falar sobre algum cientista que fez a diferença pra você. Como eu adoro constrained writing, resolvi participar (na verdade, eu adoro constrained anything, por isso que vivo criando programas em uma linha, programas que rodam em computadores de 8 bits, e assim por diante).
Eu já falo de cientistas aqui todo o tempo. Olhando no histórico, eu já falei do Knuth, do Erdös, do prof. Routo, do prof. Henrique, e de vários outros. Em comum, todos eles foram cientistas que eu conheci depois de adulto. Achei apropriado então que eu falasse de um cientista que fez a diferença quando eu era criança, e pra isso vamos ter que rebobinar até a década de 80.
Se você perguntar pra alguém sobre revistas de computador na década de 80, invariavelmente irá ouvir sobre a Micro Sistemas ("a primeira revista brasileira sobre microcomputadores"). A Micro Sistemas era muito legal, mas o que eu gostava mesmo era de outra revista, menos conhecida, chamada Microhobby.
A diferença da Micro Sistemas pra Microhobby era mais ou menos a diferença de Informática pra Computação. Na primeira, nós ficávamos encantados com as notícias da maravilhosa terra além da reserva de mercado (onde aprendíamos que a Apple planejava lançar um novo computador chamado McIntosh, que vinha com um periférico estranho e esquisito chamado mouse), enquanto que na segunda aprendíamos a calcular geodésicas e a usar o método de Bolzano para achar raízes de uma equação.
Mas o diferencial mesmo da Microhobby eram as colunas escritas pelo Renato da Silva Oliveira. Uma googlada rápida revela que o Renato é formado em Física, trabalhou nos planetários de São Paulo, Campinas, Vitória e Tatuí, e atualmente trabalha em uma empresa que vende planetários infláveis (how cool is that?!). Mas é claro que eu não sabia disso na época, o que eu sabia era que ele contava historinhas!
Foi lendo as historinhas do Renato que eu descobri que era possível escrever sobre ciência e computação, com clareza e bom humor. Pena que isso ainda não é muito difundido, a julgar pela quantidade de crianças que ainda acham que ciência é uma coisa chata :(
As historinhas que ele escrevia sempre tinham o mesmo formato: um certo sr. Nabor Rosenthal, em suas viagens pelo mundo, deparava-se com alguma situação que sugeria uma análise matemática (os tópicos eram os mais variados e iam de teoria dos grafos até contatos com extraterrestres). Depois de ponderar sobre o problema sem conseguir resolvê-lo, o Nabor tomava uma dose do raro Suco de Ramanujan, que o colocava num transe que ampliava suas capacidades analíticas, e conseguia solucionar o problema.
Mas a coluna sempre acabava antes que o Nabor mostrasse qual a solução! Ao invés disso, o leitor tinha um mês pra conseguir resolver o problema, e só no mês seguinte a solução era apresentada. Na década de 80 ainda não tinha spoj, então as colunas do Renato eram o que bombava pra quem gostava de puzzles computacionais.
Um dos puzzles apresentados foi o segundo puzzle mais difícil da minha vida, eu levei mais de dez anos pra conseguir resolver. Em uma das Microhobby, o Nabor entrou em transe após tomar o Suco de Ramanujan, e durante o transe ele sonhou "com um método para calcular o número pi, usando apenas o gerador de números aleatórios de seu micro" (essa era a época onde o micro mais avançado era o TK-82C, com 2kb de RAM).
Na época eu pensei muito e não consegui solucionar, achando que ia precisar de alguma matemática que eu ainda não tinha aprendido. Eu nunca consegui achar a revista seguinte com a solução, tive que passar pelo primário, pelo colégio técnico, e só no meio da faculdade é que caiu a ficha (e eu percebi que poderia ter solucionado ainda no primário, se tivesse insistido o suficiente :)
O truque é o seguinte: você vai fazer N experimentos, cada um consistindo no sorteio de dois números aleatórios escolhidos uniformemente entre 0 e 1. Se soma dos quadrados dos números for menor ou igual a 1, incremente um contador (digamos, M). Ao final dos experimentos, pi=4*M/N. O script abaixo implementa esse algoritmo:
Script em python para calcular pi usando números aleatórios
O funcionamento é bem simples e baseia-se na figura ao lado. Você começa inscrevendo um quarto de círculo num quadrado de lado 1. Os dois números que você sorteia a cada iteração podem ser interpretados como um ponto dentro do quadrado, e o teste feito é equivalente a testar se o ponto está dentro do círculo ou não. Como a distribuição dos pontos é uniforme, espera-se que a razão M/N seja igual à razão entre as áreas da figura. A área do quadrado é 1, a área do círculo é pi*r2. Como o raio é unitário, então a área do quarto de círculo é pi/4. Isolando pi, chega-se em pi=4*M/N, QED.
A pergunta que deve ser feita ao encontrar qualquer algoritmo novo é: qual é sua complexidade? Infelizmente, esse método aleatório é bem ruim. No fundo, o que estamos fazendo é aproximar pi por uma fração, cujo denominador é N. Então a precisão máxima que podemos obter é 1/N, e se você quer calcular n dígitos de pi, esse método converge, no melhor caso, em O(10n), e na prática em bem menos que isso, porque os seus geradores de números aleatórios não são perfeitamente uniformes.
Eu nunca soube qual o método que o Nabor usou pra calcular o pi. Como ele tinha o Suco de Ramanujan e eu não, espero que tenha sido um método melhor que o meu :)
Update!
Em março de 2011 comemoram-se os 30 anos do ZX81, então decidi revisar esse post com o conteúdo original da época. Esse aqui era o enunciado original da Microhobby, clique para ampliar:
Eu implementei a solução usando o BASIC de um TK85 (real, não é emulador):
Se não der para ler na foto, o programa em BASIC é o abaixo:
1 FAST 4 LET I=0 5 LET M=5000 6 LET T=0 10 FOR N=0 TO M 20 LET X=RND 30 LET Y=RND 40 LET T=T+1 50 IF X*X+Y*Y<1 THEN LET I=I+1 60 NEXT N 70 PRINT 4*I/TO resultado, após 5000 iterações, foi 3.1137773. Conseguimos duas casas de precisão, após dez minutos de processamento de uma tecnologia de 30 anos atrás :)
Outro update!
O Leonardo Roman da Rosa achou a solução original da revista, saiu na Microhobby #25. No final, a solução original é exatamente a mesma que descobri! Até as constantes foram coincidemente iguais :)
quarta-feira, 11 de junho de 2008
Programa de milhagem
Olhando no Google Analytics, eu descobri que alguém chegou aqui no blog procurando por "transformar kilometros em arquivo binarios". Se você é essa pessoa, desculpe, mas eu não entendi sua dúvida. Se você não é essa pessoa, puxe uma cadeira pra ver como até uma pergunta sem sentido pode ser desenvolvida num tema interessante :)
Uma coisa que me incomoda toda vez que venho pra Califórnia é ter que lidar com milhas e libras. Eu cresci com o sistema métrico: se você me disser que a distância de São Paulo pra Florianópolis é de 700km, eu sei o que isso significa. Agora, se você me disser que a distância de San Francisco pra Mountain View é de 100 milhas, minha intuição falha, e eu vou precisar converter pra km pra poder ter noção da distância.
Como esperado, converter mentalmente de milhas pra km é algo que faço todo o tempo por aqui. Existem várias maneiras de converter de cabeça, mas a filosofia Ricbit dita que, de várias maneiras equivalentes, o correto é escolher a mais bizarra! Sendo assim, vou mostrar a conversão utilizando números de Fibonacci.
Todo mundo conhece os números de Fibonacci, eles estão em todo lugar. Em minha época de estudante, uma das minhas diversões secretas era entrar escondido no andar superior da biblioteca do IME só pra ler a Fibonacci Quartely. Os leitores assíduos da revista conhecem um monte de fatos curiosos sobre os Fibonacci, e eu vou usar dois deles aqui.
O primeiro é que os números da seqüência de Fibonacci podem ser calculados diretamente, sem precisar fazer toda a recursão F(n+2)=F(n+1)+F(n). Pra isso, basta calcular qual o inteiro mais próximo da expressão abaixo:
Na expressão acima, o phi é a conhecida razão áurea. A demostração dessa fórmula é elementar e fácil de encontrar na web.
O segundo é que qualquer número pode ser escrito como uma soma de números distintos de Fibonacci. Por exemplo, 100 pode escrito como 89+8+3. Daí, se você enfileirar os números de Fibonacci, e atribuir a cada número 1 se ele é usado na soma, e 0 se não é, você pode atribuir a qualquer inteiro uma string de zeros e uns que funciona como uma espécie de base binária alternativa (o povo chama isso de base de Fibonacci).
Fazendo o processo com o número 100, chegamos em 10000101000. Essa base tem algumas propriedades curiosas, por exemplo, ela não é bijetora (de fato, você pode escrever 100 de outras maneiras, como 1100101000). Uma base numérica que tem representação múltipla tem utilidades bastante curiosas em design de circuitos elétricos (mas isso é uma história pra outro dia :).
Além disso, cada número possui pelo menos uma representação onde não há nenhuma seqüência com dois uns consecutivos (pela própria definição de Fibonacci, se houver dois uns em algum ponto, você pode apagá-los e trocar por um único 1 na posição seguinte). Esse fato é explorado em alguns tipos de sinalização, para fazer detecção de erro: se você receber dois uns seguidos, certamente recebeu um erro.
Mas qual a relação disso com milhas e quilômetros? É simples, para converter de milhas para km, basta fazer um shift de Fibonacci!
Uma milha equivale a 1.609344 km. A razão áurea é 1.61803399. Os dois números, apesar de não serem relacionados, são muito parecidos, e esse é o truque que vamos usar pra converter. Na base binária tradicional, um shift para a esquerda é equivalente a multiplicar por dois; na base de Fibonacci, um shift para a esquerda equivale a multiplicar pela razão áurea. Então, se você tiver um valor em milhas na base de Fibonacci, um shift irá transformar o valor para o equivalente em quilômetros.
Vamos conferir: 100 é 10000101000. Com um zero extra no final, fica 100001010000, que é 144+13+5=162. Se, ao invés disso, você converter diretamente, teria 160.9km, ou seja, o método realmente aproxima muito bem a conversão! O gráfico abaixo mostra a porcentagem do erro do método em relação ao ideal, e mais abaixo está o programa que converte a milhagem para quilometragem e gera o gráfico:
Script que plota o gráfico acima, em python
Como pode ser visto, o método esquenta bem rápido, pra distâncias superiores a 10 milhas o erro é inferior a 5%, e acima disso praticamente desaparece. Em comparação, o método naive de aproximar por 1.5 (multiplicar de cabeça por 3 e dividir por 2), tem um erro constante de mais ou menos 8%.
Uma coisa que me incomoda toda vez que venho pra Califórnia é ter que lidar com milhas e libras. Eu cresci com o sistema métrico: se você me disser que a distância de São Paulo pra Florianópolis é de 700km, eu sei o que isso significa. Agora, se você me disser que a distância de San Francisco pra Mountain View é de 100 milhas, minha intuição falha, e eu vou precisar converter pra km pra poder ter noção da distância.
Como esperado, converter mentalmente de milhas pra km é algo que faço todo o tempo por aqui. Existem várias maneiras de converter de cabeça, mas a filosofia Ricbit dita que, de várias maneiras equivalentes, o correto é escolher a mais bizarra! Sendo assim, vou mostrar a conversão utilizando números de Fibonacci.
Todo mundo conhece os números de Fibonacci, eles estão em todo lugar. Em minha época de estudante, uma das minhas diversões secretas era entrar escondido no andar superior da biblioteca do IME só pra ler a Fibonacci Quartely. Os leitores assíduos da revista conhecem um monte de fatos curiosos sobre os Fibonacci, e eu vou usar dois deles aqui.
O primeiro é que os números da seqüência de Fibonacci podem ser calculados diretamente, sem precisar fazer toda a recursão F(n+2)=F(n+1)+F(n). Pra isso, basta calcular qual o inteiro mais próximo da expressão abaixo:
Na expressão acima, o phi é a conhecida razão áurea. A demostração dessa fórmula é elementar e fácil de encontrar na web.
O segundo é que qualquer número pode ser escrito como uma soma de números distintos de Fibonacci. Por exemplo, 100 pode escrito como 89+8+3. Daí, se você enfileirar os números de Fibonacci, e atribuir a cada número 1 se ele é usado na soma, e 0 se não é, você pode atribuir a qualquer inteiro uma string de zeros e uns que funciona como uma espécie de base binária alternativa (o povo chama isso de base de Fibonacci).
Fazendo o processo com o número 100, chegamos em 10000101000. Essa base tem algumas propriedades curiosas, por exemplo, ela não é bijetora (de fato, você pode escrever 100 de outras maneiras, como 1100101000). Uma base numérica que tem representação múltipla tem utilidades bastante curiosas em design de circuitos elétricos (mas isso é uma história pra outro dia :).
Além disso, cada número possui pelo menos uma representação onde não há nenhuma seqüência com dois uns consecutivos (pela própria definição de Fibonacci, se houver dois uns em algum ponto, você pode apagá-los e trocar por um único 1 na posição seguinte). Esse fato é explorado em alguns tipos de sinalização, para fazer detecção de erro: se você receber dois uns seguidos, certamente recebeu um erro.
Mas qual a relação disso com milhas e quilômetros? É simples, para converter de milhas para km, basta fazer um shift de Fibonacci!
Uma milha equivale a 1.609344 km. A razão áurea é 1.61803399. Os dois números, apesar de não serem relacionados, são muito parecidos, e esse é o truque que vamos usar pra converter. Na base binária tradicional, um shift para a esquerda é equivalente a multiplicar por dois; na base de Fibonacci, um shift para a esquerda equivale a multiplicar pela razão áurea. Então, se você tiver um valor em milhas na base de Fibonacci, um shift irá transformar o valor para o equivalente em quilômetros.
Vamos conferir: 100 é 10000101000. Com um zero extra no final, fica 100001010000, que é 144+13+5=162. Se, ao invés disso, você converter diretamente, teria 160.9km, ou seja, o método realmente aproxima muito bem a conversão! O gráfico abaixo mostra a porcentagem do erro do método em relação ao ideal, e mais abaixo está o programa que converte a milhagem para quilometragem e gera o gráfico:
Script que plota o gráfico acima, em python
Como pode ser visto, o método esquenta bem rápido, pra distâncias superiores a 10 milhas o erro é inferior a 5%, e acima disso praticamente desaparece. Em comparação, o método naive de aproximar por 1.5 (multiplicar de cabeça por 3 e dividir por 2), tem um erro constante de mais ou menos 8%.
sábado, 7 de junho de 2008
Primos aleatórios
Dia desses a Alice me perguntou se era possível criar um gerador de números aleatórios que só retornasse números primos. Eu respondi que sim, mas que provavelmente ela não iria gostar da resposta:
Eu sabia que o que ela queria na verdade era uma fórmula bonitinha; então, como esperado, ela não gostou :) Mas a verdade é que esse algoritmo é bem melhor que as alternativas!
Antes de mostrar porque isso é verdade, precisamos formalizar um pouco o problema. É claro que não existem algoritmos que geram números aleatórios: se você quiser aleatoriedade real, precisa pegar alguma fonte física, como o decaimento radiativo. Assumindo então que existe uma fonte física que gera uma distribuiçao uniforme sobre algum intervalo, para criar o algoritmo que retorna números primos aleatórios, basta criar uma função bijetora que leve naturais para primos. Ou seja, uma função que, para um dado um número n, retorne o n-ésimo primo.
O problema é que não existe nenhuma fórmula fechada que calcule isso de maneira eficiente. Você pode calcular alguma constante irracional que resolva o problema, no estilo da constante de Mills, só que mais cedo ou mais tarde a precisão vai te limitar. Você pode calcular o n-ésimo primo com base em alguma outra distribuição, como a função de Möbius, mas aí você só está empurrando o problema com a barriga, porque a outra função é tão difícil de calcular quanto a original.
Uma maneira sem as desvantagens acima é usar o teorema de Wilson pra chegar na seguinte fórmula:
Mas mesmo essa fórmula ainda está longe do ideal, primeiro porque você vai ter que lidar com números enormes nela (pra n=10 os valores intermediários ficam tão grandes que estouram o limite do que cabe num float), segundo porque, mesmo que você use uma lib para long floats, a complexidade é O(2n), ou seja, mais lento que os programadores do Duke Nukem Forever. Se ainda assim você quiser testar, minha implementação em python é a abaixo:
Implementação em python da fórmula acima
Sendo assim, quão melhor era a implementação original por tentativa e erro? Pra avaliar isso, precisamos calcular a complexidade daquele algoritmo. Não é difícil ver que a complexidade do algoritmo como um todo é a complexidade do is_prime() multiplicado pelo valor esperado do número de iterações do loop.
Se você estiver trabalhando numa faixa pequena de primos, pode tabelar todos os primos no intervalo e fazer um is_prime() que seja O(1), mas aí também não tem necessidade da tentativa e erro, você pode indexar seu número aleatório direto na tabela. O caso legal é quando você não pode tabelar, nesse caso você pode implementar o is_prime() usando, por exemplo, o algoritmo AKS, cuja complexidade é O((log n)10.5).
O que resta então é calcular o valor esperado do loop. Lembrando que E[x]=sum(x*p(x)), o que precisamos é calcular qual é a probabilidade de ter uma iteração, duas iterações, e assim por diante. Ora, o teorema dos números primos nos garante que a quantidade de números primos menores que n é assintoticamente igual a n/log(n), então a chance de um número ser primo, num conjunto com n elementos, é 1/log(n). Vamos chamar isso de "p" só pra ficar mais fácil, e o complemento disso é q=1-p, ou seja, a chance de um número não ser primo.
Vejamos então: pra você acertar o primo de primeira, a chance é p. Se você acertar o primo na segunda, a chance é pq. Na terceira, é pq2, na quarta pq3 e assim por diante. Então o valor esperado é:
X = 1p + 2pq + 3pq2 + 4pq3 + ...
X = p (1 + 2q + 3q2 + 4q3 + ....)
Quem tem prática com a transformada z sabe calcular isso de cabeça, mas dá pra calcular também só com matemática elementar. Se você isolar q na soma, fica com:
X = p (1 + q(2 + 3q + 4q2 + ....))
Agora você tira da cartola y=1+q+q2+q3+... e substitui:
X = p (1 + q(2 + 3q + 4q2 + ....))
X = p (1 + q(y + 1 + 2q + 3q2 + ....))
X = p (1 + q(y + X/p)) = p + pqy + pXq/p = p(1+qy) + Xq
X - Xq = p (1 + qy)
X (1-q) = Xp = p (1 + qy)
X = 1 + qy
Mas y é só a soma de uma PG, e isso nós sabemos que vale y=1/(1-q)=1/p. Então:
X = 1 + q/p = (p+q)/p = 1/p
Como p=1/log(n), então o valor esperado que nós queríamos é tão somente X=log n (vocês também não se impressionam quando tudo simplifica no final?)
É claro que eu não iria resistir à tentação de implementar uma simulação pra ver se o valor bate mesmo. A nossa fórmula diz que, para a faixa de 10 milhões de números, o valor esperado tem que ser da ordem de log(107)=16.1. A simulação abaixo retorna 15.2, bem próximo do valor que foi predito.
Simulação monte carlo do valor esperado, em C++
No fim das contas, a complexidade do algoritmo com tentativa e erro é apenas O(log n), se você tiver um tabelão de primos. Na prática, esse é o método usado por todos que precisam de primos aleatórios: a libgcrypt usada no gpg, por exemplo, utiliza esse método na função gen_prime(), com vários truques pra tornar o teste de primalidade bem rápido.
int random_prime(int n) {
int x;
do {
x = random(n);
} while (!is_prime(x));
return x;
}
Eu sabia que o que ela queria na verdade era uma fórmula bonitinha; então, como esperado, ela não gostou :) Mas a verdade é que esse algoritmo é bem melhor que as alternativas!
Antes de mostrar porque isso é verdade, precisamos formalizar um pouco o problema. É claro que não existem algoritmos que geram números aleatórios: se você quiser aleatoriedade real, precisa pegar alguma fonte física, como o decaimento radiativo. Assumindo então que existe uma fonte física que gera uma distribuiçao uniforme sobre algum intervalo, para criar o algoritmo que retorna números primos aleatórios, basta criar uma função bijetora que leve naturais para primos. Ou seja, uma função que, para um dado um número n, retorne o n-ésimo primo.
O problema é que não existe nenhuma fórmula fechada que calcule isso de maneira eficiente. Você pode calcular alguma constante irracional que resolva o problema, no estilo da constante de Mills, só que mais cedo ou mais tarde a precisão vai te limitar. Você pode calcular o n-ésimo primo com base em alguma outra distribuição, como a função de Möbius, mas aí você só está empurrando o problema com a barriga, porque a outra função é tão difícil de calcular quanto a original.
Uma maneira sem as desvantagens acima é usar o teorema de Wilson pra chegar na seguinte fórmula:
Mas mesmo essa fórmula ainda está longe do ideal, primeiro porque você vai ter que lidar com números enormes nela (pra n=10 os valores intermediários ficam tão grandes que estouram o limite do que cabe num float), segundo porque, mesmo que você use uma lib para long floats, a complexidade é O(2n), ou seja, mais lento que os programadores do Duke Nukem Forever. Se ainda assim você quiser testar, minha implementação em python é a abaixo:
Implementação em python da fórmula acima
Sendo assim, quão melhor era a implementação original por tentativa e erro? Pra avaliar isso, precisamos calcular a complexidade daquele algoritmo. Não é difícil ver que a complexidade do algoritmo como um todo é a complexidade do is_prime() multiplicado pelo valor esperado do número de iterações do loop.
Se você estiver trabalhando numa faixa pequena de primos, pode tabelar todos os primos no intervalo e fazer um is_prime() que seja O(1), mas aí também não tem necessidade da tentativa e erro, você pode indexar seu número aleatório direto na tabela. O caso legal é quando você não pode tabelar, nesse caso você pode implementar o is_prime() usando, por exemplo, o algoritmo AKS, cuja complexidade é O((log n)10.5).
O que resta então é calcular o valor esperado do loop. Lembrando que E[x]=sum(x*p(x)), o que precisamos é calcular qual é a probabilidade de ter uma iteração, duas iterações, e assim por diante. Ora, o teorema dos números primos nos garante que a quantidade de números primos menores que n é assintoticamente igual a n/log(n), então a chance de um número ser primo, num conjunto com n elementos, é 1/log(n). Vamos chamar isso de "p" só pra ficar mais fácil, e o complemento disso é q=1-p, ou seja, a chance de um número não ser primo.
Vejamos então: pra você acertar o primo de primeira, a chance é p. Se você acertar o primo na segunda, a chance é pq. Na terceira, é pq2, na quarta pq3 e assim por diante. Então o valor esperado é:
X = 1p + 2pq + 3pq2 + 4pq3 + ...
X = p (1 + 2q + 3q2 + 4q3 + ....)
Quem tem prática com a transformada z sabe calcular isso de cabeça, mas dá pra calcular também só com matemática elementar. Se você isolar q na soma, fica com:
X = p (1 + q(2 + 3q + 4q2 + ....))
Agora você tira da cartola y=1+q+q2+q3+... e substitui:
X = p (1 + q(2 + 3q + 4q2 + ....))
X = p (1 + q(y + 1 + 2q + 3q2 + ....))
X = p (1 + q(y + X/p)) = p + pqy + pXq/p = p(1+qy) + Xq
X - Xq = p (1 + qy)
X (1-q) = Xp = p (1 + qy)
X = 1 + qy
Mas y é só a soma de uma PG, e isso nós sabemos que vale y=1/(1-q)=1/p. Então:
X = 1 + q/p = (p+q)/p = 1/p
Como p=1/log(n), então o valor esperado que nós queríamos é tão somente X=log n (vocês também não se impressionam quando tudo simplifica no final?)
É claro que eu não iria resistir à tentação de implementar uma simulação pra ver se o valor bate mesmo. A nossa fórmula diz que, para a faixa de 10 milhões de números, o valor esperado tem que ser da ordem de log(107)=16.1. A simulação abaixo retorna 15.2, bem próximo do valor que foi predito.
Simulação monte carlo do valor esperado, em C++
No fim das contas, a complexidade do algoritmo com tentativa e erro é apenas O(log n), se você tiver um tabelão de primos. Na prática, esse é o método usado por todos que precisam de primos aleatórios: a libgcrypt usada no gpg, por exemplo, utiliza esse método na função gen_prime(), com vários truques pra tornar o teste de primalidade bem rápido.
quinta-feira, 5 de junho de 2008
Aproximações
Ao longo da vida, encontramos aproximações a todo momento. Na escola, a gravidade é aproximadamente 10 m/s2, e a velocidade da luz é aproximadamente 3x108 m/s. Segundo a Bíblia, pi é aproximadamente três (1 Reis 7:23). Mas a minha aproximação predileta é uma que os computeiros usam a todo momento: infinito é aproximadamente oito!
De fato, essa aproximação é o que permite aos computadores trabalhar com números negativos, através do complemento de dois. Lembrando a regra, para calcular o oposto de um número qualquer, basta inverter os bits e somar um. Vamos ver na prática como isso funciona, calculando o oposto de 5.
Cinco em binário é 101b, e pode ser escrito também como uma soma de potências de dois: 1+4. Para inverter os bits de 5, precisamos lembrar que há infinitos zeros na frente do 101b, então o inverso vai ter infinitas potências de dois:
Esse valor, y=~x+1, não se parece com -5. Mas as coisas ficam mais claras se você multiplicar y por dois, parcela a parcela, ...
... e depois subtrair esse valor do original:
Todos as parcelas maiores que 16 cancelam, e do lado de cá, 2y menos y é o próprio y. Então y=-5, QED. Quando calculamos complementos de dois no computador, usualmente fazemos as contas apenas em um byte, mas, no fundo, é a mesma coisa: ao invés de fazer a conta com infinitas parcelas, você aproxima o valor por apenas 8 parcelas.
Seu professor de Cálculo 4 certamente deve ter te avisado dos horrores de manipular seqüências divergentes, que invariavelmente levam a paradoxos (como somar apenas parcelas positivas e obter um resultado negativo). No entanto, às vezes até os paradoxos têm utilidades práticas :)
De fato, essa aproximação é o que permite aos computadores trabalhar com números negativos, através do complemento de dois. Lembrando a regra, para calcular o oposto de um número qualquer, basta inverter os bits e somar um. Vamos ver na prática como isso funciona, calculando o oposto de 5.
Cinco em binário é 101b, e pode ser escrito também como uma soma de potências de dois: 1+4. Para inverter os bits de 5, precisamos lembrar que há infinitos zeros na frente do 101b, então o inverso vai ter infinitas potências de dois:
x = 1 + 4
~x = 2 + 8 + 16 + 32 + ...
~x+1 = 1 + 2 + 8 + 16 + 32 + ...
Esse valor, y=~x+1, não se parece com -5. Mas as coisas ficam mais claras se você multiplicar y por dois, parcela a parcela, ...
y = 1 + 2 + 8 + 16 + 32 + ...
2y = 2 + 4 + 16 + 32 + ...
... e depois subtrair esse valor do original:
2y = 2 + 4 + 16 + 32 + ...
y = 1 + 2 + 8 + 16 + 32 + ...
2y-y = -1 + 4 - 8
y = -5
Todos as parcelas maiores que 16 cancelam, e do lado de cá, 2y menos y é o próprio y. Então y=-5, QED. Quando calculamos complementos de dois no computador, usualmente fazemos as contas apenas em um byte, mas, no fundo, é a mesma coisa: ao invés de fazer a conta com infinitas parcelas, você aproxima o valor por apenas 8 parcelas.
Seu professor de Cálculo 4 certamente deve ter te avisado dos horrores de manipular seqüências divergentes, que invariavelmente levam a paradoxos (como somar apenas parcelas positivas e obter um resultado negativo). No entanto, às vezes até os paradoxos têm utilidades práticas :)