Mapas de símbolos proporcionais e a meia integralidade da cobertura de vértices

Julho 31, 2011

O problema que ataquei no mestrado consistia em achar uma ordem de empilhamento entre símbolos em um mapa, de modo a maximizar uma dada função objetivo. Não se sabe a complexidade desse problema e por isso usei programação linear inteira.

Há uma variante desse problema, que consiste em maximizar uma outra função objetivo, para a qual existe um algoritmo polinomial.

Depois de obter as soluções para o meu problema, precisei fazer uma comparação entre as soluções desses dois problemas. Aconteceu que a diferença visual entre as duas soluções não é tão perceptível, sendo que apenas alguns discos estão em posição diferente. Para se ter uma idea, veja um dos casos


Jogo dos 7 erros ;) (clique na imagem para aumentar)

Decidi fazer um programa para destacar apenas os discos que mudaram de lugar. A minha primeira ideia foi marcar todos os discos que estavam em níveis diferentes entre as duas soluções.

Isso não funciona bem quando temos duas componentes disjuntas de discos, por exemplo na figura abaixo, com duas componentes disjuntas {A, B, C} e {D, E}. Note que as ordens {A > B > C > D > E} e {D > E > A > B > C} resultam visualmente na mesma figura, embora nenhum disco esteja exatamente na mesma posição na ordenação e portanto o programa vai marcar todos os discos.

Desenho que pode ser representado por mais de uma ordem dos discos

Uma segunda ideia foi marcar, considerando todos os pares de discos que se interceptam, apenas aqueles em que a ordem relativa do par foi invertida. Com essa estratégia, nenhum disco do exemplo acima seria destacado, o que está de acordo com nossas expectativas.

Porém, essa abordagem ainda apresenta um efeito indesejável. Se tivermos uma pilha de discos, onde todos se interceptam e um disco que está no topo em uma solução aparece na base na outra, sendo que todos os outros discos permaneceram intactos, como por exemplo {A > B > C > D} e {B > C > D > A}, ainda assim todos os discos serão destacados pois todos os pares formados com A tiveram sua ordem relativa invertida. Note que seria suficientemente informativo destacar apenas o disco A.

Assim, para todo par de discos que se intercepta, basta destacarmos um deles. Como então destacar o menor número possível para diminuir a poluição visual?

Se olharmos para o grafo de discos desse problema, existem algumas arestas que queremos destacar (cobrir), usando pelo menos um dos vértices mas tentando minimizar a quantidade desses vértices. Soa familiar?

Sim, reduzimos para o problema da cobertura de vértices (vertex cover).

Cobertura de vértices

O problema da cobertura de (arestas por) vértices é sabidamente NP-difícil. Assim, é pouco provável um algoritmo polinomial para resolvê-lo de maneira ótima. No meu trabalho, usei uma heurística simples: percorra os vértices começando por aqueles de maior grau. Se todas as arestas incidentes a esse vértice já foram cobertas por algum vértice da cobertura, ignore-o. Caso contrário adicione-o à cobertura.

Aplicando esse algoritmo nas figuras do começo do post, ficamos com o seguinte resultado:

Ficou mais fácil encontrar as diferenças entre as soluções (clique na imagem para aumentar)

Esse método ainda tem um problema que é quando um par de discos tem sua ordem relativa trocada, mas a parte que foi visualmente modificada está encoberta por um outro disco. Um exemplo disso pode ser visto no canto inferior direito da figura acima. Para consertar esse problema teríamos que fazer um pré-processamento e só considerar vértices cujos discos têm alguma parte do seu borda diferente entre uma solução e outra.

Modelo de programação linear

Como não poderia deixar de ser, apresento aqui o modelo de programação linear inteira da versão mais geral desse problema, na qual os vértices têm um custo associado. Dado um grafo G = (V, E) com uma função de custo c : V \rightarrow \mathbb{R}^+, definimos a variável binária x_v para cada v \in V, que vale 1 se o vértice v está na cobertura e 0 caso contrário. O modelo então fica:

Sujeito a

(1) Para todo (u, v) \in E,

(2) Para todo v \in V,

Meia integralidade

O modelo apresentado tem uma propriedade interessante chamada de meia-integralidade. Isso quer dizer que, dada a relaxação linear desse modelo, existe uma solução ótima onde os valores das variáveis x_v são 0, 1/2 ou 1. A prova é relativamente curta e está descrita na Seção 14.3 de [1].

Assim, dada uma solução fracionária do modelo, podemos usar a técnica de arredondamento, criando uma solução x'_v com x'_v = 1 sempre que x_v \ge 1/2. É possível mostrar que essa solução é uma cobertura de vértices válida e que seu custo não é maior do que 2 vezes o valor da solução ótima, resultando em um algoritmo 2-aproximado.

Referências

[1] Approximation Algorithms – Vijay V. Vazirani (2003)

Anúncios

Shader OSL Multithread

Julho 24, 2011

Essas últimas semanas eu estava tentando implementar suporte a múltiplas threads do shader OSL. O BRL-CAD usa pthreads para paralelizar o trabalho.

A primeira vez que eu pus o shader OSL pra rodar, tive problemas, provavelmente devido à leitura e escritas simultâneas a dados na memória. Eu tinha adiado esse problema e estava trabalhando apenas com uma thread.

O que a minha função osl_render faz é preencher uma estrutura chamada RenderInfo com dados como ponto de interseção, a normal, um ponteiro para o shader dessa superfície, etc. Preenchida essa estrutura, ela faz uma chamada ao método QueryColor da classe OSLRender. O objeto dessa classe é global, o que significa que está em uma posição de memória compartilhada pelas threads. Como a chamada do método QueryColor modifica o objeto, temos o problema do acesso simultâneo a esse método.

Primeira tentativa

A primeira coisa que eu tentei, foi proteger a chamada dessa função com semáforos. O BRL-CAD tem uma biblioteca interna de utilitários, entre eles a implementação de um semáforo.

Basicamente, o que eu fiz foi o seguinte:

OSLRender oslr;
...
bu_semaphore_acquire(BU_SEM_SYSCALL);
oslr.QueryColor();
bu_semaphore_release(BU_SEM_SYSCALL);

BU_SEM_SYSCALL na verdade é uma macro representando um inteiro. Quando uma thread chama a função bu semaphore acquire(BU_SEM_SYSCALL), ela bloqueia o acesso a outras threads para o mesmo inteiro.

Com isso o código passou a funcionar para mais de uma thread! Porém, fui conservador demais ao bloquear a chamada do QueryColor. Creio que 90% do trabalho da função osl_render seja feito lá. Logo, colocando um semáforo nessa função, a paralelização vai por água abaixo. De fato, fiz testes com 1 e 2 processadores e os tempos para renderizar uma imagem de teste foram exatamente os mesmos.

Segunda tentativa

Fui estudar um pouco mais o código do OSL e descobri que ele tem suporte a aplicativos multithread. Basicamente, cada thread deve manter uma estrutura local chamada ThreadInfo. Ao chamar o sistema de shaders do OSL, passamos essa informação, que será modificada.

Gambiarra!

Para manter uma ThreadInfo para cada thread, precisamos identificá-las de alguma forma. O problema é que não temos acsso a essa informação diretamente. O único jeito que encontrei sem ter que mudar a estrutura interna do raytracer do BRL-CAD foi a seguinte: toda vez que a função osl_render é chamada, uma estrutura chamada struct application é passada. Essa estrutura tem um campo chamado resource, que é um ponteiro para um pedaço de memória exclusivo para cada thread.

Esse pedaço é alocado na criação das threads e não muda ao longo de suas vidas. Assim, usar o endereço desse pedaço de memória é um jeito (feio) de identificá-las. Outro problema é que, embora haja uma função de inicialização do shader (chamado osl_setup), as threads são criadas após isso, o que me obriga a inicializá-las na função osl_render mesmo.

A solução na qual cheguei foi a seguinte:

/* Every time a thread reaches osl_render for the 
   first time, we save the address of their own
   buffers, which is an ugly way to identify them */
std::vector<struct resource *> visited_addrs;
/* Holds information about the context necessary to
   correctly execute a shader */
std::vector<void *> thread_infos;

int osl_render(struct application *ap, ...){
    ...
    bu_semaphore_acquire(BU_SEM_SYSCALL);
    
    /* Check if it is the first time this thread is
       calling this function */
    bool visited = false;
    for(size_t i = 0; i < visited_addrs.size(); i++){
        if(ap->a_resource == visited_addrs[i]){
            visited = true;
            thread_info = thread_infos[i];
            break;
        }
    } 
    if(!visited){
        visited_addrs.push_back(ap->a_resource);
        /* Get thread specific information from 
           OSLRender system */
        thread_info = oslr->CreateThreadInfo();
        thread_infos.push_back(thread_info);
    }
    bu_semaphore_release(BU_SEM_SYSCALL);
    ...
}

Classe thread-safe

Com as ThreadInfo's, podemos nos despreocupar com o acesso concorrente ao objeto da classe OSLRender. Para ter uma garantia de que não será feita escrita nesse objeto, podemos declarar os métodos como const, como no exemplo abaixo:

struct A {
    int x;
    // OK
    void set_x(int _x) { x = _x; };
    // ERRO, tentando modificar membro do objeto 
    // num método const
    void const_set_x(int _x) const { x = _x; };
};

Testes

Com essa nova implementação, os tempos melhoraram bastante com mais processadores. A tabela abaixo mostra os tempos de execução para renderizar uma cena de exemplo, usando até 4 processadores.

Gráfico 1: Tempo x no. de processadores

Os ganhos foram praticamente lineares, sendo que para 4 processadores a medida de paralelização foi de 1.25 (se fosse perfeitamente linear seria 1.0).

Próximos passsos

Falei com meu mentor sobre o desenvolvimento daquele modo de visualização incremental, mas parece que ele não gostou muito da ideia. Isso exigiria uma mudança bastante grande no sistema, pois o renderizador do BRL-CAD foi desenvolvido para ray-tracing e não path-tracing.

Por enquanto estou estudando maneiras de adaptar o código OSL para suportar ray-tracing, mas não sei se isso é viável. A propósito, eu ainda confundo bastante os termos ray-tracing, path-tracing, photon mapping e todos esses algoritmos de iluminação e pretendo em breve escrever um post cobrindo esses tópicos bem por cima.


GDB Gnu Debugger

Julho 17, 2011

Acho que fui apresentado ao GDB logo que aprendi a programar, mas nunca me acostumei a usá-lo. Talvez por causa da maratona, sempre fui adepto do debug com printf, exceto quando há falha de segmentação. Mas mesmo nesse caso eu acho o valgrind mais eficaz.

Entretanto, dizem que essa é uma má prática e que o GDB pode ser uma ferramenta essencial na hora de encontrar causas de bugs. Pessoalmente, não me lembro de nenhuma vez em que o GDB me ajudou a achar um bug que eu não teria achado de outra forma, tão ou mais rapidamente quanto.

Mas isso deve ser porque eu não aprendi a usá-lo direito :P Em uma nova tentativa de melhorar meu uso com essa ferramenta, decidi organizar informações em um post de forma que eu possa vir a aprender novas funcionalidades.

Exemplo Básico

Para poder rodar código com GDB, devemos passar a flag -g ou -ggdb.

Dado um exemplo teste.c

#include <stdlib.h>
int main (int argc, char **argv){

    int i, v[5], n = atoi(argv[1]);

    for(i=n; i>=0; i++)
        v[i] = i;

    return 0;
}

Podemos compilá-lo com:

$ gcc teste.c -g

Gerando o executável ./a.out. Aí invocamos o gdb com o comando:

$ gdb ./a.out

Sem os possíveis argumentos do programa. Esses argumentos serão passados quando formos rodar o programa dentro do gdb, através do comando run.

(gdb) run 4

O programa acima tem um bug que causará falha de segmentação (ou algum erro relacionado ao acesso indevido de memória). No meu caso a saída foi:

Starting program: /home/kunigami/workspace/c/a.out 

Program received signal SIGSEGV, Segmentation fault.
0x00000000004004d9 in main () at teste.c:6
6	        v[i] = i;

Ele acusa o erro na linha 6 e mostra o trecho de código correspondente. Nesse caso, podemos analisar o código com a função print (ou simplesmente ‘p’) do gdb:

(gdb) p v[i]
Cannot access memory at address 0x7ffffffff000
(gdb) p i
$2 = 832

Vemos, ao tentar imprimir v[i], que tal posição da memória é inválida e isso porque o índice i está além das 5 posições alocadas.

Breakpoints

breakpoint (b) — coloca um marcador em uma linha do código para que o programa pare. A partir daí podemos executar o código passo-a-passo para uma análise mais detalhada.

Há varias maneiras de se colocar esse marcador. Uma é você dizer diretamente o número da linha, como por exemplo:

(gdb) b 123 

Colocará um breakpoint na linha 123. Se há vários arquivos fontes que formam o executável, dá pra especificar em qual arquivo quer se colocar o breakpoint. Se por exemplo quisermos colocar no código foo.c, fazemos:

(gdb) b foo.c:123

Também dá para passar o nome da função.

(gdb) b foo.c:bar

Para listar os breakpoints setados,

info b

Associado a cada breakpoint está um identificador, que pode ser usado se quisermos remover um dado breakpoint através de

delete <id do breakpoint>

Breakpoints condicionais

Se uma dada função é chamada muitas vezes e sabemos em quais condições um erro acontecerá, poderíamos fazer a seguinte gambiarra no código:

...
if(<condição>){
/* linha L do código */
}
...

E colocar o breakpoint na linha L. Ao invés disso, podemos colocar um breakpoint condicional. Dado o id do breakpoint, a sintaxe é a seguinte:

condition <id do breakpoint> <condição>

A condição é exatamente a mesma que você colocaria dentro do if.

Andando pelo código

Uma vez que o breakpoint fez a execução do programa parar, temos vários comandos para analisar detalhadamente as próximas linhas que serão executadas.

* step (s) — vai para a próxima instrução. Se for uma função, entra nela (step into).

* next (n) — vai para a próxima instrução. Se for uma função, não entra nela (step over).

* continue (c) — continua a execução normal do código (irá para a próxima vez que encontrar um breakpoint)

* finish (sem atalho) — sai da função atual (step out)

GDB e multi-threading

Quando lidamos com código multi-thread, o GDB começa a fazer mais falta. Vamos analisar o seguinte exemplo, adaptado daqui, usando pthreads:

threads.c

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

void *foo(void *id)
{
    int i, s = 0;
    for(i = 0; i < 10000000; i++)
        s = (s + i) % 1234567;
    printf("[%d] s = %d\n", *((int*)id), s);
}

main()
{
     pthread_t t1, t2;
     int  id1 = 1, id2 = 2;

     /* Create independent threads each of which
        will execute function */
     (void) pthread_create(&t1, NULL, foo, (void*) &id1);
     (void) pthread_create(&t2, NULL, foo, (void*) &id2);

     /* Wait till threads are complete before main
        continues. Unless we wait we run the risk
        of executing an exit which will terminate   
        the process and all threads before the 
        threads have completed. */
     pthread_join(t1, NULL);
     pthread_join(t2, NULL); 

     exit(0);
}

Nesse código a thread inicial cria duas novas threads e bloqueia. Essas duas threads chamam a função foo, onde ficam em um loop. Quando ambas saem da função, a thread inicial desbloqueia e termina o código. Para compilá-lo,

gcc threads.c -lpthread -g

Vamos rodar o gdb e colocar um breakpoint na função foo, mas apenas para uma das threads. Podemos usar o breakpoint condicional.

(gdb) b foo
Breakpoint 1 at 0x400680: file teste.c, line 7.
(gdb) info b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   0x0000000000400680 in foo at teste.c:7
(gdb) condition 1 (*(int*)id) == 1

Aí rodamos o código

(gdb) run
Starting program: /home/kunigami/workspace/c/threads/a.out 
[Thread debugging using libthread_db enabled]
[New Thread 0x7ffff783c700 (LWP 30676)]
[Switching to Thread 0x7ffff783c700 (LWP 30676)]

Breakpoint 1, foo (id=0x7fffffffe2fc) at teste.c:7
7	    int i, s = 0;

Aqui ele avisa que trocou para outra thread (“Switching to Thread….“). Podemos visualizar as threads através do comando info threads.

(gdb) info threads
[New Thread 0x7ffff703b700 (LWP 30677)]
  3 Thread 0x7ffff703b700 (LWP 30677)  clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:84
* 2 Thread 0x7ffff783c700 (LWP 30676)  foo (id=0x7fffffffe2fc) at teste.c:7
  1 Thread 0x7ffff7fcf700 (LWP 30673)  clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:84

O que a tabela acima mostra é que a thread com id 1 é a thread inicial que está sendo clonada para criar a thread 3. A thread 2 está na função foo, onde ativou o breakpoint e é ela quem o GDB está analisando, indicado pelo asterisco na frente do id. Note que essa tabela provavelmente vai ser diferente a cada execução. Depende de como o sistema operacional está escalonando a execução das tarefas, por exemplo.

Vamos dar uns steps até que a variável i do loop valha, por exemplo, 3.

(gdb) p i
$2 = 3

Agora, vamos trocar para a thread 3 através do comando thread 3 e imprimir o valor de i.

(gdb) thread 3 
[Switching to thread 3 (Thread 0x7ffff703b700 (LWP 30779))]
#0  foo (id=0x7fffffffe2f8) at teste.c:8
8	    for(i = 0; i < 100000000; i++)
(gdb) p i
$7 = 2

Novamente, note que o estágio em que a thread 3 se encontra dependerá muito. Se ela ainda não tiver entrado no loop, dê alguns steps.

Conclusão

A possibilidade de se analisar várias threads ao mesmo tempo torna o GDB uma ferramenta muito mais atraente para se debugar código multi-thread do que o uso de printf‘s.

Referências

[1] http://www.delorie.com/gnu/docs/gdb/gdb_25.html


O shader OSL

Julho 10, 2011

Nessa última semana consegui avançar de maneira satisfatória no projeto. Primeiramente, implementei o shader sh_osl que é chamado pelo aplicativo rt e usa os serviços do sistema do OSL.

No post anterior mencionei que o sistema OSL exigia muitos raios por pixel e isso seria um problema ao usar o rt, mas descobri que uma opção de hypersampling que faz exatamente isso.

De maneira simplista, temos então que o aplicativo rt atira raios através da função rt_shootray e que chama uma callback chamada shadeview toda vez que um objeto é atingido. A função shadeview chama outra callback que está associada ao objeto e corresponde ao shader dele.

Por exemplo, se o objeto atingido tem o shader sh_glass, então esse objeto possui um ponteiro para a função glass_render. O que a função shadeview faz é chamar a função referenciada por esse ponteiro.

O que a função glass_render ou qualquer outra função xyz_render deve fazer é essencialmente retornar a cor do ponto de consulta (passado como parâmetro). Vou descrever então como implementei a função osl_render, correspondendo ao shader sh_osl.

Reflexão

Primeiramente, implementei a reflexão através de uma chamada recursiva da função rt_shootray, que atira um novo raio em uma dada direção. Essa parte foi fácil pois o código que eu tinha desenvolvido para o raytracer independente que eu havia implementado fazia do mesmo jeito.

A figura a seguir é a renderização da cena da caixa de cornell usando um shader osl de espelho para a caixa maior e um shader brl-cad de xadrez para a caixa menor.


Figura 1: Teste de integração de um shader OSL e um shader BRL-CAD

No final das contas não tive que tomar nenhum cuidado adicional ao misturar shaders OSL e BRL-CAD, ao contrário do que eu havia dito no post anterior.

Refração

Implementar a refração foi mais complicado. Primeiramente, descobri que o detector de colisões do BRL-CAD sempre calcula dois pontos de interseção para cada superfície. Um, chamado inhit, é o ponto da superfície no qual o raio bate inicialmente. O outro, chamado outhit, supõe que o raio foi atirado de dentro da superfície.

A função osl_render já é chamada com o ponto de interseção P equivalente a inhit, pois quem seta isso é o shadeview. Entretanto quando o raio é de refração (interno) eu quero que P seja o outhit.

Portanto, tive que escrever uma nova callback para tratar raios refratados. Sempre que um raio for retornado pelo OSL, verifico se ele é refratado. Para isso, basta fazer um produto escalar entre a normal e o novo raio para ver se eles apontam em direções “opostas”. Em caso positivo, mudo a callback que será chamada a próxima vez que um objeto for atingido no rt_shootray. O que essa callback faz é setar P como outhit e chamar xyz_render.


Figura 2: Teste com o shader vidro

Implementando novos shaders para o OSL

Para verificar se a interface sh_osl está computando os dados corretamente, decidi criar novos shaders que utilizam esses dados. Um exemplo é o shader xadrez, que usa as coordenadas do mapeamento uv.

Aproveitei para implementar suporte aos grupos de shaders, comentado em um post anterior, além da possibilidade de setar os parâmetros via a interface mged.

Agora é possível definir um grupo de shaders através da própria interface. O grupo de shaders é então construído a partir de shaders primitivos. Por exemplo a descrição do shader

shadername gen_color#layername#c1#base#color#1.0#0.0#1.0
shadername gen_color#layername#c2#base#color#0.0#1.0#0.0
shadername checker#layername#chk#K#float#4.0
join c1#Cout#chk#Cw
join c2#Cout#chk#Cb

foi usada no caixa menor, como ilustra a imagem abaixo:


Figura 3: Teste com shader xadrez verde-magenta

A sintaxe ficou meio feia mas está funcional. A primeira linha descreve um shader de cor genérica inicializado com a cor magenta enquanto o da segunda linha possui a cor verde. O terceiro shader descrito é o shader xadrez propriamente dito, que usa a saída de dois outros shader’s para compor a cor dos quadrados.

A quarta e quinta linhas tratam de ligar a saída do primeiro e segundo shader à entrada do shader xadrez.

A vantagem de descrever shaders dessa maneira, é que se por exemplo eu quiser compor o shader xadrez com uma cor verde e com um shader de espelho é só mudar para

shadername mirror#layername#c1
shadername gen_color#layername#c2#base#color#0.0#1.0#0.0
shadername checker#layername#chk#K#float#4.0
join c1#Cout#chk#Cw
join c2#Cout#chk#Cb

Que o resultado será:


Figura 4: Teste com shader xadrez verde-espelho

Além do mais, a parede do fundo da Figura 3 usa o shader “nuvem” adaptado de um shader BRL-CAD. Não parece nem um pouco com nuvem, mas é uma textura procedural, então não dá pra fazer milagre.

Próximos passos

Preciso resolver a questão do multi-threading. Semana que vem pretendo estudar como funcionam as funções de aquisição e liberação de recursos providas pelo BRL-CAD e qual recurso do sistema OSL eu preciso garantir acesso exclusivo.

Mais pra frente gostaria de implementar um modo de visualização incremental. Atualmente se quisermos visualizar a imagem enquanto ela é renderizada, temos que esperar todas as amostragens serem feitas para cada pixel antes de vê-lo. Minha ideia é que a cada amostragem, toda a imagem seja atualizada, de forma que inicialmente vejamos um monte de pixels dispersos e conforme mais e mais amostragens sejam feitas, a imagem vá convergindo para uma sem ruídos.

Idealmente, gostaria também de implementar uma interface para fazer a composição do grupo de shaders OSL. A GUI do BRL-CAD é escrita em Tcl, linguagem que eu teria que estudar antes de mais nada. Creio que não conseguirei fazer isso antes do término do programa, mas pretendo fazer isso algum dia.


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.