Loading ...
Sorry, an error occurred while loading the content.

16052Re: [PrimeNumbers] Sieving primo-factorials

Expand Messages
  • Décio Luiz Gazzoni Filho
    Feb 8 9:24 PM
    • 0 Attachment
      On Wednesday 09 February 2005 01:59, you wrote:
      > Can you combine this fast caching method with using the prime factorization
      > of n as you suggest (except for initialization if you don't start at n=1)?
      > I considered the factorization method but thought my suggestion would be
      > much faster.

      I coded the sieve (code reproduced below). Despite the fact that it segfaults
      with gp2c, so that I had to run it interpreted instead of compiled, it's
      still a speed demon. With 10 minutes of CPU time (P4 2.6 GHz here) I
      eliminated nearly 4000 composites. I shudder to think that this could
      probably have cost me somewhere in the vicinity of 10 CPU-days.

      I'm sure I commited a few egregious errors down there, not to mention some
      ugly hacks, so please refrain from flaming me (:

      I wonder, at what point should I switch back to trial factoring?

      Décio

      ---
      primorial(x) =
      {
      local(a);
      a = 1;
      forprime(p=2,x,a*=p);
      return(a);
      }

      \\ Reverse indexes the primes -- the inverse of PARI's prime() function
      init_invprime(nmax) =
      {
      local(last,k);
      global(invprime);
      invprime = vector(nmax);
      last = 1;
      k = 0;
      forprime(p=2,nmax,
      k++;
      for(i=last,p-1,
      invprime[i] = k-1;
      );
      last = p;
      );
      }

      \\ Computes the prime factorization of n!/n#
      legendre_factorial(n) =
      {
      local(v,i,k);
      \\ v[] is initialized with -1 to simulate division by n#
      v = vector(invprime[n],i,-1);
      i = 0;
      forprime(p=2,n,
      i++;
      k = 1;
      while((cont = floor(n/p^k)) > 0,
      v[i] += cont;
      k++;
      );
      );
      return(v);
      }

      sieve(nmin,nmax,pmax) =
      {
      local(v,i,j,a,b,apb,amb);
      \\ Integers of the form n! +/- n# +/- 1 are not divisible by any
      \\ prime up to n, and since the minimum value of n that we're
      \\ considering is nmin, then we'll start there. Sure, we'll waste
      \\ a little time checking whether the primes from nmin to a given n
      \\ divide n! +/- n# +/- 1, for n > nmin, but that's not too costly
      \\ anyway.
      j = 0;
      forprime(p=nmin,pmax,
      if(bitand(j++,255)==0,
      print("Sieving "p"...");
      );
      v = legendre_factorial(nmin);

      \\ a == nmin!/nmin# (mod p)
      a = Mod(1,p);
      i = 0;
      forprime(q=2,floor(nmin/2),
      a *= Mod(q,p)^v[i++];
      );

      \\ b == nmin# (mod p)
      b = Mod(1,p);
      forprime(q=2,nmin,
      b *= Mod(q,p);
      );

      \\ Now a == nmin! (mod p)
      a *= b;

      for(i=nmin,nmax,
      apb = a + b;
      amb = a - b;

      if(apb == Mod(p-1,p),
      write("composite_"nmin"_"nmax, i"! + "i"# + 1, "p);
      );
      if(apb == Mod(1,p),
      write("composite_"nmin"_"nmax, i"! + "i"# - 1, "p);
      );
      if(amb == Mod(p-1,p),
      write("composite_"nmin"_"nmax, i"! - "i"# + 1, "p);
      );
      if(amb == Mod(1,p),
      write("composite_"nmin"_"nmax, i"! - "i"# - 1, "p);
      );

      a *= Mod(i+1,p);
      if(isprime(i+1),
      b *= Mod(i+1,p);
      );
      );
      );
      }


      [Non-text portions of this message have been removed]
    • Show all 9 messages in this topic