Bill Allombert on Tue, 23 Apr 2024 17:00:57 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Computing Sathe-Selberg


On Tue, Apr 23, 2024 at 08:04:33AM -0400, Charles Greathouse wrote:
> The Sathe-Selberg theorem uses a function G to improve the range of the
> classical Landau estimates for the density of k-almost primes. It is an
> infinite product over primes; its truncation is
> 
> G(z,lim=1e4)=prodeuler(p=2,lim,(1+z/p)*(1.-1/p)^z)/gamma(z+1)
> 
> What is a better way to compute this?

You can try this:

G(z)=
{
  my(a=max(2,ceil(abs(z))));
  my(lim=4*a);
  my(e=getlocalbitprec());
  localbitprec(2*e);
  my(S=sum(n=2,e,-(z+(-z)^n)/p^n/n));
  prodeuler(p=2,lim,(1+z/p)*(1.-1/p)^z)*exp(sumeulerrat(S,,lim+1))/gamma(z+1);
}

If z is a small integer, you can use prodeulerrat:
Gi(z)=prodeulerrat((1+z/p)*(1-1/p)^z)/gamma(z+1)

if z=a/b is a rational of small height, you can also do:
Gr(z)=
{
  my(a=numerator(z),b=denominator(z));
  prodeulerrat((1+z/p)^b*(1-1/p)^a)^(1/b)/gamma(z+1)
}

Cheers,
Bill.