Algoritmo de Karatsuba

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.