Review: Logicomix

Agosto 12, 2012

Logicomix é uma estória em quadrinhos sobre a busca do fundamento da matemática. É um livro fantástico porque mistura romance, história, mitologia e matemática. Além disso é uma meta-estória, pois relata o processo de desenvolvimento dos quadrinhos na própria estória.

Um dos autores (e personagem) do livro é Christos Papadimitriou, atualmente professor de teoria da computação em Berkeley. Durante minha iniciação científica cheguei a estudar um de seus livros escrito com Ken Steiglitz sobre Otimização Combinatória, mas achei muito difícil e acabei usando outro :P

Teóricos da computação fazem muito uso da lógica, especialmente os que trabalham com complexidade computacional e computabilidade. A estória dá a entender que Papadimitriou foi convidado como o especialista em lógica, embora o autor principal, Apostolos Doxiadis tenha formação matemática também.

Spoiler! O restante deste post contém um breve resumo do livro, além de uma análise pessoal.

Resumo

Tractatus Logico-Philosophicus

Bertrand Russell é o personagem central do livro que descreve sua saga em busca dos fundamentos da matemática.

Russell aparece em dois planos: em um está no presente, palestrando a um grupo de americanos contrários à adesão dos Estados Unidos à Segunda Guerra Mundial. Durante essa palestra ele conta a história de sua vida, que acontece em um outro plano. A grande maioria dos acontecimentos acontece nesse segundo plano e o Russell do plano presente serve mais como um narrador.

Além desses dois planos temos o plano em que os autores do livro estão elaborando a estória ao mesmo tempo em que ela acontece. Eles geralmente aparecem para justificar alguma discrepância entre a história e a estória, além de explicar conceitos mais técnicos.

Logo criança, Russell torna-se órfão e passa a viver com seus avós na mansão Pembroke. É apresentado aos trabalhos de Euclides por um professor particular de matemática.

Russell segue para o Trinity College, onde conhece Alfred Whitehead. Casa-se com Alys Smith e então fazem uma viagem ao continente, passando primeiramente pela Alemanha, onde conhecem Gottlob Frege e George Cantor (embora o livro deixe claro que Russell não os conheceu pessoalmente).

Na França, conhecem a esposa de Whitehead, Evelyn Wade. A estória retrata um interesse de Russell em Evelyn, embora em nenhum momento do livro pareça ter havido algo mais entre eles. Depois, Russell participa do congresso internacional de matemáticos onde David Hilbert divulga a famosa lista de 23 problemas abertos em matemática.

Neste ponto Russell descobre um paradoxo na teoria dos conjuntos, que abala o mundo da matemática. Whitehead e Russell decidem então começar a escrita do Principia mathematica, que se mostra uma tarefa desgastante.

Russell passa a dar aulas em Cambridge e é quando Ludwig Wittgenstein torna-se seu estudante de doutorado. Evelyn, a esposa de Whitehad acha que está prestes a morrer e pede para Russell consolar seu filho, Eric. Cansado de seus trabalhos com a lógica, Russell sente-se reconfortado por exercitar sentimentos como compaixão e amor para com Eric.

Começa a primeira Guerra Mundial. Russell se posiciona contra a adesão do Reino Unido, o que acaba o levando à prisão. Wittgenstein se alista ao exército Autro-Húngaro enquanto Eric, filho de Whitehead, faz o mesmo pelo lado britânico. Wittgenstein sobreviveu à guerra, fazendo-o mudar sua maneira de ver o mundo. Por outro lado, Eric é abatido em combate.

Após o fim da primeira guerra, Wittgenstein estava impedido de ir à Inglaterra por ser austríaco e portanto se encontra com Russell em território neutro: Haia, Holanda. Nesse encontro Wittgenstein contrapõe diversas ideias de Russell.

Russell casa-se pela segunda vez com Dora Black, tendo o seu primeiro filho. Seu casamento era bastante liberal, sendo que ele aceitava que amantes frequentassem sua casa. Seus pais também tiveram um casamento desta forma. Juntos decidiram fundar uma escola progressista para crianças, a Beacon Hill. À mesma época, Wittgenstein começa a dar aulas em escolas primárias.

O matemático Moritz Schlick passa a organizar o Círculo de Viena, uma associação de filósofos que se encontravam em torno da Universidade de Viena. Em uma dessas reuniões, Russell conhece Kurt Gödel e John Von Neumann e foi nesse encontro que Gödel anunciou seu primeiro teorema de incompletude, destruindo o sonho de Russell de definir um conjunto de axiomas que pudesse demonstrar todas as verdades matemáticas. Em seguida, as influências nazistas começaram a se espalhar pelo continente e em 1936, Schlick foi assassinado por um fanático nazista, o que pôs fim ao círculo e é nesse ponto que Russell conclui sua narrativa.

A estória termina voltando para o plano dos autores, que vão assistir à participação de Anne Bardy em uma encenação da peça Oresteia.

Matemática

Das diversas referências matemáticas presentes no livro, gostaria de comentar sobre as três que achei mais interessante.

Paradoxo de Russell

Em 1901, Russel apresenta um paradoxo da teoria dos conjuntos. Esse paradoxo por ser melhor entendido através de um exemplo, o paradoxo dos barbeiros.

Considere uma cidade onde existe apenas um barbeiro e uma lei que diz que o barbeiro deve fazer a barba de todos os homens que não se barbeiam sozinhos, mas não pode fazer a barba de quem se barbeia sozinho.

Consideremos duas possibilidades: o barbeiro se barbeia sozinho ou não. Se ele se barbeia sozinho, está infringindo a lei, porque ele não pode barbear quem se barbeia sozinho. Por outro lado, se ele não se barbeia, a lei o obriga a fazer sua própria barba. Conclui-se que a lei é contraditória.

A generalização para a teoria dos conjuntos é a seguinte: seja R conjunto de todos os conjuntos que não contém a si mesmos. Pergunta-se: R está contido em si mesmo? Vamos analisar duas possibilidades. Se R está contido em si mesmo, então R é um conjunto que contém a si mesmo e não poderia estar em R. Por outro lado, se R não está contido em si mesmo, ele é um conjunto que não contém a si mesmo e portanto deveria ser incluído em R :)

O Infinito e a Teoria dos Conjuntos

Em Logicomix, Russell descreve o paradoxo do hotel de Hilbert que exemplifica o quão contra-intuitivo é trabalhar com o conceito de infinito. A ideia por trás do paradoxo pode ser usada para provar, por exemplo, que existe uma bijeção entre os números naturais pares e os números naturais!

Isso está relacionado com o trabalho de George Cantor, que criou os conceitos de conjuntos enumeráveis e não enumeráveis, sendo os enumeráveis aqueles com mesma cardinalidade dos números naturais, ou seja, para os quais existe um mapeamento um pra um de seus elementos para com os números naturais.

Assim como os números naturais pares têm mesma cardinalidade dos números naturais, Cantor demonstrou que os números inteiros são enumeráveis assim como o conjunto dos números racionais! O conjuntos dos números reais não é enumerável.

Os 23 problemas de Hilbert

David Hilbert divulgou uma lista de 23 problemas não resolvidos à época de 1900. Há três problemas sobre os quais gostaria de comentar.

O principal problema da lista, e que até hoje não foi resolvido, é o problema 8, que é demonstrar a Hipótese de Riemann. O livro The music of primes, fornece uma introdução bem acessível sobre a hipótese, a função zeta de Riemann e a relação com números primos.

Um problema importante também é o número 2, provar que os axiomas da aritmética são consistentes. O segundo teorema da incompletude de Gödel mostra que não é possível provar esse teorema usando a própria aritmética. Entretanto, não há um consenso se isto é uma solução ou não para o problema.

Finalmente, um problema que eu acho particularmente interessante é o número 10: encontrar um algoritmo que determina se uma equação Diofantina polinomial com coeficientes inteiros tem uma solução.

Gosto desse problema porque é fácil de enunciar, está relacionado à área que pesquisei no mestrado (programação linear inteira) e porque foi resultado de um trabalho longo e de contribuição de diversas pessoas, Martin Davis, Yuri Matiyasevich, Hilary Putnam e Julia Robinson.

Mitologia

Algumas referências mitológicas aparecem na estória. A primeira é a um mito indiano sobre a tartaruga que segura o mundo. Russell faz a analogia como o mundo sendo a matemática e a tartaruga seus fundamentos.

Em certo ponto Whitehead mostra a Russell o quadro “As Danaides” de John Waterhouse. As Danaides são as 50 filhas de Dánao, que estavam prometidas aos 50 filhos de Egipto, seu irmão gêmeo. Todas exceto Hipernestra assassinaram seus esposos na noite de núpcias e como punição foram obrigadas por Hades a encherem um jarro de água com furos durante toda a eternidade.

As Danaides: pintura a óleo de John Waterhouse

Whitehead usou o mito como metáfora para o trabalho do Principia Mathematica. Russell queria encontrar um fundamento para a matemática, mas sempre via necessidade de fundamentar cada base que encontrava. Como se disse no livro, ele estava procurando a base da tartaruga que segurava o mundo, mas encontrou uma torre infinita de tartarugas.

Finalmente, temos a encenação peça Oresteia, uma trilogia escrita pelo dramaturgo grego Ésquilo.

O livro retrata a terceira peça, Eunêmides, no qual as Erínias, deusas da vingança decidem punir Orestes. Atena intervém e decide que o futuro de Orestes deve ser decidido por um júri popular.

Literatura e Teatro

Alguns clássicos da literatura e teatro aparecem na estória. O avô de Russell possuia uma biblioteca onde Bertie descobriu a Divina Comédia de Dante Alighieri, que foi chamado de o livro proibido.

Em Cambridge lê Pais e Filhos, de Ivan Turguniev. Segundo a Wikipédia, o protagonista do livro é Bazárov, um jovem intelectual materialista, que nega o amor, a arte, a religião e a tradição, e diz acreditar apenas em verdades cientificamente comprovadas pela experiência. Essa negação da religião e tradição podem ser encontradas na personalidade de Russell na estória.

A estória também retrata a peça de teatro Espectros, do norueguês Henrik Ibsen. Há um destaque para a frase “os pecados dos pais se repetem nos filhos”. Talvez seja uma referência ao estilo de vida dos pais de Russell, que acabou se repetindo em seu casamento com Dora.

Russell passa a recitar versos do poema Alastor, de Percy Shelley. Segundo a Wikipedia, Alastor conta a história de um poeta que persegue a parte mais obscura da natureza em busca de “verdades estranhas em terras desconhecidas”. Podemos enxergar um paralelo aqui com a busca pelos fundamentos da matemática por porte de Russell.

Ao conhecer Alys, Russell cita diversos personagens do livro Alice no País das Maravilhas, de Lewis Carroll, como o Tweedledee, Tweedledum, Gato de Cheshire e o Chapeleiro maluco.

Após o fim da Primeira Guerra Mundial Russell demonstra seu descontentamento com o dadaísmo e cita versos do poema The Second Coming, the William Butler Yeats, que refletem sua opinião.

Conclusão

Gostei bastante deste livro e recomendo a quem gosta de lógica e história da matemática. Eu particularmente não sou muito fã de história em quadrinhos, achei que se fosse em prosa a estória poderia conter mais detalhes. Lendo sobre a vida dos principais personagens na Wikipedia, penso que faltou um pouco de coesão nos fatos relatados.

Achei excelente inserir referências a clássicos da literatura para ajudar a caracterizar a personalidade de Russell. Fico feliz de ter escrito esse resumo e análise porque acho que consegui enxergar alguns detalhes que tinham passados despercebidos na primeira leitura.

Interpretei o fato da estória ser uma meta-estória como uma alusão ao paradoxo de Russell, que tem como base conjuntos que contêm a si mesmos. No caso do Logicomix é uma estória que contém a si mesma.

A encenação de Orestéia me pareceu desconexa com o restante da estória. Talvez a ideia fosse mesmo contar duas estórias independentes (a de Bertrand e a dos autores).

Leitura adicional

Alguns livros e links que eu selecionei para o aprofundamento no conteúdo do Logicomix. Na verdade só li o primeiro e é o único que posso recomendar. Os outros estão na minha pilha de leitura.

[1] The Music of the Primes de Marcus du Sautory – Já li, e gostei muito. Fala sobre os 23 problemas de Hilbert e a hipótese de Riemmann de uma maneira bastante acessível.
[2] Alice no país das maravilhas de Lewis Carroll – Lewis Carroll era um matemático e lógico. Este livro é mencionado no Logicomix.
[3] Scott Aaronson – Filosofia e Ciência da Computação Teórica, Scott Aaronson é um professor de Teoria de Computação no MIT e conhecido por seu trabalho em Computação Quântica.

Anúncios

Paul Graham e Combinatores em Haskell

Agosto 5, 2012

Paul Graham é um ensaísta, programador e investidor. É graduado em Cornell e possui doutorado por Harvard [1].

Como ensaísta, escreveu diversos artigos sendo o mais famosos deles o “A plan for spam” sobre o uso de um filtro Bayesiano para combater spam.

Como programador, Graham é conhecido por seu trabalho com Lisp, tendo escrito livros sobre esta linguagem e vem desenvolvendo um dialeto chamado Arc. Também desenvolveu a primeira aplicação web chamada Viaweb, posteriormente adquirida pela Yahoo!

Como investidor, fundou em 2005 a empresa Y Combinator que faz investimentos em start-ups. Algumas das empresas financiadas pela Y-Combinator são: Dropbox, Airbnb, reddit e Posterous.

Sobre o nome da empresa, eis uma tradução livre do FAQ deles:

Porque vocês escolheram o nome “Y Combinator?”

O combinador Y é uma das ideias mais legais em ciência da computação e é também uma metáfora para o que nós fazemos. É um programa que roda programas; é uma empresa que ajuda a abrir empresas.

Neste post gostaria de apresentar algumas notas de estudos sobre combinadores em Haskell, descrevendo os principais tipos de combinadores incluindo o tipo Y.


Combinadores em Haskell

Esse post teve origem durante os estudos do Capítulo 9 do Real World Haskell.

Em certo ponto do capítulo o autor menciona o termo de combinadores, não deixando muito claro o que são e para que servem. Fui pesquisar um pouco sobre o assunto e decidi escrever um post.

Um combinador pode ser definido como uma função lambda que não tem variáveis livres. Uma variável livre é uma variável que não foi passada como parâmetro na função lambda. Alguns exemplos de funções lambda em Haskell:

1) A função abaixo é um combinador pois a variável x está no parâmetro

\x -> x

2) No exemplo a seguir y é uma variável livre e portanto não é um combinador.

\x -> x y

3) Aqui ambos x e y foram passadas por parâmetro, então temos um combinador.

\x y -> x y

4) O exemplo abaixo não é um combinador porque f não foi passado como parâmetro

\x y -> f (x y)

5) A seguir temos uma pegadinha. Observe que o operador + é na verdade uma função que recebe dois argumentos, como no exemplo (4). Porém como o operador não foi passado como parâmetro, também não é um combinador.

\x y -> x + y 

Simplificadamente, podemos pensar em um combinador como uma função que recebe funções e retorna outra função resultado da combinação das funções de entrada.

Tipos de combinadores

Nas próximas seções, vamos apresentar alguns dos combinadores mais conhecidos. Esses combinadores têm como nome uma única letra maiúscula. Porém, no livro To Mock a Mocking Bird [2], Raymond Smullyman apelida os diferentes tipos de combinadores com nomes de pássaros. Por isso, cada seção terá o nome do combinador e entre parênteses o apelido dado por Smullyman.

Combinador I (Identity Bird ou Idiot Bird)

Este é o combinador mais simples, a identidade. É uma função que retorna o seu próprio parâmetro. Em Haskell podemos escrever:

i x = x

Essa é justamente a definição da função id da biblioteca Prelude.

Combinador K (Kestrel)

Esse é o combinador constante, ele recebe dois argumentos e ignora o segundo. Em Haskell podemos escrever:

k x y = x

Essa é justamente a definição da função const da biblioteca Prelude.

Combinador S (Starling)

O combinador S é uma função que pode ser definida da seguinte maneira em Haskell:

s x y z = x z (y z)

Podemos escrever o combinador I usando apenas os combinadores S e K como

i = s k k

A demonstração é a seguinte:

Aplicando a função (s k k) sobre um argumento x temos

(s k k) x = s k k x

Substituindo s por sua definição,

(s k k x) = k x (k x)

Agora seja f = (k x) uma função qualquer. Temos que k x f = x por definição, então

k x (k x) = k x f = x = i x

É possível mostrar que qualquer combinador pode ser escrito como função dos combinadores S e K!

Combinador B (Bluebird)

O combinador B pode ser escrito em função de S e K como

b = s (k s) k

Vamos ver a cara dessa função aplicada a alguns argumentos. Seja f1 = (k s). Usando currying vemos que essa é uma função que recebe um argumento e sempre retorna s. Aplicando b sobre um argumento x, temos

b x = s f1 k x = f1 x (k x) = s (k x)

Seja f2 = (k x). Temos que f2 é uma função que recebe um argumento qualquer e sempre retorna x. Assim, temos b x = s f2. Vamos aplicá-la sobre mais dois parâmetros y e z:

(s f2) y z = s f2 y z = f2 z (y z) = x (y z)

Ou seja

b x y z = x (y z)

O que esta função faz é a composição das funções x e y e é justamente a definição do operador (.) do Haskell.

Combinador C (Cardinal)

Esse combinador pode ser escrito em função dos combinadores K, S e B como

c = s (b b s) (k k)

Pode-se mostrar que esse combinador equivale a

c f x y = f y x

Essa é justamente a definição do operador flip do Prelude, que retorna uma função com a ordem dos parâmetros trocada.

Combinador Y (Sage Bird)

Simplificadamente o combinator Y é um combinador que transforma uma função não recursiva em uma função recursiva!

Antes de descrevê-lo, vamos analisar como faríamos isso sem o combinador Y. Para isso, vamos usar como exemplo uma função que calcula o fatorial [8]. Fazendo uso de recursão, podemos escrever:

fat = \n -> if n == 0 then 1 else n * fat (n - 1)
fat 5 -- imprime 120

A expressão lambda acima não é um combinador porque depende de fat, que não foi passado como parâmetro.

Vamos tentar passar a função por parâmetro então

fat' = \f n -> if n == 0 then 1 else n * f (n - 1)

Observe que fat’ não é uma função recursiva. Agora fazemos o seguinte:

fat = fat' fat
fat 5 -- imprime 120

Estamos definindo fat como sendo o resultado de fat' quando passamos fat como parâmetro. Parece mágica, mas em linguagens com avaliação preguiçosa isso vai funcionar. Por quê?

Vamos analisar alguns exemplos. A função fat pode ser desmembrada em

fat = 
  fat' fat = 
      \n -> if n == 0 then 1 else return n * fat (n - 1)

Para n = 0, a função retornará 1 e não se dará ao trabalho de calcular fat(-1).

Vamos agora desmembrar a função duas vezes:

fat = 
  fat' fat = 
      fat' (fat' fat) = if n == 0 then 1 else return n * 
      (if (n - 1) == 0 then 1 ele return (n-1)* fat (n - 2))

Para n = 1, é possível ver que uma avaliação preguiçosa não precisa desmembramentos do que isso.

Ok, conseguimos transformar uma função não recursiva fat' em uma função recursiva. Vamos criar uma função auxiliar para que a função fat' seja passada como parâmetro

aux f = f (aux f)
fat = aux fat'
fat 5 -- imprime 120

A função aux faz o que queremos que o combinador Y faça: recebe uma função não recursiva (por exemplo fat') e retorna uma recursiva (nesse caso fat). Entretanto, aux não é um combinador porque o aux que é chamado recursivamente é uma variável livre.

Vamos apresentar então o combinador Y. Sua definição teórica é dada por

y = \f -> (\x -> f (x x)) (\x -> f (x x))

Ainda não consegui entender essa definição de maneira clara o suficiente para arriscar uma explicação aqui. Minha sugestão é ler [8], onde Mike Vaniel deriva a função acima em Scheme.

De qualquer maneira, se formos implementá-lo em Haskell, teremos um erro de compilação devido a tipo infinito [3]. Para resolver esse problema, [4] apresenta a seguinte solução:

newtype Rec a = In { out :: Rec a -> a }

y :: (a -> a) -> a
y = \f -> (\x -> f (out x x)) (In (\x -> f (out x x)))

Para evitar problemas com tipo infinito, encapsulamos o segundo termo (\x -> f (x x)) em um novo tipo, o Rec a, através do construtor In. As funções lambda recebem x que é do tipo Rec a e por isso não podemos aplicar x sobre x diretamente. Antes devemos “desencapsular” a função contida nele usando out.

Lembrando, o newtype define uma estrutura semelhante a um tipo de dado algébrico, com a diferença que só se permite um construtor, no caso definido pelo In, e um “getter”, especificado pelo out.

Por exemplo, vamos definir um tipo chamado NovoTipo que encapsula uma função que recebe um tipo genérico e retorna esse mesmo tipo genérico, ou seja, com tipo a -> a.

newtype NovoTipo a = Construtor { getter :: a -> a }

Podemos encapsular uma função que recebe e retorna um inteiro. Por exemplo:

somaUm::Int -> Int
somaUm v = v + 1
let x = Construtor somaUm

Podemos verificar o tipo dessa variável no ghci:

> :type x
x :: NovoTipo Int

O getter extrai a função que foi encapsulada, ou seja, podemos fazer:

> (getter x) 9
10

Outros combinadores

Este site [5] tem uma lista mais extensa de combinadores. Existe um pacote do Haskell chamado data-aviary [6] contendo a implementação de diversos desses combinadores.

Reg Braithwaite discute alguns combinadores na prática usando Ruby no seu blog/projeto homoiconic [7].

Conclusão

Fiquei satisfeito em desviar um pouco da leitura do livro Real World Haskell para aprender mais sobre combinadores. Graças a isso acabei descobrindo coisas interessantes como a origem do nome da empresa do Graham e que algumas funções muito usadas em Haskell são combinadores.

Além disso, embrora o combinador Y não seja tão útil na prática, as ideias por trás dele são muito interessantes (e, na minha opinião, complicadas!).

Referências

[1] Paul Graham – Biografia
[2] How to Mock a Mocking Bird – Smullyman
[3] StackOverflow – Y Combinator in Haskell
[4] [Haskell] How to define Y combinator in Haskell
[5] Combinator Birds
[6] Hackage – The data-aviary package
[7] Github – Homoiconic
[8] Mike’s World-O-Programming – The Y Combinator
[9] Y-Combinator em Javascript


Equações de Pell

Fevereiro 12, 2012

Equações diofantinas são equações polinomiais que só admitem soluções inteiras. Neste post, vou comentar sobre três casos particulares e me concentrar em um deles, conhecido como equações de Pell.

Triplas de Pitágoras

O caso mais conhecido de equações diofantinas é da forma:

x^2 + y^2 = z^2

As soluções inteiras para essa equação são conhecidas como triplas pitagóricas, por causa do teorema de Pitágoras, sobre a relação entre os comprimentos dos lados de um triângulo retângulo.

Para o caso mais geral, da forma x^n + y^n = z^n com n > 2, o famoso último teorema de Fermat diz que não existem soluções inteiras.

Indentidade de Bézout

Outro caso particular conhecido são as equações diofantinas lineares de duas variáveis, por exemplo, x e y, da forma

ax + by = c

Onde a, b, c são constantes inteiras. No caso em que c é o máximo divisor comum de a, b, temos a identidade de Bézout. É possível encontrar um par x, y inteiro usando o algoritmo de Euclides estendido.

Note entretanto que há infinitas soluções inteiras para essa equação, pois se x, y é uma solução, x' = x + b e y' = y - a é também uma solução.

Equação de Pell

Só recentemente ouvi falar sobre um outro caso particular, conhecido como equações de Pell, com a forma:

(1) x^2 - ny^2 = 1

Onde n é um inteiro positivo. Se n não possui raíz exata, então existem infinitas soluções inteiras x, y (Se n tiver raíz exata dá pra mostrar que a única solução é x = \pm 1 e y = 0). Vamos apresentar um algoritmo para encontrar as soluções para esse caso particular de n.

Frações continuadas

Um conceito usado no algoritmo é o de frações continuadas, que podem ser usadas para aproximar um número irracional através de frações. A fração continuada possui a seguinte forma:

Dado um número d, os coeficientes de sua fração continuada, a_0, a_1, \cdots, a_n são positivos e podem ser calculados através das recorrências, para i \ge 1:

r_i = \frac{1}{r_{i-1} - a_{i-1}}

a_i = \lfloor r_i \rfloor

e com a_0 = \lfloor d \rfloor.

Podemos representar uma fração continuada por seus coeficientes, ou seja, [a_0,a_1,a_2,\cdots]. Incrivelmente, no caso particular em que o número irracional é da forma \sqrt{n}, prova-se que os coeficientes da sua fração continuada são periódicos, ou seja \sqrt{n} = [a_0, \overline{a_1,a_2,\cdots,a_{r-1},a_{r}}] e a_{r} = 2a_0.

Como um exemplo, para n = 14 temos \sqrt{14} = [3,\overline{1,2,1,6}]

É possível calcular explicitamente o numerador p_i e o denominador q_i da fração continuada com i termos, através das recorrências:

\begin{array}{lcl}   p_0 & = & a_0\\  p_1 & = & a_1 a_0 + 1\\  p_n & = & a_n p_{n-1} + p_{n-2}   \end{array}

e

\begin{array}{lcl}   q_0 & = & 1\\  q_1 & = & a_1\\  q_n & = & a_n q_{n-1} + q_{n-2}   \end{array}

A fração p_i / q_i é chamada de n-ésimo convergente. Voltando ao exemplo para n = 14, temos:

\begin{array}{llcl}   i: & p_i / q_i  &  & \\  0:  & 3/1 & = & 3.0\\  1: & 4/1 & = & 4.0\\  2: & 11/3 & = & 3.66666666667\\   3: & 15/4 & = & 3.75\\  4: & 101/27 & = & 3.74074074074  \end{array}

Sendo que \sqrt{14} é aproximadamente 3.74165769645

Uma solução para a equação de Pell

Considere os coeficientes da fração continuada de \sqrt{n} e r o índice a partir do qual os coeficientes ficam periódicos.

Se r é par, seja x = p_{r-1} e y = q_{r-1}. Caso contrário, seja x = p_{2r-1} e y = q_{2r-1}. Surpreendentemente, existe um teorema que diz que x, y representam a menor solução inteira positiva para (1).

Algoritmo

Dadas essas propriedades, um algoritmo para encontrar uma solução para (1), consiste em iterar sobre os coeficientes p_i e q_i, até encontrar r tal que a_r = 2a_0. Segundo Lenstra [2], r \in O(\sqrt{n} \log n), ou seja, tem complexidade pseudo-polinomial. Assim, se d é o número de bits necessários para representar n, temos que r \in O(2^{d/2}d), ou seja, mesmo ignorando as complexidades das operações artiméticas no código, o nosso algoritmo é exponencial no número de bits.

Gerando mais soluções

Dada uma solução fundamental (x_1,y_1), Lenstra [2] afirma que se ordenarmos as soluções por magnitude, então a k-ésima solução (x_k, y_k) é tal que

(2) x_k + \sqrt{n}y_k = (x_1 + \sqrt{n}y_1)^k

Usando a ideia de [1], como ambos (x_1,y_1) e (x_k, y_k) são solução para (1), temos

x_k^2 - ny_k^2 = (x_1^2 - ny_1^2)^k = 1

Fatorando temos,

(x_k + \sqrt{n}y_k)(x_k - \sqrt{n}y_k) = (x_1 + \sqrt{n}y_1)^k(x_1 - \sqrt{n}y_1)^k

Por (2), concluímos que,

(3) x_k - \sqrt{n}y_k = (x_1 - \sqrt{n}y_1)^k

Resolvendo para (2) e (3), chegamos em:

\begin{array}{lcl}   x_k & = & \dfrac{(x_1 + y_1\sqrt{n})^k + (x_1 - y_1\sqrt{n})^k}{2}\\  \\  y_k & = & \dfrac{(x_1 + y_1\sqrt{n})^k - (x_1 - y_1\sqrt{n})^k}{2\sqrt{n}}  \end{array}

Implementação

Com base na teoria apresentada acima, é simples escrever um algoritmo. Adicionei uma implementação em python à minha biblioteca pessoal de teoria dos números, disponível no github.

Encontrei um post bastante interessante sobre equações de Pell e Haskell aqui [3]. Para praticar, resolvi implementar a minha versão antes de olhar o código do post. Apanhei bastante para lidar com o fato do Haskell não fazer conversão de tipos ponto flutuante e inteiro automaticamente. É preciso fazer isso explicitamente e por isso a função fromIntegral no código deste link.

Leitura complementar

[1] Wolfram – Pell Equation
[2] Solving the Pell Equation – H. W. Lenstra Jr.
[3] A Foray into Number Theory with Haskell


Topcoder: SRM 523

Novembro 27, 2011

No SRM 523 demorei bastante para fazer o problema fácil (50 min!) e por isso sobrou pouco tempo para fazer o médio. Cheguei a bolar uma solução, mas descobri no final que ela estava errada.

Depois de muito tempo voltei acertar um challenge! O erro em questão era um típico caso de Off-by-one error no problema fácil.

Depois de passada a competição, consegui encontrar o erro da minha ideia original e resolvi o problema médio. Neste post vou explicar minha solução.

BricksN

Considere que temos um número infinito de blocos de tamanho 1 a k. Queremos construir uma torre sobre uma base de largura w e altura máxima h usando esses blocos, com a seguinte restrição: só podemos posicionar um bloco em uma dada posição, se todas as posições logo abaixo desse bloco estiverem preenchidas (sem buraco). A pergunta é: quantas torres diferentes podemos formar seguindo essas restrições?

Exemplo de configuração proibida

(Observação: se duas torres têm mesma “forma”, mas existem blocos em posições diferentes, elas também são consideradas diferentes)

Soluções consideradas diferentes

Descrição do problema original.

Solução

Primeiramente, notei que este era um problema de contagem e portanto havia uma grande chance da solução ser por programação dinâmica.

Para encontrar a recorrência em programação dinâmica, é interessante pensar em casos particulares. Vamos por exemplo pensar em uma torre de altura 1 (uma linha) e largura w. Além disso, vamos supor que não há buracos. De quantas maneiras podemos formar tal torre usando blocos de tamanho até k?

Seja C_k(w), o número de maneiras de construir uma linha de comprimento w, permitindo blocos de tamanho até k. Considere todas as soluções com essas características. Podemos particionar essas soluções em grupos nos baseando no tamanho do último bloco utilizado.

Considere um grupo onde o tamanho do último bloco é i. Se removermos esse bloco, teremos linhas de comprimento (w-i). Todas as possíveis soluções diferentes com comprimento (w-i) continuarão diferentes ao adicionarmos o bloco de tamanho i de volta.

Fazendo isso para todos os blocos de tamanho 1 a k, obteremos todas as possíveis soluções contadas em C_k(w). Note que não estamos contando soluções repetidas pois duas soluções com o último bloco diferente, serão diferentes independentemente do que tem anteriormente. Note também que não estamos deixando de contar nenhuma solução, já que qualquer solução deverá conter como último bloco algum entre 1 e k.

Logo, temos a seguinte recorrência:

C_k(w) =  \sum_{i=1}^k C_k(w-i)

Para i \ge  w e com C_k(0) = 1.

Agora vamos generalizar um pouco mais e considerar linhas que podem ter buracos. A receita para encontrar uma recorrência é a mesma: precisamos particionar todas as possíveis soluções e para cada grupo, tentar reduzir para um subproblema.

Seja B_k(w) o número de maneiras de se formar uma linha de comprimento w, permitindo buracos.

Nesse caso, podemos particionar as soluções pela posição onde um buraco primeiro aparece, ou seja, temos um pedaço contíguo de blocos de tamanho i, seguido de pelo menos um buraco. Uma vez que estamos considerando um grupo de soluções para um dado i, qualquer solução diferente formada depois do buraco (posição i+2 em diante), continua sendo diferente ao adicionarmos o pedaço contínuo de comprimento i mais o buraco.

Também temos que contar os casos em que não há buracos. Dessa forma teremos a recorrência

B_k(w) = C_k(w) + \sum_{i=0}^{w-1} C_k(i) \cdot B_k(w-i-1)

Aqui também temos B_k(0) = 1 e podemos argumentar que soluções para i’s diferentes são diferentes (não estamos contando repetições) e também que todas as soluções estão sendo contempladas (inclusive a vazia).

Finalmente, vamos generalizar para torres de altura h. Seja A_k(w,h) o número de soluções respeitando as condições do problema original, sobre uma base w e altura h.

Observamos que devido às restrições do problema, não podemos fazer “pontes” com buracos em baixo. Assim, se tivermos um pedaço contínuo de comprimento w’ no primeiro nível de uma torre com altura máxima h, a construção que pode estar sobre ele são todas as torres sobre uma base w’, com altura máxima h-1.

Com essa observação, podemos continuar com nossa recorrência de B, mas agora temos que considerar todas as possíveis torres que vão sobre o primeiro bloco contínuo. Note que não precisamos nos preocupar com o restante dos blocos depois do buraco, pois qualquer diferente solução obtida continuará diferente após juntarmos com a nossa primeira torre.

Assim, ficamos com a recorrência:

A_k(w, h) = C_k(w) \cdot A_k(w, h-1) +
\; \sum_{i=0}^{w-1} C_k(i) \cdot A_k(i, h-1) \cdot A_k(w-i-1, h)

A implementação dessa ideia é direta. Coloquei no gitbhub essa solução.


Topcoder: SRM 517

Outubro 30, 2011

Já faz um tempo que ocorreu o SRM 517, mas só agora resolvi estudar a solução. No dia resolvi só problema mais fácil (250) e ainda errei um challenge :(

Neste post vou apresentar apenas a solução do problema médio, que estava mais difícil do que o normal, o que estava sugerido por sua pontuação (600 ao invés de 500).

Adjacent swaps

Descrição

Seja um vetor de N elementos v = {0, 1, 2, …, N-1} e uma permutação de p = {p[0], p[1], …, p[N-1]} de v. Definimos como troca-completa o conjunto das N-1 possíveis trocas entre elementos nas posições i e e i+1 (i variando de 0 a N-2) executadas exatamente uma vez, em uma ordem arbitrária.

O problema pede para encontrar quantas trocas-completas levam a v a p.

Solução

A primeira observação é que ao invés de considerar trocas-completas que levem v a p, consideramos aquelas que levam p a v. Existe uma bijeção entre uma troca-completa do problema original e do novo problema, que é simplesmente inverter a ordem das trocas e portanto a resposta é a mesma para os dois problemas.

A solução proposta no editorial é resolver esse problema através de programação dinâmica. Para tanto, definimos

u(i,j, k) — como o número de trocas-completas que ordenam o subvetor p[i…j] = {p[i], p[i+1] …, p[j]}, usando apenas as trocas entre elementos consecutivos entre i e j-1, com a restrição adicional de que a última troca efetuada seja entre k e k+1. (0 <= i <= k < j <= N-1).

Seja t(i,j) = \sum_{k = i}^{j-1} u(i,j,k). A solução que estamos procurando é então t(0, N-1). Além do mais, para o caso base i=j, temos um vetor de um único elemento ordenado e já não resta movimentos, portanto, t(i,i) = 1 para i=0, N-1.

Vamos ver como podemos calcular u(i,j,k).

Seja q[i…j] o subvetor p[i…j] ordenado e q'[i..j] o subvetor q[i..j] antes de fazer a troca entre k e k+1. Assim, q'[k] = q[k+1] e q'[k+1] = q[k], ou seja q'[i…j] = {q[i], q[i+1], …, q[k-1], q[k+1], q[k], q[k+2], …, q[j]}.

Observação 1: (q[i] < q[i+1] < … q[k-1] < q[k+1]) e (q[k] < q[k+2] < … < q[j]).

Por definição, q[i…j] e q'[i…j] são permutações de p[i..j]. Isso implica que todo elemento de q'[i…k] pertence a p[i…k] ou p[k+1…j]. Porém, como por hipótese a troca k,k+1 não foi feita ainda, um elemento de q'[i…k] não pode ter vindo de p[k+1…j]. Logo, concluímos que q'[i…k] deve ser uma permutação de p[i…k]. Temos assim duas condições a serem satisfeitas:

Condição 1: O subvetor p[i…k] deve ser uma permutação de q'[i…k].

Condição 2: O subvetor p[k+1…j] deve ser uma permutação de q'[k+1…j].

Na verdade podemos mostrar que as condições 1 e 2 são equivalentes. Supondo que as condições são satisfeitas, pela Observação 1, temos que q'[i…k] está ordenado, logo q'[i…k] = q[i…k], que é uma ordenação de p[i…k]. É possível argumentar o mesmo para q'[k+1…j] e p[k+1…j].

Por definição, t(i,k) conta o número de trocas-completas que leva p[i…k] em q[i…k] (o mesmo para t(k+1,j)). Assim, a princípio poderíamos combinar cada troca-completa em t(i,k) com cada troca-completa em t(k+1,j), levando à recorrência u(i, j, k) = t(i, k)*t(k+1, j).

Entretanto, há um número bem maior de combinações. Considere uma dada troca-completa A de t(i,k) e outra B de t(k+1, j). Embora a ordem das trocas em A e em B sejam fixas, não há restrição de ordem entre uma troca feita em A e outra feita em B. Por exemplo, poderíamos fazer todas as trocas de A primeiro e depois todas as de B ou qualquer combinação de intercalação entre os elementos de A e B.

Podemos mostrar que o número total de intercalações possíveis é dada por um coeficiente binomial C(n,r). No nosso caso em particular, n é o número total de trocas disponíveis, ou seja, j-i-1 (k-i de A e j-(k+1) de B) e r é o número de trocas em A (ou em B), ou seja, C(j-i-1, k-i). Para ver o porquê, considere que temos posições de 1 a j-i-1. Queremos escolher k-i dessas posições para atribuir a elementos de A e o restante das posições fica para B. Agora basta pensarmos nessas posições como a ordem em que as trocas são efetuadas.

Isso dá a recorrência

u(i, j, k) = t(i, k)*t(k+1, j)*C(j-i-1, k-i)

No caso de as condições 1 e 2 não serem satisfeitas, u(i, j, k) = 0.

Podemos usar a recorrência C(n,r) = C(n-1, r-1) + C(n-1, r) para pré-computar os valores dos coeficiente binomiais em O(N^2). Também podemos pré-computar se para uma dada tripla i, j, k, as condições 1 e 2 serão satisfeitas, em O(N^3) (usando ordenação linear para encontrar q).

Dado isso, podemos computar cada u(i, j, k) em O(1) e cada t(i, j) m O(N). A complexidade total das recorrências fica O(N^3), que também é a complexidade total do algoritmo.

Decidi colocar minhas implementações de SRM’s também no github.

Conclusão

Achei a solução desse problema bastante complicada e como de costume fiquei surpreso de tanta gente ter conseguido resolvê-lo durante os 75 minutos de competição.

Programação dinâmica é uma técnica que eu acho difícil aprender, porque em geral não consigo encontrar uma estrutura no problema que me indique a forma de atacá-lo. Minha esperança é continuar estudando suas aplicações e quem sabe aprender por “osmose” :P

O uso do coeficiente binomial para contar intercalações é bem engenhoso! Lembro que alguma outra vez estudei um problema em que a contagem por coeficiente binomial também não era trivial de ver.


Google Code Jam 2011

Setembro 18, 2011

O Google Code Jam é uma competição de algoritmos anual promovida pelo Google, que acontece desde 2007 se não me engano. Esse ano eu fui até o round 2, mas acabei não me classificando entre os 500 primeiros para avançar ao round 3:( Entretanto, acabei ganhando um prêmio de consolação por ter ficado entre os 1000 primeiros, que é a camiseta do evento.

Camiseta (frente)

Embora a competição tenha acabado em Julho, a camiseta só chegou a algumas semanas. O desenho na parte de trás são linhas de código em diversas linguagens de programação e tem o formato do Japão que é onde foi a final mundial, na qual o japonês Makoto Soejima (rng_58 no topcoder) sagrou-se vencedor.

Detalhe da parte de trás

Com a chegada da camiseta, lembrei que eu tinha ficado de estudar a solução dos problemas do round 2 que eu não havia conseguido resolver. Felizmente os problem setters apresentam a solução, mas também vou descrevê-las aqui do meu jeito.

C. Expensive Dinner

O enunciado desse problema está aqui. A solução está baseada em três proposições principais.

Proposição 1: Dado um grupo de amigos x_1, x_2, \dots, x_i, o valor da conta é sempre LCM(x_1, x_2, \dots, x_i) não importa a ordem em que eles cheguem. LCM representa o mínimo múltiplo comum.

Prova: Seja y o valor da conta.

Parte 1 – mostrar que y \ge LCM(x_1, x_2, \dots, x_i)
Essa parte é direta, dada a definição de mínimo múltiplo comum e que a conta tem que ser divisível por todo x_i.

Parte 2 – mostrar que y \le LCM(x_1, x_2, \dots, x_i)
Quando um amigo x_i pede alguma coisa, o valor da conta vai para o próximo múltiplo de x_i. Como LCM(x_1, x_2, \dots, x_i) é múltiplo de todos, o valor da conta não poderá “pular” esse valor e uma vez que esteja nesse valor, todos estarão satisfeitos e a conta não irá mais aumentar.

Juntando as partes 1 e 2, provamos a proposição.

Seja M o número de vezes que o garçon é chamado e M_{\max} e M_{\min} o maior e o menor valor possível de M, respectivamente. O problema pede M_{\max} - M_{\min}. Sejam p_1, p_2, \dots, p_P os P primos menores ou iguais a n (o número de amigos). Além do mais, seja e_i o maior inteiro tal que p_i^{e_i} \le n.

Proposição 2: M_{\max} = 1 + \sum_{i=1}^P{e_i}

Prova: vamos quebrar em duas partes novamente.

Parte 1 – mostrar que M_{\max} \ge 1 + \sum_{i=1}^P{e_i}
Para isso basta exibir um caso onde M = 1 + \sum_{i=1}^P{e_i}. Esse é o caso em que os amigos chegam na ordem 1, 2, \dots, n. Considere a chegada de x_i tal que x_i = p_i^k. É fácil ver que LCM(x_1, x_2, \dots, x_{i-1}) não é divisível por p_i^k e portanto o garçon terá de ser chamado. Quantos x_i existem satisfazendo x_i = p_i^k? Justamente \sum_{i=1}^P{e_i}. Além disso, quando x_1 = 1 chegar a primeira vez, o garçon também deverá ser chamado. Portanto, o número de vezes que o garçon foi chamado é pelo menos 1 + \sum_{i=1}^P{e_i}

Parte 2 – mostrar que M_{\max} \le 1 + \sum_{i=1}^{e_i}
Pela Proposição 1, o valor final da conta deverá ser LCM(x_1, x_2, \dots, x_n), que é igual a \prod_{i=1}^P p^{e_i}. Para qualquer ordem de amigos, temos LCM(x_1, x_2, \dots, x_i) = k \cdot LCM(x_1, x_2, \dots, x_{i-1}) para um inteiro k. É fácil ver que se um garçon for chamado, então k > 1, ou seja, estamos aumentando a potência que algum fator primo até chegarmos em \prod_{i=1}^P p^{e_i}, logo o garçon não pode ser chamado mais do que o número total de fatores na conta final, ou seja, \sum_{i=1}^{e_i} (+1 por causa da primeira chamada).

Juntando as partes 1 e 2, provamos a proposição.

Proposição 3: M_{\min} = P

Prova: vamos quebrar em duas partes novamente.

Parte 1 – M_{\min} \le P
Para isso basta exibir um caso onde M = P. Considere uma ordem onde os P primeiros amigos a entrarem são exatamente aqueles iguais a p_i^{e_i} para todo i = 1, \dots, P. É fácil ver que depois que eles entratem, o valor da conta será já LCM(x_1, x_2, \dots, x_n) e todo mundo que entrar depois já estará satisfeito, tendo o garçon sido chamado apenas P vezes.

Parte 2 – M_{\min} \ge P
Vamos por contradição: suponha que existe uma ordem de P-1 ou menos amigos tal que LCM(x_1, x_2, \dots, x_{P-1}) = LCM(x_1, x_2, \dots, x_n) = \prod_{i=1}^P p^{e_i}. Note que, para cada p_i, deve existir um amigo com fator p_i^{e_i} (caso contrário o LCM seria menor). Pelo princípio das casas dos pombos, como existem apenas P-1 amigos e P primos, há de haver algum amigo múltiplo de p_i^{e_i} e p_j^{e_j}. Entretanto, ambos p_i^{e_i} e p_j^{e_j} são maiores que \sqrt{n}. Isso implica que esse amigo tem valor maior que n, uma contradição.

Juntando as partes 1 e 2, provamos a proposição.

Dadas essas proposições, o problema se reduz a calcular (1 + \sum_{i=1}^P {e_i}) - P para todo primo P \le n. Considerando que podemos calcular e_i em O(\log(n)), levando a um algoritmo O(n \log(n)). Porém, como n pode ser até 10^{12}, temos que fazer melhor.

A sacada é reescrever a equação como 1 + \sum_{i=1}^P ({e_i} - 1). Isso quer dizer que só estamos interessados em fatores primos com potências maiores que 1. Logo, podemos nos restringir a primos menores ou iguais a \sqrt{n}, levando a um algoritmo O(\sqrt{n} \log(n)).

Depois dessa teoria toda, o algoritmo fica bem simples. Aqui está minha versão em C++, usando o crivo de eratóstenes para encontrar os primos.

D. A.I. War

Este problema se reduz para o seguinte problema: Dado um grafo G(V, E), encontrar o menor caminho entre s e t em um grafo não-direcionado e não-ponderado. Se houver várias soluções, escolha o caminho com maior vizinhança. Consideramos a vizinhança de um caminho como a união das vizinhanças dos vértices desse caminho, excluindo os vértices do próprio caminho. O problema pede para imprimir o comprimento desse caminho e o tamanho de sua vizinhança.

Vamos denotar por N o número de vértices nesse grafo e M o número de arestas.

Primeiramente encontramos as menores distâncias a partir de s, para todos os vértices, através de uma busca em largura. Denotamos por dist(v) a menor distância de s a v. O tamanho do menor caminho é dist(t).

Construímos então um grafo direcionado G' = (V, A), com mesmo conjunto de vértices do grafo de entrada e um arco (u, v) \in A se a dist(u) + 1 = dist(v). Podemos argumentar que todo caminho mínimo s-t no grafo original corresponde a um caminho direcionado de s a t nesse novo grafo.

Solução para o caso pequeno: podemos fazer uma força bruta a partir de s e, sempre que chegarmos em t, calculamos o tamanho da vizinhança desse caminho. Dessa forma estamos enumerando explicitamente todos os possíveis menores caminhos. Fiz uma estimativa e me pareceu que no pior caso, haverá cerca de O(\sqrt{N}^{\sqrt{N}}) caminhos. No caso fácil, P=36, o que dá uma ordem de 50k operações. No caso difícil, P=400, o que torna a abordagem inviável. De qualquer maneira, eis aqui minha solução para o caso fácil.

Solução para o caso geral: uma ideia é usar programação dinâmica para guardar o tamanho da vizinhança para um dado caminho, mas note que não basta representar um caminho por seu último vértice, pois ao adicionar um vértice a este caminho, sua contribuição para o tamanho da vizinhança pode depender de mais vértices nesse caminho.

A observação chave é que na verdade, a contribuição de um novo vértice c só depende do último e penúltimo vértice desse caminho. Para entender porque, suponha que dist(c) = D + 1 e sejam b e a o último e penúltimo vértices de um caminho, respectivamente. Então dist(b) = D e dist(a) = D-1. Se c possuísse vizinhança comum com algum vértice x tal que dist(x) \le D-2, então o menor caminho até c seria menor ou igual a D, uma contradição.

Logo, é suficiente representar um caminho por seus dois últimos vértices. Seja f(b, c) o tamanho da maior vizinhança de um caminho terminado em b \rightarrow c. Podemos definir a seguinte recorrência:

f(b, c) = \max \{ (a, b) \in A(G') | f(a, b) + g(c, a, b) \}

Aqui, g(c, a, b) é o número de vértices na vizinhança de c que não está contido nem na vizinhança de a nem na de b. O caso base é f(s, c), onde o valor é simplesmente o tamanho da união da vizinhança de s e c.

Se pré-computarmos g(c, a, b), podemos resolver a recorrência acima em O(N) e como cada par (b, c) deve pertencer a A(G'), o algoritmo tem complexidade O(NM).

Porém, calcular g(c, a, b) é a parte mais demorada. A solução mais simples seria enumerar todas as triplas de vértices e processar cada vértice v vizinhança de c, levando a um algoritmo O(N^4), o que dá na ordem de 25 \times 10^9, que é muito lento.

A sacada é usar o fato que c e v devem ser adjacentes, assim como os vértices a e b. Assim, para cada par de arestas e_1 e e_2, supomos e_1 = (c, v) e e_2 = (a, b). Usando uma matriz de adjacência, podemos verificar se v pertence à vizinhança de a e b em O(1), levando a uma complexidade total de O(M^2) ~4M operações. O editorial do Google Code Jam apresenta essa e outras 3 alternativas para pré-computar g(c, a, b). Minha implementação da solução apresentada se encontra aqui.


Algoritmo de Karatsuba

Julho 3, 2011

Depois de bastante tempo, decidi voltar ao problema Scott New Trick, do Facebook Hacker Challenge desse ano. Eu analisei uma solução em um post anterior, no qual também disse que ia estudar a solução usando o algoritmo de Karatsuba.

Entretanto, descobri que o algoritmo de Karatsuba é usado apenas para realizar a multiplicação de polinômios rapidamente, em O(nˆ(log_2(3))) ou aproximadamente O(nˆ1.58). A multiplicação ingênua leva O(nˆ2).

Primeiro vou explicar a solução e depois descrever o algoritmo. Baseei minha análise no código do tomek.

Para relembrar o que o problema pede, cito um parágrafo do post anterior:

Dados dois vetores A e B de inteiros aleatórios, com tamanhos N e M, respectivamente, um primo P e um inteiro 1 ≤ L ≤ P contar quantos índices 1 ≤ i ≤ N e 1 ≤ j ≤ M satisfazem (A[i] * B[j]) mod P < L.

Lembrando que N e M são da ordem de 10,000,000 e P é da ordem de 250,000.

Raíz primitiva módulo P

Dado um inteiro P, dizemos que G é uma raíz primitiva módulo P, se para todo coprimo i de P (i.e. i < P tal que mdc(i, P)=1), existe uma potência de G módulo P que é igual a i, ou seja, existe um inteiro k, tal que Gˆk % P = i, onde k é dito logaritmo discreto de i.

Em particular, se P é primo e G é uma raíz primitiva módulo P, para todo número i de 1 a P-1 existe uma potência de G módulo P igual a i.

Uma maneira de determinar a raíz primitiva módulo P é através do seguinte teorema: um número G é raíz primitiva módulo P se a ordem multiplicativa de G módulo P é .

Em outras palavras, isso quer dizer que se o menor k positivo tal que Gˆk % P = 1, é k = , então G é uma raíz primitiva módulo P. Lembrando que para P primo, o código abaixo determina uma raíz primitiva G, além de calcular gpow[k] = Gˆk % P e gpowinv[k] = x, onde x é tal que Gˆx % P = k.

vector<int> gpow, gpowinv;
void CalcGen(int P){
    gpow.resize(P);
    while(1){
        long long G = rand() % (P - 1) + 1;
        gpow[0] = 1;
        bool ok = true;
        for(int i=1; i<P-1; i++){
            gpow[i] = long long(gpow[i-1])*G % P;
            if(gpow[i] == 1){
                ok = false;
                break;
            }
        }
        if(ok) break; 
    }
    gpowinv.resize(P);
    for(int i=0; i<P; i++) gpowinv[gpow[i]] = i;
}

Esse algoritmo probabilístico fica chutando candidatos a raízes primitivas até encontrar algum que tenha ordem multiplicativa igual a P-1. Não sei ao certo qual a sua complexidade, mas deve ser algo em torno de O(P√P) igual ao algoritmo do post anterior.

Representação por polinômio

Antes de mais nada, vamos supor que os elementos de A e B são menores do que P. Se não for, basta fazer A[i] = A[i] % P que a contagem não mudará.

Vimos no post anterior que podemos guardar os elementos de A e de B em um vetor de frequências de tamanho P de acordo com seu módulo P, como por exemplo fazendo fA[A[i]]++ para todo elemento A[i] de A. Entretanto vimos também que isso leva a um algoritmo O(Pˆ2), o que não passa no tempo.

A sacada é representar A[i] por Gˆx % P, onde x = gpowinv[A[i]] e B[j] por Gˆy % P, onde y = gpowinv[B[j]]. Note que A[i]*B[j] é equivalente a Gˆ(x+y) % P.

Que tal se fizermos nosso vetor de frequências baseado em x e y, por exemplo fA[x]++ e fB[y]++? Então, fA[x]*fB[y] conta quantos pares de elementos A[i] e B[j] existem em A e B, tal que A[i]*B[j] = Gˆ(x+y) % P = gpow[x + y]. Então, se gpow[x + y] < L, existem fA[x]*fB[y] pares que devem ser contados.

Se considerarmos fA e fB como um polinônimo, sendo fA[x] o coeficiente do termo de grau x, e fAB sendo o resultado de fA*fB, então fAB[z] é a soma de todos os coeficientes fA[x]*fB[y] tais que x+y=z, ou seja, fAB[z] conta todos os pares A[i] e B[j] tais que A[i]*B[j] = Gˆz % P.

O problema se reduz a contar a quantidade fAB[k] para todo k tal que gpow[k % (P-1)] < L. Observe que o módulo deve ser P-1 pois o período de Gˆk % P é P-1, ou seja, gpow[k] = gpow[k + P-1].

Algoritmo de Karatsuba

O algoritmo de Karatsuba é baseado em divisão e conquista. Dado um polinômio na forma:

Dado m = n/2, podemos representá-lo como:


ou de forma compacta por . Multiplicar dois polinômios x e y é então dado por:

Chamando de

(i)
(ii)

é possível mostrar a seguinte relação:

(iii)

Podemos multiplicar recursivamente para calcular (i), (ii) e a primeira multiplicação de (iii). A complexidade da recorrência é 3O(n/2) + cn, o que no total, usando o teorema mestre, é .

A multiplicação resultante é dada por:

Segue uma implementação em C++. Ela supõe que o vetor de entrada tem o tamanho que é uma potência de 2. Isso não é uma restrição muito forte, pois podemos completar o polinômio com zeros até que o tamanho satisfaça a restrição.

/**  
 * Karatsub algorithm. Fast multiplication of the first n
 * positions of fa and fb.
 * PS. Assumes n is power of 2 
 **/
void Karatsuba(int n, i64 *fa, i64 *fb, Poly &fab)
{
    int m = n/2;
    /* Base */ 
    if(n <= 16){
        fab.assign(2*n - 1, 0);
        for(int i=0; i<n; i++)
            for(int j=0; j<n; j++)
                fab[i+j] += fa[i]*fb[j];
        return;
    }
    Poly z0, z1, z2, tmp;
    /* Find z0 and z2 recursively */
    Karatsuba(m, &fa[0], &fb[0], z0);
    Karatsuba(m, &fa[m], &fb[m], z2);
    /* Find z1 recursively */
    Poly fai(m), fbi(m);
    for(int i=0; i<m; i++){
        fai[i] = fa[i] + fa[i+m];
        fbi[i] = fb[i] + fb[i+m];
    }
    Karatsuba(m, &fai[0], &fbi[0], z1);
    for(int i=0; i<z1.size(); i++)
        z1[i] -= (z0[i] + z2[i]);
    /* Merge z0, z1 and z2 in fab */
    fab.assign(2*m + z2.size(), 0);
    for(int i=0; i<z0.size(); i++){
        fab[i] += z0[i];
        fab[i+m] += z1[i];
        fab[i+2*m] += z2[i];
    }
}
/* Example of calling */
Karatsuba(fa.size(), &fa[0], &fb[0], fab);

O código completo encontra-se aqui.

Conclusão

No final das contas o algoritmo de Karatsuba era apenas uma maneira de resolver um problema conhecido (multiplicação de polinômios) de maneira rápida.

Achei a solução muito interessante de se transformar o domínio do problema para que ele possa ser resolvido de maneira mais eficiente. Isso me lembra a transformada de Fourier, mas como não sei muito sobre o assunto, posso estar enganado.