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

Anúncios

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.


Facebook Hacker Cup – Round 2 (Continuação…)

Março 27, 2011

A final do facebook hackercup aconteceu dia 11 de Março em Palo Alto. O campeão foi o Petr Mitrichev.

Eu, como mero mortal (:P), ainda estou tentando entender as soluções dos problemas do round 2. Nesse post vou comentar a solução do problema Scott’s New Trick.

A descrição do problema é bem simples: 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.

Solução:

Bom, é trivial resolver esse problema em O(NM), enumerando todos os índices. O problema é que 2 ≤ N, M ≤ 10.000.000, então temos que pensar em uma solução melhor.

Uma segunda ideia é criar vetores de frequências fA e fB, ambos de tamanho P. A entrada fA[i] guarda quantos elementos v de A existem tal que v % P = i. O vetor fB é definido analogamente para B. Em O(P^2) percorremos todos os i e j. Se i * j % P < L, somamos fA[i]*fB[j] na conta final. Temos que P < 250.000, mas ainda é inviável rodar 20 testes em menos 5 minutos (que era o tempo para submeter as respostas).

Minha ideia durante a competição foi tentar paralelizar esse programa usando OpenMP em uma máquina com 4 cores. Consegui reduzir o tempo de execução para 1 minuto no pior caso. Porém, eu ia precisar de umas 5 ou mais máquinas dessas para conseguir submeter a tempo :( Ouvi dizer que alguém passou esse problema com esse algoritmo paralelizando em um cluster!

Parece que a maior parte das pessoas que passaram o problema, implementou o algoritmo de Karatsuba (o autor do algoritmo é russo e não japonês como o nome sugere, hehe). Vou estudar essa solução posteriormente, mas agora quero apresentar a solução bastante inteligente do meret, que melhora o desempenho do código usando teoria dos números e programação dinâmica.

Vamos continuar com os vetores fA e fB. Notamos que a soma que estamos buscando é dada por

(1)

O primeiro somatório conta os elementos onde os elementos de A são 0 (% P). Logo a multiplicação deles por quaisquer elementos de B vai ser sempre 0, e menor do que L, já que L ≥ 1.

Para entender o segundo somatório, considere um dado j. O valor inv(j) é o inverso modular de j, ou seja, é um número x tal que j * x ≡ 1 (mod P). Ele conta os números com produto igual a j*(k*inv(j)) % P, como multiplicação modular é comutativa e associativa, esse produto é simplesmente k. Por isso k vai até L-1.

Note que a restrição de P a primo é importante para que haja exatamente um único inverso modular para todo j de 0 a P-1. Além disso, todas as contas de índice são módulo P.

Agora, considere inteiros ‘c’ e ‘i’ tais que c ≡ i*j (mod P), equivalentemente, c * inv(j) ≡ i (mod P). Dados dois inteiros x e y, podemos escrever x como

(2)

Pois é o quociente e x % y é o resto da divisão. Sendo o índice k (da equação (1)) e c ambos inteiros, podemos escrevê-los na forma da equação (2). Além disso, multiplicando ela por inv(j) ficamos com

Fazendo módulo P e substituindo c * inv(j), temos

(3)

Na equação (1) podemos colocar para fora do somatório de k o termo fA[j]. Fixando um dado j, temos

Agora vamos substituir cada fB[k * inv(j)] pelo equivalente em (3). Para k = 0 até c-1, , então os índices serão

.

Para k = c até 2c-1,

e assim por diante, até que , onde os índices serão


Vamos reorganizar esse somatório agrupando os índices de forma que (k % c) sejam iguais, ou seja, soma os termos da primeira coluna, depois os termos da segunda coluna e por aí vai. O primeiro grupo fica,

(4)

O segundo grupo fica (5)

O k-ésimo grupo fica (6)


Note que são c grupos dessa forma. Introduzimos diversas variáveis, rearranjamos os termos, mas o somatório continua o mesmo. O intuito dessas complicações é permitir o cálculo rápido dos grupos acima. Para isso, dado um i, considere o vetor V, onde V[k] = i*k % P para k = 0, …, P-1.

A sacada da solução é notar que todo grupo de índices definido acima, corresponde a um intervalo de V. É fácil ver que o primeiro grupo de índices (eq. (4)), corresponde ao intervalo 0 a de V.

O segundo grupo é mais difícil de enxergar. Será que existe um intervalo de V correspondente a (5)? O primeiro termo é inv(j). Queremos encontrar um k’ tal que k’i ≡ inv(j) (mod P). Para tal, multiplique essa igualdade por j, ficando com k’*i*j ≡ 1 (mod P). Substituindo i*j por c, temos que k’c ≡ 1 (mod P) ou seja, k’ é o inverso modular de c!

Agora os índices do segundo grupo aumentam de i em i, então é fácil ver que eles formam um intervalo em V (possivelmente dando a “volta” no vetor). Como são termos, dá pra descobrir o intervalo de V, que será de inv(c) a inv(c) + .

Para um k genérico, é possível mostrar que k’ = k * inv(c). O intervalo fica de k*inv(c) a k*inv(c) + .

Vamos pré-computar um vetor com a soma acumulada dos elementos de fB indexados por V, denotado S, ou seja

Para todo n de 0 a P-1. Esse vetor pode ser calculado em O(P). Para cada grupo de índices, correspondente ao intervalo (lb, ub) em V, calculamos a soma

Dado um j, a complexidade fica O(P) para construir S e O(c) para computar as somas dos grupos de índices.

A primeira abordagem é, para cada j, escolher um i que minimize c. Claramente o i que estamos procurando é inv(j), pois nesse caso c=1. Aí podemos calcular as somas dos grupos em O(1), mas a complexidade de calcular S é O(P). Isso torna o algoritmo O(P^2), o que é tão ruim quanto o algoritmo mais simples.

No outro extremo, como S só depende de i, podemos escolher um i fixo para todo j e calcular c em função de i e j. Aí calculamos S apenas uma vez em O(P) e aí para cada j calculamos a soma em O(c). A complexidade fica O(P + P*C), onde C é o maior valor de c considerando todos os j’s. O problema é que C pode ser da ordem de P.

A ideia é usar uma abordagem híbrida. Escolhemos um i que minimize c, mas limitamos o valor de i a √P. Calculamos então o vetor imin, onde para cada j, imin[j] contém i (1 a √P) que minimize c. Aí então calculamos S apenas √P vezes. A complexidade fica O(P√P + P*C). Porém C ainda pode ser da ordem de P. Podemos apertar a análise da complexidade como O(P√P + ∑(c)), onde ∑(c) é a soma dos c’s para cada j.

Aí entra uma conjectura/teorema de que ∑(c) é O(P√P). (No final das contas isso estraga a beleza da solução. Seria legal ver a prova disso…). Com tudo isso, temos a complexidade O(P√P), que é suficiente para passar no tempo.

Conclusão

Apesar da conjectura, a solução do meret tem várias ideias interessantes. Vale a pena ver a implementação dele (dá pra baixar no site do Facebook). Não fiquei muito empolgado para implementar essa solução, mas quero ainda estudar as soluções usando o algoritmo de Karatsuba. Fica para um próximo post.


Facebook Hacker Cup – Round 2

Março 11, 2011

Recentemente teve uma competição promovida pelo Facebook, chamada Facebook Hacker Cup. Por mais que o nome sugira outra coisa, essa competição é nos moldes da Maratona de Programação e muito similar ao Google Code Jam.

Teve uma fase de qualificação, onde bastava fazer pelo menos um problema. Depois teve o Round 1 que era dividido em três partes: A, B e C. Você podia tentar qualificar para o Round 2 em qualquer uma delas desde que não tivesse conseguido em uma anterior.

Essas duas primeiras fases foram relativamente tranquilas, mas o Round 2 foi muito difícil! Os 300 primeiros dos ~3000 classificados iriam ganhar uma camiseta, mas só 150 pessoas fizeram pelo menos um problema e eu não estava entre elas :/

Os 25 primeiros colocados vão competir na final onsite lá em Palo Alto. Até onde eu sei não há brasileiros.

Vou apresentar apenas o problema mais fácil dessa etapa. Os outros eu ainda não sei como fazer :P Se eu descobrir faço um novo post.

Bonus Assignments

O descrição do problema se encontra aqui. No problema você deve distribuir uma certa quantidade de dinheiro para N trabalhadores para eles comprarem refrigerante, com a restrição de que o trabalhador que ganhar menos não deve ganhar menos do que A ou mais do que B. O trabalhador que ganhar mais não deve ganhar menos do que C ou mais do que D. Uma restrição adicional é que pelo menos um trabalhador não vai conseguir torrar todo o seu dinheiro em refrigerante. Você não sabe o preço exato do refri, só que é mais do que $1. Queremos saber quantos jeitos temos de distribuir dinheiro, satisfazendo essas restrições. A resposta deve ser módulo 1,000,000,007.

A primeira observação é que, para que sempre sobre algum trabalhador que não consiga torrar todo o seu dinheiro com refri, independente do custo desse, deve existir pelo menos dois trabalhadores que recebam montantes x e y primos entre si. Caso contrário, se todo par tivesse montante com algum divisor comum mdc > 1, então para um refri com esse valor todos conseguiriam torrar seu dinheiro.

O problema se reduziu a contar quantas combinações existem para N elementos, de modo que o valor do menor elemento esteja entre [A, B] e do maior entre [C, D] e exista pelo menos um par de elementos com valores primos entre si.

Solução:

Vamos por partes. Primeiro, como contar quantas combinações existem de N elementos, sendo que o valor deles deve estar entre [A, D]? Por exemplo, N = 3, A = 1, D = 5. Cada elemento pode assumir o valor {1, 2, 3, 4, 5}. Como eles são independentes, há 5^3 combinações. De maneira mais geral, temos (D-A+1)^N.

Porém, queremos que o menor elemento tenha valor de no máximo B. Por exemplo se B = 2, na conta acima estamos contando o vetor [3, 3, 5], cujo menor elemento é 3. Então queremos jogar fora os vetores cujo menor elemento sejam maior do que B. Esse número é justamente quantas combinações de N elementos temos com valor entre [B+1, D]. Ou seja, devemos descontar (D – B)^N.

De igual maneira, ainda estamos contando vetores cujo maior elemento seja menor do que C. Por exemplo, se C=4, temos [1, 1, 1] na conta acima. Vamos jogar fora aqueles vetores cujo maior elemento seja menor do que C, ou seja, aqueles com N elementos de valores entre [A, C-1], dado por (C – A)^N.

Ficamos com um problema. Observe que o vetor [3, 3, 3] é um vetor com valores entre [A=1, C-1=3] e [B+1=3, D=5]. Ou seja, removemos ele duas vezes! De modo geral, removemos todos os vetores que têm valores entre [B+1, C-1] duas vezes. Devemos repô-lo com (C – B – 1)^N. Observe que essa interseção só ocorre se C > B+1.

Bom, sem a restrição de ter pelo menos dois valores primos relativos a conta é basicamente:

(D-A+1)^N – (D – B)^N – (C – A)^N + (C – B – 1)^N,

onde o último termo é zero se C <= B+1.

Nessa contagem estamos contando tanto os vetores com pelo menos dois primos relativos, quanto aqueles em que nenhum par é primo relativo, ou seja, os vetores para os quais todos os valores tem um divisor comum maior do que 1. Vamos tentar encontrar aqueles com mdc=2. Por exemplo, para N = 3, A = 3, B = 4, C = 4, D = 7. Um exemplo é [4, 6, 4]. Claramente podemos dividir esse vetor por 2, resultando em [2, 4, 2]. Fica claro que podemos dividir todos os limites por 2 e obter um subproblema que é justamente a quantidade de combinações onde todos os valores são pares.

Assim, dividimos A, B, C e D por 2 para obter A’=2, B’=2, C’=2, D’=3. Note que usamos teto para os limites inferiores A e C, mas chão para os superiores B e D, pois não queremos considerar um intervalo maior do que o original. Note também que ao dividir por 2 e contar, estamos considerando os casos com mdc=4, mdc=8, etc. Ou seja, só precisamos remover as divisões por primos!

Bom, basta então dividir por todos os primos para os quais os intervalos resultantes tenham algum vetor válido. Um limite suficiente é tentar dividir por todos os primos menores ou iguais que D. Porém, temos novamente um problema de remoção repetida. Ao contar os que têm 2 como mdc e os que têm 3 como mdc, estamos contando duas vezes aqueles que têm 6 como mdc! Por isso precisamos usar o princípio da inclusão-exclusão. Dá para fazer isso recursivamente.

Acho que essa teoria é suficiente para resolver o problema. Tem o detalhe de achar primos de forma eficiente, que pode ser feita com o crivo de Eratóstenes e fazer a exponenciação modular de forma rápida.

Tem uma pegadinha adicional, que eu só descobri depois que implementei: nada impede que B > D ou C < A. Algumas das contas podem dar errado por causa disso, mas note que o resultado não muda se fizermos B = min(B, D) e C = max(C, A) nesse caso.

Minha implementação em C++ está aqui. Baseei minha análise e código na solução do kalinov.