Fatorial revisitado II

Uma propriedade interessante dos números surgiu nas minhas investigações do fatorial. O número de vezes que um número divide seus múltiplos entre duas potências suas consecutivas é igual à primeira potência menos um mais o expoente dessa primeira potência. Bastante intuitivo, eu penso.

Entre 2ⁿ e 2ⁿ⁺¹, por exemplo, temos os múltiplos 2ⁿ+2, 2ⁿ+4, 2ⁿ+6, e assim por diante. Cada um desses números pode ser dividido por uma potência de 2. E a soma dos expoentes é 2ⁿ-(1+n).

Tomemos outro exemplo com 5. Entre 5¹ e 5² temos os seguintes múltiplos de 5: 10 (2*5); 15 (3*5); e 20 (2*2*5). A soma dos expoentes de 5 é 3, que é igual a 5-2, ou, 5ⁿ-(1+n)=5¹-(1+1)=3.

Essa propriedade podemos usar para acelerar o cáculo do fatorial. Infelizmente, ainda precisamos manter o método anterior para calcular os expoentes de cada primo em cada múltiplo posterior à última potência anterior a n (quando estivermos calculando n!).

Então, o algoritmo ficou divido em 3 partes:

1. Para cada primo, marcar seus múltiplos menores que n+1;
2. Para cada potência de cada primo p, calcular p^m-(1+m);
3. Para os múltiplos de p superiores à última potência inferior a n, calcular os expoentes de p.

Em Java, montei o seguinte:

  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) {
        //Marcar todos os múltiplos
        for(int j=i+i; j<m; j+=i) {
          p[j]=-1;
        }
        //Calcular os expoentes entre potências consecutivas
        int k=i;
        int power=1;
        for(; k*i<=m; k*=i) {
          p[i]+=k-power-1;
          p[i]+=power;
          power++;
        }
        //Calcular os expoentes para os múltiplos que seguem a última potência
        if(k<m) {
          p[i]+=power;
          for(int l=k+i; l<m; l+=i) {
            p[i]+=in(i, l);
          }
        }
      }
    }
    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;
  }

O mais surpreendente é que funciona.