| 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.