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.

Anúncios

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.