O Jogo do Pi

No último post eu falava de variações em jogos, aí lembrei de uma variação divertida em um jogo bem conhecido. Quem assistia Topa Tudo Por Dinheiro certamente deve se lembrar do Jogo do Pim. O objetivo é enumerar os naturais pelo maior tempo possível, trocando os múltiplos de quatro pela palavra "pim":

1,2,3,pim,  5,6,7,pim, 9,10,11,pim...

Parece simples, mas na prática muita gente se confundia e errava antes mesmo de chegar ao quarenta!

Uma variação mais difícil desse jogo foi inventada pelo Juca uns anos atrás: o Jogo do Pi. Dessa vez você ainda precisa enumerar os naturais, mas precisa falar "π" toda vez que passar um múltiplo inteiro de pi. Os primeiros múltiplos inteiros são aproximadamente 3.14, 6.28 e 9.42, então a sequência começa assim:

1,2,3,π,  4,5,6,π, 7,8,9,π...

Parece o mesmo jogo né? Você conta três números, fala π, conta mais três, fala π, e assim por diante. Mas se você seguir essa estratégia, vai perder! Olha o que acontece se você continuar:

..., 16,17,18,π, 19,20,21,π, 22,23,24,25,π,...

A estratégia de contar de três em três falha! Entre 7π e 8π tem quatro naturais ao invés de só três.

Bem, será que não dá pra melhorar a estratégia? Podemos contar quantos números tem entre os múltiplos de pi. Com sorte, o padrão é 3‑3‑3‑3‑3‑3‑3‑4 e depois repete. Nesse caso, o padrão do Jogo do Pi é oito vezes maior que o padrão do Jogo do Pim. Infelizmente, se você fizer as contas, esse padrão também é quebrado após 112π.

Será que devemos procurar um padrão de padrões então? Nessa altura não já compensa fazer as contas na mão, melhor deixar o python achar o período pra gente:

Script em python para tentar achar um período

Más notícias: o script acha padrões maiores e maiores, sem parar. Não tem como concluir se existe ou não um padrão só analisando a saída do programinha.

O jeito de resolver essa dúvida, então, é usando matemática! Se você tem medo de teoria dos números, pule o quadro azul:

Antes de mais nada, vamos colocar nosso problema na notação correta. Nós queremos descobrir se a quantidade de números entre múltiplos inteiros consecutivos de pi formam uma sequência periódica. Vamos primeiro escrever a sequência explicitamente:

Isso gera a sequência que queremos. Agora vamos supor que essa sequência tem um período p. Das duas uma: ou a gente acha o valor de p, ou batemos em uma contradição pelo caminho.

Se nós somarmos todos os primeiros n termos dessa sequência, temos uma soma telescópica e um resultado curto:

Isso vale para todo n, inclusive no caso onde n é um múltiplo inteiro do período p:

Mas, nesse caso, podemos fazer a somatória de outro jeito. Ao invés de somar direto de 0 a kp, podemos somar k vezes o período p. Aí, como a sequência é periódica, um dos termos some:

O termo mais interno da somatória não varia mais com i, então essa somatória está somando k vezes a mesma coisa:

A conclusão é que, quando p é o período da sequência, podemos jogar um k inteiro de dentro para fora do piso. Podemos agora abrir a definição de piso para kpπ:

Dividindo todo mundo por kp, temos:

Agora, para todo real ε > 0, podemos escolher um k grande o suficiente tal que ε seja maior que 1/kp. Logo:

Pelo teorema do sanduíche, concluímos que:

Note que o numerador dessa função é um piso, logo é um inteiro. Sabemos também que p é um inteiro, logo a fração é um racional. Mas pi é irracional, portanto chegamos a uma contradição, e a sequência não tem período.

Ou seja, concluímos matematicamente que o Jogo do Pi não tem período. Outra maneira de dizer a mesma coisa é que o período do Jogo do Pi é infinito, e portanto o Jogo do Pi é infinitamente mais díficil que o Jogo do Pim!

Agradecimentos ao povo do Math Stack Exchange pela ajuda na demonstração

Estimule sua Criatividade #2

Nesse segundo post sobre criatividade, vou contar outra história do tempo do colégio técnico. As aulas regulares do curso de Eletrônica começavam à tarde, mas duas vezes por semana tinha aula de Educação Física às 7am. Entre essas aulas tinha um intervalo de várias horas, onde muita gente voltava para casa para almoçar.

Eu ficava lá pelo colégio mesmo. Nesses intervalos sem aula eu adquiri vários hábitos, como ler ficção científica da Editora Aleph, gibis do Conan e dos X-Men, Diários de Bordo da Frota Estelar Brasileira, livros-jogos do Steve Jackson, jogar RPG (desde Tagmar até AD&D), ir ao fliperama da esquina jogar Street Fighter 2, enfim, aprendi tudo que não presta.

E quando não tinha mais nada para fazer, eu jogava Jogo da Velha.

O Jogo da Velha

Qualquer um que tenha se dedicado um pouquinho ao Jogo da Velha sabe que é um jogo chato. Se o primeiro jogador começar no centro, ele pode forçar um empate. Mas nós éramos sobreviventes do Atari, então conhecíamos o jogo 3-D Tic-Tac-Toe, criado pela Carol Shaw (a programadora mulher do sexo feminino que também foi responsável pelo clássico River Raid).

3-D Tic-Tac-Toe no Atari 2600

O Jogo da Velha 3D é como o tradicional, mas agora você joga em um tabuleiro cúbico. O objetivo ainda é completar uma linha antes do seu oponente. Note que o jogo do Atari era 4x4x4, e não 3x3x3. No jogo 3x3x3, o primeiro jogador sempre pode forçar a vitória, então é uma variação sem graça.

Depois de algum tempo jogando o Jogo da Velha 3D, tivemos uma idéia: se a versão 3D é mais divertida que a 2D, será que o Jogo da Velha 4D não seria mais legal ainda? Nós desenhamos um tabuleiro 4D no papel e, realmente, era divertido mesmo! Nós desenhávamos o tabuleiro a caneta e jogávamos com lápis, assim era só apagar com borracha para jogar de novo. No fim, ficamos tão viciados que os papéis rasgaram de tanto apagar!

Hora da segunda idéia: os alunos de Eletrônica tinham acesso à sala dos computadores, cheio de PC-XTs rodando DOS. Eu aproveitei que tinha acabado de aprender Pascal, e fiz um tabuleiro eletrônico para jogar direto no computador, sem ter que ficar apagando o tabuleiro a cada jogo. (Você pode pegar esse source no github: 4dvelha.pas)

Jogo da Velha 4D, rodando no PC-XT

Essa versão nos divertiu um monte. Eventualmente tivemos outra idéia: por que não jogar o Jogo da Velha 5D? Essa idéia não foi para frente, depois de algumas partidas nós percebemos que o 5D tinha o mesmo problema do 3D, o primeiro jogador conseguia forçar a vitória.

E o 6D? Um dos amigos usou o plotter do pai para desenhar um monstruoso tabuleiro de Jogo da Velha 6D, que só coube em um papel A2. Mas esse também não tinha graça: o jogo era muito esparso, sempre acabava antes que conseguíssemos utilizar todas as seis dimensões.

No fim, o Jogo da Velha 4D era o ponto ótimo, e jogamos essa variação até o fim do colégio. Mas quando eu entrei na faculdade, topei com um problema diferente.

O Jogo da Velha Minimax

A biblioteca da Poli era enorme e cheia de cheia de títulos legais, como o livro de Inteligência Artificial do Peter Norvig. Ali eu conheci o algoritmo minimax para jogos competitivos, e naturalmente bateu a vontade de implementá-lo no meu Jogo da Velha 4D. (A versão com minimax também está no github: 4dvelha2.pas)

Infelizmente, o branching factor do Jogo da Velha 4D é muito alto: começa em 256, quase a mesma dificuldade do jogo de Go. Nos computadores da época (isso foi por volta de 1994), eu conseguia no máximo descer um único nível da árvore. Eu fiquei bastante frustrado com isso. Será que eu conseguiria bolar alguma variação do Jogo da Velha que tivesse um branching factor menor?

Alguns anos depois eu consegui uma solução: ao invés de mudar o número de dimensões, eu mudei a geometria do tabuleiro! Se você pegar um Jogo da Velha 3D, esticar as pontas, e depois torcer o tabuleiro sobre si mesmo, o resultado é um Jogo da Velha 3D não-Euclideano.

Essa versão eu implementei como um applet Java, então você pode jogar online, ou pegar o source no github. Com o branching factor menor, deu para colocar três níveis de dificuldade, desde o Easy (abre um nível da árvore), até o Hard (abre três níveis da árvore).

Jogo da Velha não-Euclideano

Ganhar do computador no Hard não é impossível, mas você vai precisar de muito treino para conseguir :)

A Criatividade

Dessa vez quem cantou a moral da história foi o Douglas Hofstadter. Em um dos seus livros ele conclui que "variações sobre um tema são o ponto crucial da Criatividade". Pode parecer anti-intuitivo para a maioria: ser criativo não é o mesmo que ser original? Então como pode a criatividade estar baseada na variação?

O problema é que quem está de fora só vê uma pessoa girando um botãozinho: esse Jogo da Velha está na posição 2D, eu giro o botão aqui e tenho 3D, giro mais e tenho 4D. Mas criatividade não é girar o botão, criatividade é achar um botão que não estava lá! Quando eu notei que tinha um botão que girava a geometria, aí sim eu estava sendo criativo.

Em resumo, para expandir a sua criatividade você precisa achar variáveis onde antes só haviam constantes. E isso vale tanto para arte quanto para computação :)

Estimule sua Criatividade #1

Um tempo atrás me perguntaram se existe algum método para ser mais criativo. Isso é bem difícil: para testar se um método funciona ou não, você precisaria medir a criatividade antes e depois de usar o método, e ninguém sabe como mensurar criatividade.

Dito isso, eu conheço dois truques que parecem funcionar bem. Neste post eu vou falar do primeiro método, e pra isso eu vou contar um história que aconteceu comigo lá nos antigos tempos de colégio técnico.

A Federal

Em 1991 eu comecei o curso técnico em eletrônica na ETFSP, popularmente conhecida como a Federal (e que hoje em dia é IFSP). A primeira coisa que nos passaram foi a lista de material para comprar. Além dos básicos livros e cadernos, a lista também incluía itens mais específicos, como caneta nanquim, normógrafo e calculadora científica.

A calculadora foi um problema. Naquela época, calculadoras científicas eram caras (ou melhor, até hoje em dia elas são caras). Como eu não tinha condições de comprar uma, resolvi encarar o desafio de fazer o curso usando só uma calculadora de quatro operações.

O primeiro ano foi fácil: no início do curso, o único componente usado era o resistor, e contas com resistores você consegue fazer só com as quatro operações. No segundo ano apareceram os capacitores, e aí a coisa complicou: agora as contas envolviam exponenciais e funções trigonométricas, que a minha calculadora simples não fazia.

O que fazer então? Dado que eu não conhecia ninguém para me emprestar uma calculadora científica, nem tinha dinheiro para comprar uma, o jeito foi apelar para a herança do meu avô.

A herança do meu avô não era dinheiro. Melhor que isso, era conhecimento! Mais especificamente, era uma estante enorme que ia de parede a parede, com um monte de livros sobre todo tipo de assunto. Os livros eram antigos, da década de 60, e iam desde o Tesouro da Juventude até a Gramática do Jânio Quadros.

O livro que me salvou foi uma curiosa coleção de Matemática para Seminaristas. Eu não sei exatamente por que um padre precisa saber matemática, mas o fato é que o livro era bem completo, começava com produtos notáveis e avançava até o Cálculo. E no meio do caminho, tinha o truque que eu precisava para usar minha calculadora!

A Exponencial

Os problemas que caíam na prova eram mais ou menos assim: dado o circuito RC-série abaixo, calcule a tensão no capacitor após 2 segundos.

No colégio técnico você não precisava resolver a equação diferencial. Ao invés disso, eles davam a fórmula final; você só precisava decorar e aplicar:

O meu problema era a exponencial, que só tem em calculadoras científicas. Mas, usando o livro dos seminaristas, eu descobri a fórmula que resolveu meu problema: a expansão em série de Taylor:

Parece complicado, mas é simples calcular essa fórmula numa calculadora de quatro operações, desde que ela tenha os botões de memória (M+, M-, MC e MR). No exemplo dado, em que eu preciso calcular exp(-0.2), a sequência de botões é curta:

C MC 1 M+ . 2 M- × . 2 ÷ 2 = M+ × . 2 ÷ 3 = M- × . 2 ÷ 4 = M+ MR

Com 29 toques eu tenho o valor de 0.818, correto com três casas decimais (e o professor só pedia duas). Se você tiver os dedos ágeis, praticamente não leva tempo!

Esse método também é bom para calcular senos e cossenos, que aparecem no cálculo de ângulos de fasores e no cálculo do fator de potência. Mas tem um caso onde o método não funciona...

O Logaritmo

Às vezes, o que caía na prova era o problema inverso: dado aquele mesmo circuito RC-série, quanto tempo leva pro capacitor ficar com metade da tensão da fonte? Nesse caso, a conta depende de log(0.5), e não tem logaritmo na calculadora de quatro operações.

Pior, nesse caso o livro também não ajudava. A única fórmula que tinha no livro era a seguinte:

Se você tentar usar essa fórmula na calculadora, logo vai perceber que ela não é prática: são necessários muitos termos para conseguir a precisão de duas casas decimais desejada. Olhando as fórmulas, a velocidade de convergência é clara: os coeficientes da exponencial caem com a velocidade do fatorial, enquanto que os coeficientes do logaritmo caem com a velocidade da série harmônica, que é muito mais lenta.

Hoje em dia, eu provavelmente usaria o método de Newton-Raphson para calcular o logaritmo. Mas eu não entendia de cálculo numérico, então nem sabia que esse método existia. Por isso, eu tive que apelar. Ao invés de decorar uma fórmula, eu decorei quatro números:

É fácil decorar quatro números de 7 dígitos né? Certamente todo mundo sabe pelo menos quatro números de telefones de cabeça, e a dificuldade é a mesma. Com esses quatro números, eu conseguia calcular tudo que queria! Precisa de log(0.5)?

Precisa de log(4.2)? Beleza:

E se for log(2.6)? Ops, aí não dá, 26 tem um fator primo 13 que não está na minha lista. Mas eu posso aproximar:

O resultado não é exato, mas tem precisão de duas casas. Quem precisa de calculadora científica, afinal?

Isso foi sorte, ou sempre é possível aproximar? Na época eu tinha uma intuição que sempre dava, mas hoje em dia eu consigo demonstrar que isso é verdade! A demonstração está na caixa azul, pule se você não sabe teoria dos números:

A demonstração usa diretamente o teorema da equidistribuição de Weyl: o conjunto das partes fracionárias dos múltiplos inteiros de um número irracional qualquer é denso em [0,1). Ou seja, para qualquer irracional α e real q, onde 0 ≤ q < 1, e para qualquer ε positivo dado, sempre existe um inteiro n tal que a fórmula abaixo é verdadeira:

Com um pouquinho de álgebra é fácil estender o domínio de [0,1) para todos os reais. Imagine que escolhemos um real x qualquer, tal que x=p+q, onde p é a parte inteira e q é a parte fracionária. Então podemos partir da equação anterior:

Note que o piso de nα é um inteiro, e p é um inteiro. Vamos chamar de m a diferença dos dois:

Ou seja, eu posso escolher qualquer irracional α e qualquer real x, que sempre vai existir um par de inteiros m e n que satisfazem a inequação. Já que eu posso escolher qualquer número, então escolho α e x como abaixo:

Note que α é irracional, portanto podemos usar nossa inequação:

Ou seja, sempre podemos aproximar o log de um real y qualquer como a soma de múltiplos inteiros de log(2) e log(3), nem precisava ter decorado o log(5) e log(7)!

Infelizmente, essa demonstração é um exemplo de prova não-construtiva: ela garante a existência dos inteiros m e n, mas não fala como achá-los! No fim, você precisa descobrí-los pela intuição, exatamente como eu fazia :)

A Criatividade

Eu poderia terminar o post com uma moral do tipo "se a vida te deu limões, então arranje uns bastões de cobre e faça uma bateria", mas tem uma lição mais importante! A Marissa Mayer cantou a bola em uma entrevista de 2006: a criatividade adora restrições, especialmente se acompanhada de um desprezo saudável pelo impossível.

Tente se lembrar de alguma coisa que você viu e achou muito criativo. Que tal a artista que fazia retratos usando fita cassete? O cara que faz arte usando frutas? O doido que reescreveu The Raven, do Poe, de modo que o número de letras do poema batesse com os dígitos de pi?

Todos eles são exemplos onde o autor se auto-impôs uma restrição (de forma, de conteúdo, de matéria-prima). Se você quer expandir sua criatividade, impor limites costuma ser mais efetivo que retirá-los. Eu mesmo sempre tive problemas quando me pediam uma redação "com tema livre", era muito mais fácil quando o tema tinha alguma restrição.

E isso não vale apenas para arte, vale para computação também. Tente fazer aquele seu programa usar menos memória, menos CPU, tente programar a mesma coisa em menos tempo: tudo isso vai te forçar a usar soluções mais criativas.

Eu ainda tenho outro truque para estimular a criatividade, mas esse fica para o post seguinte :)

(Agradecimentos ao povo do Math Stack Exchange pela ajuda na demonstração)

O Bug mais Difícil

Essa é uma pergunta clássica em entrevistas: qual foi o bug mais difícil que você já consertou? No meu caso, o bug mais difícil que resolvi deu um trabalhão para encontrar, e por um motivo bastante curioso: o bug não era no programa, mas sim na pecinha entre a cadeira e o teclado!

A história

Essa história se passa em agosto de 1998. Nessa época, eu estava escrevendo um emulador de MSX, e tomando dois cuidados na implementação: o emulador tinha que ser rápido o suficiente para rodar em full speed nas máquinas da época, e tinha que ser preciso o suficiente para ser indistinguível de um MSX real.

O primeiro objetivo eu atingi escrevendo o emulador inteiramente em Assembly, o segundo através de uma metodologia de testes sólida. Na época eu ainda não fazia unit testing, mas, para cada feature adicionada, eu sempre fazia um programa de teste e comparava a execução dele no emulador e na máquina real.

As primeiras versões do emulador eram mudas, eu emulava só a CPU e o chip de vídeo. Quando esses dois ficaram estáveis, eu comecei a escrever a emulação do chip de audio: o AY-3-8912 da General Instrument, também conhecido como PSG. Esse chip tinha capacidade de produzir três canais de som e um canal de ruído, sendo que o ruído era usado para fazer efeitos sonoros (como tiros e explosões), e também para fazer sons de percursão (como baterias).

Para testar se a minha implementação do canal de ruído estava boa, eu fiz o seguinte programinha em MSX BASIC:

10 SOUND 7,183
20 SOUND 8,15
30 FOR I=8 TO 1 STEP -1
40 SOUND 6,I
50 PRINT I
60 FOR J=0 TO 300: NEXT J,I
70 BEEP

Os números mágicos no código assustam, mas o funcionamento é simples: o registro 7 é o mixer, com o valor de 183 eu ligo o gerador de ruído no canal A e desligo todo o resto. O registro 8 é o volume do canal A, que eu setei para 15, ou seja, volume máximo. Por fim, eu faço um loop alterando o valor do registro 6, que é a frequência do ruído.

Portanto, o que se espera desse programinha é que ele toque um barulhinho de volume constante, começando grave e ficando cada vez mais agudo. Isso na teoria. Na prática, aconteceu outra coisa!

O bug

Eu resolvi repetir o experimento da época e filmar. No video abaixo você pode ouvir esse programinha rodando em um MSX real, e depois rodando no meu emulador (eu usei o dosbox para compilar a versão de 1998, antes do bug ser corrigido):




Olha que curioso, embora eu tenha programado o volume para ficar constante, no MSX real o volume abaixa sensivelmente no final. E pior, no emulador isso não acontece!

A diferença é sutil, será que eu estava viajando? Um jeito de confirmar é capturando os dois sinais e jogando no Octave para comparar (quer dizer, na época ainda não tinha o Octave, então eu usava o Matlab mesmo):

De olho parece realmente que o turboR vai diminuindo, mas é difícil comparar. Um jeito mais eficiente é notando que a grandeza que perceptualmente nós sentimos como volume, fisicamente é causada pela potência do sinal. E potência a gente sabe calcular, é a integral do quadrado da forma de onda. Jogando isso no Octave:

Olha lá! Não era viagem não, nos dois últimos passos realmente o volume do turboR fica menor que o do emulador! Mas por que isso acontece? Talvez eu tivesse alguma pista analisando o hardware com cuidado.

O circuito

Para entender o funcionamento do PSG, não adianta usar o esquemático do MSX e nem pegar o datasheet do AY-3-8912, porque nenhum dos dois explica o que acontece dentro do chip. Em teoria daria para abrir o chip e tentar ler a pastilha com um microscópio eletrônico, mas como eu não tinha acesso a um, o jeito foi apelar para engenharia reversa mesmo.

Os canais de som do PSG foram fáceis de modelar. Ele tem três geradores de onda quadrada, e cada um funciona da seguinte forma:

O clock de entrada do contador é o clock do MSX dividido por 16, ou seja, 223506Hz. A cada pulso, o contador incrementa de um, e o resultado é comparado com o valor do registro que determina a frequência. Quando o valor é igual, o contador reseta e um flip-flop tipo T inverte de estado, gerando assim a onda quadrada desejada.

Mas e o gerador de ruído? Na época ninguém sabia como funcionava, uns chutavam que era ruído térmico de resistor, outros que era shot noise na junção de algum semicondutor.

Mas a verdade era muito mais simples, o PSG tinha um gerador pseudo-aleatório! Para ser mais exato, tinha um gerador de onda quadrada igual ao descrito acima, e a cada transição na saída ele fazia uma iteração em um Linear Feedback Shift Register.

Nós sabemos que todo LFSR gera um sinal periódico, então, capturando um período inteiro do gerador, dá para deduzir qual é o polinômio característico dele. Fazendo o experimento, o polinômio era o maximal de tamanho 17, ou seja, x¹⁷+x¹⁴+1.

Eu conhecia a forma de onda gerada, e conhecia até o exato LFSR que era usado. Ainda assim, isso não explicava por que o volume abaixava. Será que eu tinha implementado errado?

A implementação

Um dos meus requisitos era rodar o emulador a toda velocidade nas máquinas da época. Em 1998, isso significava emular esse circuito inteiro em um Pentium 1, e ainda sobrar tempo para emular a CPU e o VDP. A coisa era tão apertada que simplesmente escrever em Assembly não era suficiente, eu tinha que garantir que o inner loop rodasse inteiramente em registros.

O problema é que os meros 7 registros de uso geral do x86 não eram suficientes. Minha primeira solução foi apelar feio: durante a emulação do PSG, eu desligava as interrupções e usava o stack pointer como registro! Na minha casa isso funcionava bem, mas em casa eu usava DOS. A maioria dos meus usuários à época usava Windows 95, e aí essa abordagem de mexer no esp fazia a máquina travar.

Vamos para a solução dois então. Tinha um monte de variáveis que mudavam a cada nota tocada pelo PSG, mas eram constantes durante o inner loop. Sendo assim, era uma boa oportunidade para usar código auto-modificável: eu "recompilava" o inner loop cada vez que alguém mudasse um registro do PSG. Essa solução funcionava até no Windows e era rápida o suficiente.

Porém, ela era rápida, mas bastante complexa. Você pode ver esse código no github, não é a coisa mais legível do mundo, nem a mais fácil de entender. Então, seria de se esperar que o volume não estivesse abaixando por conta de algum erro de implementação.

Mas após pensar um pouco, eu concluí que a minha implementação deveria estar correta. Se ela estivesse errada, então esse erro deveria aparecer só no meu emulador. Mas na época exisitiam outros emuladores, e em todos os outros esse mesmo problema aparecia.

Se várias implementações diferentes tinham o mesmo problema, então o erro não era no código, tinha que ser conceitual. E se é conceitual, então talvez eu conseguisse alguma pista analisando matematicamente o sinal.

A transformada

Hora de fazer a modelagem do sinal então! Vamos considerar que os bits de saída do LFSR formam uma sequência discreta {d_i}, e que a cada bit 1 na entrada nós mandamos para saída uma forma de onda s(t). Dessa maneira, se a largura do pulso for T, então o ruído do PSG pode ser computado como:

Da teoria de processos estocásticos, nós sabemos que a densidade espectral de potência do sinal, nesse caso, vale:

...onde G_d(f) é a densidade espectral de potência da sequência {d_i}, e S(f) é a transformada de Fourier de s(t). Se admitirmos que o ruído em questão seja ruído branco, então o G_d(f) é uma constante. Já a forma de onda s(t) é um pulso retangular, cuja transformada é a função sinc. Portanto, a d.e.p. do sinal tem a forma de um sinc ao quadrado. 

O caso que dá problema no emulador é quando o registro 6 vale 1, ou seja, quando T=8.98µs. Assim, a densidade espectral de potência da saída é a abaixo:

Novamente, aquilo que perceptualmente nós percebemos como volume é a potência do sinal, e, pelo teorema de Parseval, a potência é a área abaixo do gráfico acima.

Depois de olhar um tempo para esse gráfico, finalmente caiu a ficha!

A solução

No fim, o problema todo estava entre a cadeira e o teclado! Mais especificamente, no fato de que tinha um ser humano entre a cadeira e o teclado! O sistema auditivo humano só consegue ouvir freqüências até 20kHz, e esse sinal estava espalhando potência para faixas inaudíveis! De fato, um humano só consegue perceber como volume essa parcela da potência:

A integral do sinc ao quadrado não dá para resolver analiticamente, mas dá para calcular numericamente. Eu posso, por exemplo, fazer uma tabela com a potência audível dividida pela potência total, para cada um dos valores de frequência que eu usei no programinha em BASIC:

Registro 6	Porcentagem da potência audível
8	0.99
7	0.99
6	0.99
5	0.99
4	0.99
3	0.96
2	0.83
1	0.50

Se você comparar com o gráfico de potência do emulador versus turboR, dá pra notar que a porcentagem acima é exatamente o que falta para tornar a volume do emulador igual à máquina real! Adicionando essa atenuação no emulador, ele finalmente ficou fiel ao original como eu queria.

Leitores que estejam prestando atenção, a essa altura devem estar se perguntando: "Mas peraí, se o problema é puramente sensorial, então porque o gráfico de potência medido no Octave mostra a perda de volume?". A resposta é que isso também é um artefato do sistema auditivo humano, mas indiretamente.

Quando eu pluguei o turboR na entrada de microfone do meu notebook, a placa de som fez a captura do audio em qualidade de CD, ou seja, com amostragem de 44100Hz. E os projetistas da placa de som são espertos, eles sabem que antes de amostrar você precisa passar o sinal em um filtro passa-baixas com corte na frequência de Nyquist, para evitar o aliasing. E essa frequência foi calculada exatamente para jogar fora o que os humanos não ouvem, gerando o mesmo problema.

Se você quiser refazer o experimento, coloquei no github os scripts em Octave que usei para gerar as imagens e tabelas do post:

Scripts em Octave para análise de ruído

Por fim, esse bug deu um trabalhão para resolver, e no final desconfio que a diferença era tão sutil que ninguém além de mim notou a diferença :)

A Fórmula de Bhaskara

Um tempo atrás o Karlisson teve uma idéia muito boa: criar uma coletânea com os 1001 algoritmos que você deve implementar antes de morrer. O objetivo é fazer uma implementação de referência, em python, para todos os algoritmos clássicos, desde os simples como o Bubblesort até os complexos como o Edmonds-Karp.

Se você também quiser participar, é só fazer um fork do repositório no github. É um ótimo exercício para quem é iniciante, e para os veteranos é uma boa oportunidade de praticar o hábito de fazer code reviews. Já que no post anterior eu falava de parábolas, vamos aproveitar pra revisar um dos códigos do projeto: uma função que acha as raízes reais da equação de segundo grau:

def bhaskara(a, b, c):
    delta = b ** 2 - 4 * a * c
    if delta < 0:
        return None
    else:
        raizes = []
        m1 = math.sqrt(delta)
        r1 =(-b + m1) / (2 * a)
        raizes.append(r1)
        r2 =(-b - m1) / (2 * a)
        raizes.append(r2)
        return raizes

A função é uma implementação direta da fórmula:

Conferindo com a fórmula, o código parece correto. Porém, olhando com cuidado, ele tem três problemas!

Primeiro problema: API inconsistente

O problema mais imediato do código é que ele faz uma divisão por zero quando a é igual a zero. Mas isso é um bug de verdade? A intenção do código era resolver equações do segundo grau, mas, se a é nulo, então a equação não é de segundo grau, é de primeiro!

Eu interpreto isso como um problema na API. Se a idéia do código é resolver somente equações de grau 2, então, antes de sair fazendo conta, ele deveria verificar se as entradas são válidas (ou seja, se a é não-nulo). Por outro lado, se a intenção é resolver equações de grau 2 ou inferior, então ele poderia retornar [-c/b] quando a for zero. Nos dois casos, faltou sinalizar exatamente qual o propósito da API.

A API ainda tem uma segunda falha, que é o valor de retorno. O que exatamente ele vai retornar? Eu poderia, por exemplo, dizer que a função sempre retorna uma lista com as soluções existentes. Nesse caso, quando o delta é negativo, o correto seria retornar uma lista vazia, ao invés de retornar None.

Além disso, o que devemos retornar quando o delta é exatamente zero? As duas possibilidades são retornar uma lista com uma única raiz, [-b/2a], ou então retornar uma lista com duas raízes iguais, ou seja, [-b/2a,-b/2a]. Qual resposta faz mais sentido?

A escolha correta vem do Teorema Fundamental da Álgebra, que diz que o número de raízes complexas de uma equação polinomial é igual ao grau da equação. O motivo de considerarmos que algumas equações tem múltiplas raízes iguais é pra fazer esse teorema ser válido em todas as situações.

Mas o teorema só vale para raízes complexas! Como estamos trabalhando com raízes reais, então eu acho que o mais correto é retornar uma lista com um único elemento mesmo.

Você também poderia argumentar que o nome da função não é apropriado, porque está descrevendo a implementação (fórmula de Bhaskara) ao invés de descrever a interface (solução da equação quadrática). Mas aí nós chegamos no segundo problema...

Segundo problema: Bhaskara quem?

Certamente você ouviu na escola que a solução da equação quadrática chama-se fórmula de Bhaskara. Mas, curiosamente, é só no Brasil que a fórmula tem esse nome!

Pode conferir: a wikipedia em inglês nem cita o Bhaskara. Talvez na wikipedia em hindi? Nada de Bhaskara lá também, na Índia eles chamam a equação de fórmula de Sridhar Acharya. Nos outros países a fórmula nem tem nome, é simplesmente "equação quadrática" ou "fórmula a-b-c".

Na verdade, os babilônicos já sabiam resolver alguns tipos de equações quadráticas há pelo menos 4000 anos. Os gregos conheciam soluções geométricas. E a fórmula em si, é do Bhaskara mesmo?

Bem, uma maneira de resolver a dúvida é procurando exatamente o que ele escreveu. O texto mais importante do Bhaskara foi o livro que dedicou à filha, Lilavati. O original era em sânscrito, mas tem uma tradução em inglês online para quem quiser ler. Eu procurei pela solução da quadrática, e o mais próximo que achei foi a proposição 63:

Uma quantidade, incrementada ou decrementada de sua raiz quadrada multiplicada por algum número, é dada. Então adicione o quadrado de metade do multiplicador da raiz ao número dado, e extraia a raiz quadrada da soma. Adicione metade do multiplicador, se a diferença foi dada; ou subtraia, se a soma for dada. O quadrado do resultado é a quantidade procurada.

Em notação moderna, o que o Bhaskara escreveu foi:

Eu não acho que seja correto creditar ao Bhaskara a solução da equação quadrática. A fórmula acima é bem parecida com a que usamos, mas note o detalhe: ela só fornece uma raiz! Para considerar a solução como correta, eu acho que ele deveria ter explicitado que algumas equações vão ter duas raízes.

Pois bem, podemos então mudar o nome da função, de bhaskara() para algo mais apropriado, como quadratic_roots(). Mas ainda resta o último problema...

Terceiro problema: Instabilidade numérica

Vamos testar o código original com a seguinte equação: x^2-2000000000x+4=0. As raízes são aproxidamadamente 2e9 e 2e-9. Porém, olhe só o resultado:

>>> bhaskara.bhaskara(1, -2e9, 4)
[2000000000.0, 0.0]

Como assim zero? Cadê a solução 2e-9? Essa raiz foi vítima de um problema que nem todo mundo presta atenção: muito cuidado com a perda de precisão em ponto flutuante. Em especial, quando você subtrai dois números de magnitude similar, a resposta pode não ter precisão suficiente para dar o resultado correto. Para a equação dada, o código original tem os valores b e m1 muito próximos, e a subtração perde a menor raiz.

A solução é evitar a subtração. Você inspeciona o sinal de b, e escolhe m1 com o mesmo sinal, de modo que os módulos são sempre somados, achando assim a primeira raiz x1. Tendo x1, você pode achar x2 usando a fórmula de Viète para o produto de raízes:

Dessa maneira a subtração problemática some, e chegamos na versão final do programa corrigido:

def quadratic_roots(a, b, c):
  """Find all real roots of a quadratic equation.

  Args:
    a, b, c: Coefficients of the quadratic, a*x*x+b*x+c=0
  Returns:
    A list with all real roots.
  Raises:
    ValueError if a==0 (the equation is not quadratic).
  """
  if not a:
    raise ValueError
  delta = b * b - 4 * a * c
  if delta < 0:
    return []
  if not delta:
    return [-b / (2 * a)]
  if b < 0:
    x1 = (-b + delta ** 0.5) / (2 * a)
  else:
    x1 = (-b - delta ** 0.5) / (2 * a)
  x2 = c / (a * x1)
  return [x1, x2]

Conferindo:

>>> roots.quadratic_roots(1, -2e9, 4)
[2000000000.0, 2.0000000000000001e-09]

Nice!

Você pode se perguntar se esse exemplo que eu dei não é artificial demais, e se coisas assim acontecem na prática. Mas acontecem sim! Um exemplo comum onde esse erro acontece é na implementação da busca binária.

Por exemplo, suponha que você quer resolver a equação ln(x)+sqrt(x)=5. A função ln() é monotônica, o sqrt() também, então dá pra achar a raiz por busca binária:

def search(epsilon):
  inf, sup = 0.0, 1e10
  while sup - inf > epsilon:
    med = (sup + inf) / 2
    x = math.log(med) + math.sqrt(med)
    if x > 5.0:
      sup = med
    else:
      inf = med
  return inf

Vamos testar essa rotina para várias precisões:

>>> binary.search(1e-4)
8.3093709690729156
>>> binary.search(1e-10)
8.3094326941933723
>>> binary.search(1e-15)

Oops! Quando epsilon vale 1e-15, a rotina trava! De novo, o problema é perda de precisão. Nesse caso, as variáveis sup e inf tem valores muito próximos, então a média med não sai do lugar, ficando presa no valor do inf para sempre.

A solução, novamente, é sumir com a subtração. Ao invés de controlar os extremos do intervalo, você usa uma variável para o valor inicial e outra para o tamanho do intervalo:

def smart_search(epsilon):
  inf, delta = 0.0, 1e10
  while delta > epsilon:
    delta /= 2
    med = inf + delta
    x = math.log(med) + math.sqrt(med)
    if x < 5.0:
      inf += delta
  return inf

Conferindo:

>>> binary.smart_search(1e-4)
8.3093709690729156
>>> binary.smart_search(1e-10)
8.309432694193374
>>> binary.smart_search(1e-15)
8.3094326942315693
>>> binary.smart_search(1e-30)
8.3094326942315693

Agora sim!

No fim das contas, é importante você ter a matemática sempre afiada, mas também é igualmente importante conhecer os limites impostos pela computação do mundo real.

(Agradecimentos ao Karlisson, ao Henrique, ao Wladimir, ao Chester e ao povo da lista eng-misc pela ajuda na pesquisa!)