Fatorial revisitado

Eu estava tendo dificuldades para criar um executável com um banco de dados embutido (com Tcl, SQLite e Freewrap), quando resolvi brincar um pouco com o Tcl. Escrevi uma pequena função para calcular fatoriais. Tcl permite gerar números bem grandes (como 1000!) com facilidade. Para tornar a coisa mais interessante, resolvi escrever uma versão sem if (abaixo).

proc fac {n} {
  set e 1;
  for {set x 2} {$x<=$n} {incr x} {
    append e "*$x"
  }
  puts $e
  return [expr $e]
}

A função simplesmente gera uma expressão (por exemplo, ela gera "2*3*4*5" para 5!) e depois calcula o resultado. Talvez inspirado pela recente releitura o livro Metamat! (de Gregory Chaitin), achei que seria interessante encontrar uma maneira de encolher essa expressão e assim tornar o programa mais curto e eficiente.

A solução mais próxima seria gerar um produto de primos, assim a multiplicação teria o menor número possível de multiplicandos. Não demorou muito para sair um algoritmo parecido com o Crivo de Eratóstenes. Se eu quiser calcular n!, faço o seguinte:

1. Enumero todos os inteiros até n;
2. Começando em 2, marco todos os múltiplos dos primos e anoto qual potência de cada primo divide cada múltiplo (e vou somando as potências);
3. Ao fim, terei uma potência para cada primo e a multiplicação de cada um elevado à sua potência produzirá o fatorial.

Por exemplo, para 6!, os passos são:

1. Começando em 2, descubro que 2=2^1, 4=2^2 e 6=2^1*3;
2. Com 3, tenho 3=3^1 e 6=3^1*2;
3. Quatro posso ignorar, porque já o visitei e anotei como inteiro composto;
4. 5 não tem múltiplos menores ou iguais a 6, então anoto 5^1;
5. Finalmente, 6 já sei que não é primo.

O resultado é 2^4*3^2*5 e isso produz 720, conforme esperado. Abaixo apresento uma solução em Java:

  private static int in(int p, int i) {
    int n=0;
    while(i%p==0) {
      n++;
      i=i/p;
    }
    return n;
  }

  public static BigInteger f(int n) {
    int m=n+1;
    int[] p=new int[m];
    for(int i=2; i<=m/2; i++) {
      if(p[i]!=-1) {
        p[i]=1;
        for(int j=i+i; j<m; j+=i) {
          if(j%i==0) {
            int v=in(i,j);
            p[i]+=v;
            p[j]=-1;
          }
        }
      }
    }
    BigInteger f=BigInteger.ONE;
    for(int k=2; k<m; k++) {
      int pk=p[k];
      if(pk>=0) {
        if(pk==0) {
          f=f.multiply(BigInteger.valueOf(k));
        } else {
          f=f.multiply(BigInteger.valueOf(k).pow(pk));
        }
      }
    }
    return f;
  }

Esse código tem uma pequena otimização. Ele percorre os primos apenas até n/2. Todos os primos acima disso só vão aparecer uma única vez no produto final (em 10!, por exemplo, 7 só aparece como 7^1). O array p começa com zeros. As posições ocupadas por primos vão recebendo as respectivas potências e as demais recebem -1. Ao fim, as posições de primos acima de n/2 vão continuar com 0 (porque não as visitamos), então a segunda parte da função (que faz as multiplicações para calcular o resultado final) ignora os valores menores que zero, troca os zeros por uns e usa os demais valores como os encontrar.

Comparando com uma função simples que multiplica instâncias de BigInteger de 2 até n, para valores pequenos, a diferença de tempo é insignificante. Em algum ponto entre 1000! e 2000!, a nova função passa a tomar um pouco mais que a metade do tempo. Ela é mais econômica em memória também.

O método f() retorna um BigInteger, mas imagino que em algumas situações possa ser mais útil guardar o produto de primos ou simplesmente imprimir a expressão.