Dez anos de MUST

É díficil achar alguém na faixa dos 30 anos que trabalhe com computadores e não conheça o MSX, um computador pessoal de 8 bits que fez um sucesso enorme por aqui. No seu auge, entre 85 e 91, estima-se que mais de 400 mil máquinas tenham sido vendidas no Brasil.

Mas embora esse período tenha sido o auge em número de usuários, não foi o ápice na parte técnica. Os mais avançados softwares brasileiros de MSX foram feitos no período de 1996 a 2004, e há dois motivos pra isso. O primeiro é que esse foi o início da internet comercial, então muitos dos manuais e datasheets que antes eram raros, passaram a ser difundidos livremente pelas listas de discussão. O segundo é que essa também foi a época em que começaram a se formar os primeiros engenheiros que cresceram com o MSX na infância, e toda essa criançada estava sedenta pra usar os novos conhecimentos na sua velha maquininha.

Eu mesmo fiz um monte de brincadeiras nessa época, talvez a mais conhecida tenha sido o meu emulador BrMSX. Mas agora, em 2008, um outro software que fiz completa dez anos, e nada mais apropriado que comemorar com um detalhado making of. Apresento a vocês o Music Station, ou como também era carinhosamente conhecido, o MUST:

Vamos voltar uma década no tempo: em 1998, a internet crescia, nasciam os primeiros grandes portais nacionais, e começava a surgir uma prática que hoje em dia é extremamente comum: o compartilhamento de músicas em mp3. Primeiro por FTP e IRC, depois pelo Napster e Audiogalaxy, em pouco tempo os arquivos mp3 estavam por toda a parte.

Naturalmente, quem ainda brincava com o MSX também queria entrar na onda. Mais de uma vez me perguntaram "Ricardo, você não pode fazer um programinha pra tocar mp3 no MSX?". Qualquer criatura com um mínimo de bom senso iria perguntar "pra quê?". Mas bom senso nunca foi o meu forte :)

Eu já sabia de antemão que mp3 no MSX não era viável, o coitado do Z80 não iria dar conta. Mas a essa altura eu já sabia evitar um erro muito comum em desenvolvedores, o de prestar atenção no que o usuário pede, ao invés de focar no usuário quer. Sure, eles estavam pedindo mp3, mas o que eles queriam mesmo era só tocar música no MSX. Se fosse um outro formato qualquer não faria diferença. Como eu tinha acabado de cursar Processamento Digital de Sinais na Poli (com o saudoso prof. Max Gerken), achei que poderia encarar a construção de uma solução customizada pro MSX.

Hora de calibrar os objetivos: eu decidi fazer um player que tivesse a melhor qualidade possível pro MSX, que rodasse em qualquer modelo de MSX, incluindo os nacionais, e que conseguisse espremer pelo menos um minuto de música em uma MegaRAM de 256kb, o modelo mais comum. O primeiro passo, portanto, é ver se o MSX dá conta de tocar um sample com qualidade aceitável.

O MSX usa como gerador de som o chip AY-3-8910 (popularmente PSG), que é basicamente um gerador de ondas quadradas. Pra tocar um Lá central, você especifica a freqüência de 440Hz, acerta o volume apropriado, e se diverte com o barulhinho que ele produz. Se você for um mago do chiptune, dá pra fazer músicas bem legais com o PSG, mas definitivamente ele não foi feito pra reprodução de áudio digital. Para conseguir isso, nós temos entender como o chip funciona por dentro.

Internamente, o gerador de ondas quadradas é só um contador digital, que alterna o estado de um flip-flop depois do término de meio período. O chip tem três geradores que você pode ligar ou desligar de maneira independente através de um mixer. Quem tem experiência com chips TTL sabe que sinais desligados usualmente não são interpretados como nível lógico zero, mas sim como nível lógico um.

Pro nosso ouvido tanto faz, um sinal constante é silêncio sempre, independente dele estar em zero volts ou em cinco volts, mas o segredo da coisa é que o controle de volume é uma etapa analógica feita depois do mixer. Por isso, se você desligar o gerador através do mixer, e mudar bem rapidamente o controle de volume, você pode improvisar qualquer forma de onda que quiser! Como o PSG tem 16 volumes possíveis, isso significa que o MSX pode emular um PCM de 4 bits, mesmo sem ter sido projetado pra isso.

Quando o sinal está em nível lógico um, basta variar o volume!

Agora basta ligar esse PCM num timer de alta resolução e pronto. Quer dizer, exceto pelo fato do MSX não ter um timer de alta resolução. O melhor que ele tem é um timer que dispara no início do retraço vertical, a 60Hz, que é uma freqüência baixa demais pra gente. Nesse caso, o jeito é apelar pra gambiarra criatividade :)

Como o Z80 do MSX não tem cache, nem prefetching, nem nenhuma dessas firulas de processadores modernos, os opcodes sempre levam o mesmo tempo pra rodar. Assim, você pode criar um pseudo-timer simplesmente contando quantos ciclos de clock seu programa leva pra executar. Se nós precisássemos de uma rotina que lesse um byte da memória, mandasse pro PSG, e levasse exatamente 50 clocks pra rodar, uma possibilidade seria:

LD A,(HL) ; 8 clocks

LD A,(HL) ; 8 clocks

INC HL ; 7 clocks

OUT (0A1h),A ; 12 clocks

NOP ; 5 clocks

NOP ; 5 clocks

NOP ; 5 clocks

Pra fazer um timer assim, além de programar em assembly, você ainda precisa resolver um subset sum dos opcodes, o que é extremamente divertido (eu assumo que se você está lendo até aqui, deve achar essas coisas divertidas também :). Note como é importante fazer o padding correto, além de colocar 3 NOPs, eu também repliquei o carregemento da memória, só pra fazer a conta bater certinho nos 50 clocks.

Agora que temos um jeito de reproduzir o som, precisamos escolher a sampling rate. Pra conseguir o objetivo de 1 minuto em 256kb, o áudio certamente terá que ser comprimido. Então a sampling rate precisa ser baixa o suficiente pra poder suportar a descompressão em real time entre cada amostra, e alta o suficiente pra não perder muita qualidade de som. Eu escolhi 11kHz; como o clock do MSX é 3.57MHz, isso dá um total de 324 clocks por amostra. Não é muito, mas alguma compressão dá pra fazer.

Vamos calcular a compressão agora. Pra colocar 60 segundos a 11kHz em 256kb, nós precisamos de, no máximo, 256k*8/60/11k = 3.1 bits por amostra. Nosso sinal original tem 4 bits por amostra, então realmente a compressão é necessária. O Z80 não tem fôlego pra implementar transformadas, então não podemos usar nada de DCT, e nem mesmo Haar.

O jeito é modelar alguma coisa bem simples mesmo, tipo um código de tamanho variável como o código de Huffman. Esse tipo de código funciona tanto melhor quanto menos uniforme for o seu histograma. Pegando uma música típica como teste, chegamos em uma entropia de 2.6 bits; mas, por outro lado, a primeira diferença tem entropia ainda menor, apenas 2 bits!

Vamos codificar a diferença então. Com apenas 324 clocks não dá pra implementar todo o Huffman, mas tem uma maneira de simplificar, que é usando compressão com perdas. Como 99.3% das diferenças estão na faixa de -3 a 3, eu posso saturar nesses valores pra deixar o código mais simples. O código que utilizei foi o abaixo:

3 1111

2 1110

1 01

0 00

-1 10

-2 1101

-3 1100

No final a implementação desse código ficou menor do que eu esperava, apenas 169 clocks. Isso significava que eu tinha 198 clocks sobrando, com a cpu parada. Eu poderia ter aprimorado o codec com uma compressão mais eficiente, mas pensei que seria mais divertido, ao invés disso, usar esses clocks sobrando pra adicionar animações. Nesses clocks sobrando eu acabei colocando um osciloscópio com a forma de onda tocada, um banner animado na parte inferior, e no cantinho ainda sobrou cpu pra colocar um pingüim dançando!

Nesse ponto tudo que faltava era uma tela de abertura. Como no final o programa acabou ficando uma experiência multimídia, eu resolvi que ele seria uma homenagem à série Disc Station da Compile, que primava por seus excelentes demos. Eu chamei então o meu programinha de Music Station, e como os Disc Station originais sempre começavam com uma menina bonitinha, eu pedi ao meu irmão que desenhasse uma pingüinha no mesmo estilo. O desenho original que ele fez foi o abaixo:

Adaptar desenhos no MSX é tão complicado quanto adaptar código, a resolução é de apenas 256x192 e ainda tem o problema do color clash (cada bloco de 8x1 pixels pode ter no máximo duas cores). O pingüim tomando sol no canto da imagem teve que sumir, eu não tinha resolução pra isso. Já na pingüinha, eu tive que usar traços bem grossos no contorno, pra minimizar o clash, mas mesmo assim ainda sobraram algumas partes borradas. A solução foi usar sprites pra cobrir esse borramento.

Consegue achar as diferenças?

Por fim, bastou criar um logo, inspirado no logo original do Disc Station, e o MUST estava finalmente pronto! Pra quem quiser brincar com ele, abaixo tem uma versão em .dsk que pode ser carregada em emuladores, como o blueMSX ou o openMSX, e também o código fonte original.

Imagem de disco para rodar em emuladores
Código fonte original do Music Station, em assembly Z80

Na época, a comunidade MSX adorou o MUST. Como eu esperava, o pessoal começou a compartilhar músicas convertidas, do mesmo jeito que o povo do PC compartilhava mp3; e como eu publiquei a engine do MUST em GPL, um monte de gente reaproveitou o código, sendo que a engine foi usada até em joguinhos. Eu até animei a fazer uma versão com vídeo também, mas isso já é outra história, pra outra ocasião :)

Parabéns pelos dez anos, Music Station!

Agradecimentos à Ila, ao Acidx e ao Sturaro por ajudar na arqueologia digital :)

O Desentortador

Eu confesso que fiquei apreensivo quando soube que teria que passar vários meses na Califórnia. Afinal, Mountain View não é nenhuma San Francisco, eu tinha medo de morrer de tédio enquanto estivesse lá. Mas, felizmente, meu medo esvaiu-se quando eu descobri O Sebo!

O Sebo é um lugar lindo e maravilhoso onde você pode comprar The Diamond Age por apenas 50 cents, e inúmeros outros livros por preços tão baixos quanto esse. Na foto acima, você pode ver uma pequena porção do paraíso, mostrando apenas a parte de sci-fi, e apenas os autores das letras D até L. Sim, O Sebo é um lugar gigante e com diversão garantida pra várias semanas.

Depois de encher a mochila com livros, eu tive vontade de tirar uma foto deles pra mandar pros amigos. Foi então que eu tive um problema. Eu estava tentando tirar a foto à noite, e se eu colocasse a câmera exatamente sobre o livro, o flash estourava, e tudo que saía na foto era um grande borrão branco. Sem flash também não funcionava, ficava escuro demais pra câmera conseguir estabilizar a imagem.

Eu resolvi esse problema inclinando a câmera em relação ao plano do livro. Isso certamente resolve o problema do flash, mas em compensação o livro fica torto na foto. Uma solução simples seria simplesmente esperar amanhecer e tirar a foto do livro com a luz do dia. Mas é claro que optei por uma alternativa mais bizarra :) Resolvi escrever um software que desentortasse a foto do livro automaticamente!

Pra isso, vamos começar com a foto original do livro torto:


Pra começar o desentortamento, a primeira coisa que precisamos descobrir é onde estão as bordas do livro, que vão definir como vai ser a rotação que iremos fazer. Achar as bordas é fácil, basta identificar a fronteira entre livro e não-livro. Para isso, eu usei um algoritmo simples de crescimento de região:


Não ficou a coisa mais perfeita do mundo: está cheio de ruído em volta, e o algoritmo ainda vazou pra dentro do livro (a mesa era marrom, da mesma cor do urso na capa, e ele se perdeu com isso).

A boa notícia é que nenhuma dessas coisas importa! O que eu quero achar são as bordas, que tem uma direção preferencial bem determinada. Se eu utilizar a transformada de Hough, tanto o ruído quanto o vazamento devem sumir, e de fato é o que acontece. Essa é a transformada da imagem acima:


A transformada de Hough leva retas em pontos, então as quatro retas que formam a borda do livro viram quatro pontos bem brilhantes na transformada. É simples agora usar um threshold pra identificar esses pontos, e aplicar a transformada inversa pra achar as retas na imagem original:

Nessa imagem eu tracei linhas vermelhas nas retas identificadas pela transformada. Estamos quase lá, mas um problema da transformada de Hough é que ela acha retas, e não segmentos de retas. O que me interessa na verdade são os vértices do livro, então eu preciso achar as intersecções entre essas retas.

O detalhe é que quatro retas em posição geral determinam seis pontos, e não quatro. No caso da imagem acima, por exemplo, as retas verticais se encontram no ponto de fuga, bem acima do livro. A solução que eu usei foi adotar uma gambiarra heurística de considerar apenas os quatro pontos mais próximos do centro da imagem original. Além disso, também é legal converter os pontos pra coordenadas polares em relação ao centro da imagem, assim podemos ordená-los pelo ângulo.

Tendo os quatro vértices, agora é só aplicar uma transformação afim na imagem que a leve pra um retângulo e o trabalho está terminado né? Quem dera hehe. Se você fizer isso, a imagem fica com uma deformação não-linear: achatada na parte de cima, e esticada na parte de baixo.


Na década de 90, todo guri metido a programador fazia sua própria engine 3D pra imitar o Doom e o Quake; quem passou por essa época reconhece na hora a origem da deformação. O que está acontecendo é um problema de inverse perspective-correct texture mapping. A perspectiva introduz um fator do tipo 1/z na sua textura, e você precisa compensar isso quando for mapear na imagem. Vamos fazer a correção de perspectiva então:


Yatta! Finalmente temos a imagem final, desentortada como queríamos :) Na verdade, nós desentortamos a imagem em um único eixo, dá pra ver claramente que além do pitch essa imagem também precisava de correção no roll, mas isso fica como exercício pro leitor.

Eu implementei todos os passos descritos acima em python, usando a biblioteca PIL. Ele roda rapidinho, uns 3s por imagem apenas. O source é o abaixo:

Source do desentortador em python

Se alguém quiser reaproveitar o código pra fazer alguma coisa útil com ele (hehe), sinta-se à vontade. Salvo indicação em contrário, todos os programas desse blog são disponibilizados sob a Licença Crowley (Do what you want shall be the whole of the License).

PS: Ah sim, como vocês repararam, a partir de agora o blog vai ter desenhos também! Eu gostaria de dizer que se funcionou pro Chrome, deve funcionar pra mim também; mas a verdade é que a Aleph já fazia a mesma coisa lá na década de 80 :)

MapReduce is running

Quando eu era estagiário na década de 90, um dos meus primeiros projetos envolvia computação paralela. O projeto Visible Human tinha acabado de fatiar os primeiros espécimes, e nosso grupo estava doidinho pra fazer ray-tracing do volume de dados. Mas o volume de dados era grande demais pros padrões da época, e única saída foi importar uma Meiko com dez processadores, pra dar conta do trabalho.

(Meiko, aliás, que tem uma história curiosa. Apesar de ter dez processadores, só conseguíamos rodar processos em nove deles, tinha um que não funcionava de jeito nenhum. Depois de muitas ligações pro suporte, resolvemos desmontar a máquina pra ver se tinha alguma coisa errada no hardware, e, surpresa, entre o pente de memória e o conector na motherboard tinha um mosquito preso. Sim, aquela cpu tinha um bug :)

Mas programar diretamente a comunicação entre os processadores, com PVM ou MPI, é trabalho demais. A nossa solução foi escrever uma camada intermediária que abstraía a parte de comunicação e fazia o paralelismo de modo implícito. Nossa biblioteca era customizada pra visualização volumétrica paralela, mas hoje em dia a mesma abordagem é feita de maneira mais genérica, por pacotes como o Hadoop e o MapReduce, sendo que esse último eu uso bastante hoje em dia.

Mas, mesmo paralelizando o processamento, o melhor que uma abordagem dessas consegue é um ganho linear no número de processadores. Se o tamanho do seu volume de dados precisa daqueles prefixos que vêm depois do giga, então até o MapReduce pode demorar um pouco pra rodar. Nesse fim de semana eu encontrei o Randall Munroe, e pedi pra ele ilustrar como é um MapReduce na prática :)


Eu também aproveitei pra filmar o Randall enquanto ele desenhava esse sketch. Como o estilo dele é, hum, minimalista, isso o torna o único cartunista que eu conheço que consegue ser mais rápido que o Aragonés:

Eu podia estar roubando

É isso. Eu podia estar roubando, eu podia estar matando, mas estou aqui escrevendo no blog. Reclamações sobre a periodicidade do blog devem ser enviadas diretamente à Rockstar :)

Ao falar de GTA, sempre tem quem o associe ao terrorismo islâmico. Mas pra mim, o que incomoda mais é a associação imediata do islã com o terrorismo, em oposição à sociedade ocidental civilizada. Na verdade, já houve um período onde acontecia o oposto, enquanto a ciência florescia na cultura islâmica, os cristãos saqueavam e destruíam em nome da fé (mas na época chamavam isso de cruzadas, ao invés de terrorismo).

Foi nessa época em que viveu Al-Khwarizm (cujo nome deu origem às palavras algarismo e algoritmo), e também nesse período é que surgiram as primeiras demonstrações por indução finita. Mas a influência mais direta que a matemática islâmica teve em mim foi, sem dúvida, através das obras do Júlio César de Mello e Souza, que escrevia com o pseudônimo de Malba Tahan.

Malba Tahan se inspirava na cultura islâmica pra escrever seus livros de divulgação científica, sendo que o mais famoso deles é o Homem Que Calculava. Num post anterior eu tinha dito que o cálculo do Pi com método de monte carlo tinha sido o segundo puzzle que levei mais tempo pra resolver; pois o primeiro puzzle foi justamente tirado de um dos capítulos do livro: como escrever qualquer número usando apenas quatro quatros. No livro, Malba Tahan mostra as soluções até o 10, que são mais ou menos assim:

0 = 44 - 44
1 = 44 / 44
2 = 4/4 + 4/4
3 = (4 + 4 + 4) / 4
...

Eu tinha só dez anos quando li pela primeira vez, e consegui estender a seqüência até o 20. O posfácio do livro dizia que era possível escrever até o 100, mas eu levei mais de uma década até ver a solução completa!

Depois dos 20 primeiros, eu só consegui fazer algum progresso significativo quando estava no colégio técnico. Eu notei que o problema ficava mais simples se eu reduzisse o número de quatros, montando primeiro todas as soluções com um quatro, depois todas com dois quatros, e assim por diante. Eu não sabia na época, mas o que eu estava fazendo era um forward chaining manual. Mesmo assim, eu só consegui fazer os números pares até 100, minha solução tinha um monte de buracos nos ímpares.

Quando eu estava na faculdade, com os skillz mais afiados, eu larguei mão de fazer manualmente e escrevi um script pra fazer o serviço automaticamente. Eu coloquei no script não só as quatro operações básicas, mas também todas as outras que estão no Concrete Mathematics: raiz quadrada, função piso e função teto, binomiais, potências, fatorial, fatorial crescente e decrescente. A única regra é não usar nenhuma notação que envolva letras (como sin, cos e log; se você permitir log, o problema perde a graça).

Meu script conseguiu finalmente resolver todos os números até 100, e até encontrou múltiplas soluções pra eles! Atribuindo um peso a cada operação, eu mandei ele imprimir somente a soluçao mais simples (por exemplo, adição tem peso baixo, piso e teto tem peso alto). A tabela com as respostas está abaixo, junto com o programa que escrevi na época:

Solução gerada pelo script
Script original, escrito em C

O script pode resolver também os 5 cincos, os 6 seis, e qualquer outra combinação, desde que você esteja disposto a esperar o suficiente. O código original é bem ruinzinho, mas na verdade eu me orgulho de achar que meu código escrito há dez anos é ruim, significa que eu aprendi alguma coisa desde então :)

Como aprender computação

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 :)

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*r². 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(10ⁿ), 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/T

O 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 :)

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%.

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:

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(2ⁿ), 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, é pq², na quarta pq³ e assim por diante. Então o valor esperado é:

X = 1p + 2pq + 3pq² + 4pq³ + ...
X = p (1 + 2q + 3q² + 4q³ + ....)

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 + 4q² + ....))

Agora você tira da cartola y=1+q+q²+q³+... e substitui:

X = p (1 + q(2 + 3q + 4q² + ....))
X = p (1 + q(y + 1 + 2q + 3q² + ....))
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(10⁷)=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.

Aproximações

Ao longo da vida, encontramos aproximações a todo momento. Na escola, a gravidade é aproximadamente 10 m/s², e a velocidade da luz é aproximadamente 3x10⁸ 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:

 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 :)

Python one-liners são Turing-complete

Quem programa em C há décadas normalmente não se dá conta de quão ilegíveis são as expressões mais comuns da linguagem. Quando eu era um garoto recém-saído do BASIC, eu lembro de ter me assustado com coisas básicas como for(i=0; i<10; i++). Mas isso é idiossincrasia do C, outras linguagens não sofrem disso, como o Python.

Python foi planejada para ser legível. Os programadores mais experientes citam o Zen of Python, que dita que "belo é melhor que feio", e "legibilidade conta". De fato, é até difícil escrever código ilegível em Python. Mas é claro que difícil não é impossível, e se o dr. Ian Malcolm fosse um programador, ele certamente diria que "obfuscation finds a way."

Aconteceu comigo semana passada: eu olhava alguns exercícios sobre listas para iniciantes, e notei que, embora eles fossem de fato muito simples, ficariam bem mais divertidos se eu tentasse resolvê-los usando apenas uma linha em cada. Abusando de programação funcional, eu consegui fazer os dez primeiros assim:

Soluções dos exercícios em uma linha de Python cada.

Depois de brincar algum tempo com one-liners, a pergunta que naturalmente se apresenta é: será que é possível fazer qualquer programa em uma linha de Python? A dificuldade vem do fato de que o Python diferencia statements de expressions, e você só pode ter um statement por linha. Em Python, statements incluem print, if, while, for e atribuições, ou seja, um one-liner só pode usar um único desses.

Então, colocando a pergunta de outra maneira: é possível demonstrar que um programa em Python com um único statement é Turing-complete? Existem dois caminhos pra demonstrar isso, o primeiro é construir um emulador para uma máquina de Turing universal em uma linha, o segundo é mostrar que é possível converter para uma linha de Python todos os programas possíveis de um sistema que seja Turing-complete, como o cálculo lambda, ou os tag-systems.

Eu resolvi abordar o problema com a filosofia Ricbit: se existem várias maneiras equivalentes de fazer alguma coisa, escolha a mais bizarra! Assim sendo, vou demonstrar que Python one-liners são Turing-complete através de redução ao Brainfuck (cuja universalidade já foi demonstrada várias vezes).

Vamos lá então: o estado de um programa em Brainfuck pode descrito em qualquer momento por uma quádrupla (mem, p, stdin, stdout), que são respectivamente a memória, o ponteiro, a entrada e saída. Vou implementar cada operação do Brainfuck como funções que recebem quádruplas e retornam quádruplas, descrevendo assim a transição de estado.

A operação mais simples é o ponto, que só adiciona o elemento apontado na saída:

dot = lambda mem, p, stdin, stdout: (mem, p, stdin, stdout+[mem[p]])

Para implementar a vírgula, eu preciso primeiro de alguma maneira de modificar um único elemento de uma lista. Se eu pudesse usar atribuições, bastaria algo do tipo mem[p]=value, mas como atribuições em Python são statements, preciso de uma função auxiliar. Além disso, eu preciso fazer o pop() do valor frontal da lista que guarda o stdin, o que me leva à outra auxiliar:

change = lambda mem, pos, value: [value if i==pos else a for i, a in enumerate(mem)]

get = lambda s: (s[0], s[1:]) if len(s) else (0,[])

comma = lambda mem, p, stdin, stdout: (lambda now, next: (change(mem, p, now), p, next, stdout))(*get(stdin))

Tendo a função change em mãos, fazer os comandos de mais e menos é simples:

plus = lambda mem, p, stdin, stdout: (change(mem, p, mem[p]+1), p, stdin, stdout)

minus = lambda mem, p, stdin, stdout: (change(mem, p, mem[p]-1), p, stdin, stdout)

Os comandos de esquerda e direita precisam tomar o cuidado de aumentar os limites da memória se necessário, a universalidade do Brainfuck requer uma fita infinita para os dois lados:

left = lambda mem, p, stdin, stdout: ([0]+mem if not p else mem, 0 if not p else p-1, stdin, stdout)

right = lambda mem, p, stdin, stdout: (mem+[0] if p==len(mem)-1 else mem, p+1, stdin, stdout)

Agora chegamos na parte complicada, que é o operador de loop. Como for e while são statements, e lambdas recursivos precisam de uma atribuição (fat = lambda x: 1 if x<=1 else x*fat(x-1)), então a única saída é apelar pra lazy evaluation, que no Python é implementada no módulo itertools. (Incluir o módulo itertools poderia tornar o programa um two-liner, mas felizmente é possível importar um módulo usando uma expression ao invés de um statement: a função __import__).

A solução para o operador de loop é criar uma lista infinita contendo [x, f(x), f(f(x)), f(f(f(x))), ...], onde cada f é uma aplicação do conteúdo do loop. Depois, para executar o loop, basta iterar nesta lista infinita, procurando o primeiro elemento onde o elemento apontado pelo ponteiro seja nulo. Precisamos então de uma função que calcule f^n e uma que gere a lista infinita:

composite = lambda f, n: lambda x: reduce(lambda a, b: b(*a), [f]*n, x)

infinite = lambda f, x: itertools.imap(lambda n: composite(f, n)(x), itertools.count(0))

Depois, basta criar um predicado que avalie quando o loop deve parar, e pegar o primeiro elemento da lista onde o predicado é verdadeiro:

predicate = lambda mem, p, stdin, stdout: mem[p] != 0

getfirst = lambda it: [i for i in itertools.islice(it, 1)][0]

loop = lambda f: lambda *x: getfirst(itertools.dropwhile(lambda x: predicate(*x), infinite(f,x)))

Tendo todos os comandos, só precisamos de uma função extra para encadeá-los, e depois, só para o programa não ficar grande demais, um shortcut que executa strings diretamente em Brainfuck:

chain = lambda f: lambda *x: reduce(lambda y, g: g(*y), f, x)

bf = {'+':plus, '-':minus, '.':dot, ',':comma, '<':left, '>':right}

run = lambda f: chain([bf[i] for i in f])

Feito! Agora é só fazer um script para parsear o Brainfuck original e gerar o one-liner. A título de ilustração, esse é o Hello World gerado pelo script:

print ''.join(chr(i) for i in ( (lambda itertools: (lambda change, get, chain, composite: (lambda comma, dot, plus, minus, left, right, infinite, predicate, getfirst: (lambda bf, loop: (lambda run: (chain ([run ("++++++++++"), loop (run ("<+++++++<" "++++++++++<" "+++<+>>>>-")), run ("<++.<" "+.+++++++..+++." "<++.>>" "+++++++++++++++" ".<.+++." "------.--------" ".<+.<.")])) ([0],0,[],[]) )( (lambda f: chain([bf[i] for i in f])) ) )( ({'+':plus, '-':minus, '.':dot, ',':comma, '<':left, '>':right}), (lambda f: lambda *x: getfirst(itertools.dropwhile(lambda x: predicate(*x), infinite(f,x)))) ) )( (lambda mem,p,stdin,stdout: (lambda now,next: (change(mem,p,now),p,next,stdout))(*get(stdin))), (lambda mem,p,stdin,stdout: (mem,p,stdin,stdout+[mem[p]])), (lambda mem,p,stdin,stdout: (change(mem,p,mem[p]+1),p,stdin,stdout)), (lambda mem,p,stdin,stdout: (change(mem,p,mem[p]-1),p,stdin,stdout)), (lambda mem,p,stdin,stdout: ([0]+mem if not p else mem, 0 if not p else p-1, stdin, stdout)), (lambda mem,p,stdin,stdout: (mem+[0] if p==len(mem)-1 else mem, p+1, stdin, stdout)), (lambda f,x: itertools.imap(lambda n: composite(f,n)(x), itertools.count(0))), (lambda mem,p,stdin,stdout: mem[p] != 0), (lambda it: [i for i in itertools.islice(it, 1)][0]) ) )( (lambda mem,pos,value: [value if i==pos else a for i,a in enumerate(mem)]), (lambda s: (s[0],s[1:]) if len(s) else (0,[])), (lambda f: lambda *x: reduce(lambda y,g: g(*y), f, x)), (lambda f,n: lambda x: reduce(lambda a,b:b(*a),[f]*n,x)) ) )(__import__("itertools")) )[3])

QED, em uma única linha, como prometido! (Eu não prometi que seria uma linha pequena :)

O script que converte de Brainfuck para Python one-liner está abaixo, para quem quiser brincar:

Conversor para Python one-liner.

Como nota final, vale lembrar que só porque você pode escrever qualquer coisa em uma linha, não significa que você deve fazer isso. Legibilidade conta :)

Potências ótimas

Olhando no Google Analytics, eu descobri que alguém chegou aqui no blog procurando por "como implementar em c++ potências". Se você é essa pessoa, a resposta está abaixo. Se você não é essa pessoa, puxe uma cadeira que o papo é divertido :)


Calcular potências aproximadas em ponto flutuante é trivial, basta incluir a biblioteca <cmath> e usar a função pow, que internamente é implementada como exp(y*log(x)). Mas existem várias aplicações onde você precisa do valor exato da potência, como, por exemplo, durante a criptografia RSA. Nesses casos, uma primeira abordagem pode ser como no código abaixo:

int natural(int x, int n) {
  int result = 1;
  for (int i = 1; i <= n; i++)
    result *= x;
  return result;
}


Esse código funciona, mas existem maneiras mais espertas. Nós, que temos dez dedos, não estamos acostumados a pensar em binário. Mas se fôssemos como os golfinhos, que tem um cérebro maior que o nosso e só duas barbatanas, poderíamos visualizar o expoente em binário e fazer um código assim:

int binary(int x, int n) {
  if (n == 0) return 1;
  if (n == 1) return x;

  int half = binary(x, n/2);
  if (n & 1)
    return half*half*x;
  else
    return half*half;
}

A versão original era O(n), essa é O(log n), uma melhoria significativa. Mas a melhor notícia sobre esse algoritmo é que você não precisa implementá-lo, ele já está pronto na biblioteca padrão do C++. Basta incluir <ext/numeric> e usar a função power. A função da STL ainda recebe como argumento opcional um functor de multiplicação, então você pode implementar com ela exponenciação modular, ou até mesmo potências de matrizes (ela funciona mesmo que sua multiplicação não seja comutativa).

Por outro lado, esse algoritmo é prático, mas não é ótimo. Em alguns casos, existem maneiras mais rápidas de calcular a potência, o primeiro caso onde isso acontece é para n¹⁵. O algoritmo binário precisa de seis multiplicações para resolver o problema, mas existem soluções com apenas cinco:



A melhor maneira de descrever as multiplicações necessárias para calcular uma potência é através de uma addition chain. Uma addition chain é uma espécie de generalização da seqüência de Fibonacci: enquanto no Fibonacci o próximo elemento é a soma dos dois imediatamente precedentes, numa addition chain o próximo elemento é a soma de dois anteriores quaisquer, ou até mesmo a soma de um elemento anterior com ele mesmo (com a restrição de que a seqüência precisa ser estritamente crescente).

Pela regra de formação dá pra perceber que, ao contrário da seqüência de Fibonacci, existem inúmeras addition chains. E melhor ainda, você pode associar uma addition chain com a seqüência de expoentes gerados no cálculo de uma potência. Para o exemplo de n=15 acima, as duas addition chains correspondentes são [1 2 3 6 7 14 15] e [1 2 4 5 10 15].

Achar uma seqüência ótima para calcular potências, então, é equivalente a achar uma addition chain de tamanho mínimo terminada no número que queremos. Infelizmente, eu tenho duas más notícias: o Erdös mostrou que a addition chain ótima não cresce mais lentamente que O(log n), então assintoticamente ela não é melhor que o método binário; e pior, o cálculo da addition chain ótima é NP-Completo. Abaixo eu implementei a addition chain ótima em C++ (usando brute force, então está bem lento):

Implementação da addition chain ótima em C++

É interessante também dar uma olhada nos casos onde o binário perde. No gráfico abaixo, a linha vermelha é o método binário, a linha azul é a addition chain ótima, e a linha verde é log₂n:


Olhando o gráfico, um golfinho certamente perceberia que o desvio do método binário é proporcional à quantidade de dígitos 1 na representação binária do expoente (os picos são em 63, 127, 255, e assim por diante). A demonstração disso é bem simples e está no Seminumerical Algorithms, junto com várias heurísticas para aproximar a addition chain ótima.