mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Math (https://www.mersenneforum.org/forumdisplay.php?f=8)
-   -   Help with a Dickman's Rho integral (https://www.mersenneforum.org/showthread.php?t=25830)

preda 2020-08-12 04:56

Help with a Dickman's Rho integral
 
Hi, could somebody help me approximate this integral:
(where rho() is Dickman's Rho function)

integral(rho((log(x)-log(p))/log(B1)) / (x*log(x)), x = 2^N to inf)

This is what I tried myself (may be errors):

1. variable change t=log(x):
integral(rho((t - a)/b) / t, t = N to inf)
(where a==log(p), b==log(B1)).

2. tried a few things using rho()'s property:
rho'(x)=-rho(x-1)/x
but no success to integrate it based on that.

3. variable change x=(t-a)/b
integral(rho(x) / (b*x+a), x = (N-a)/b to inf)

4. approximate rho(x)=x^(-x)
integral(1/(x^x * (b*x + a)), x)

and here's where I'm at. It seems neither Sage (SageMath) nor SymPy can symbolically integrate the above. (Sage even has dickman_rho() builtin).

I'm looking for an approximation of the above integral. The approximation does not need to be precise, let's say withing a few percent (error < 5%). And my integration skill are.. let's say very elementary. A pity I couldn't get any symbolic help from Sage or SymPy (but, again, I'm very much a beginner there too).

Thanks!

R. Gerbicz 2020-08-12 12:06

[QUOTE=preda;553385]Hi, could somebody help me approximate this integral:
(where rho() is Dickman's Rho function)
[/QUOTE]

You can get the rho function's Taylor polynom for each n in the [n-1,n] interval, if you want the first n=maxn polynom, then just call ff(maxn,maxe) if you want the max degree=e in the polynoms, I've centered each polynom at n-0.5.
[CODE]
ff(maxn=50,maxe=100)={A=matrix(maxn,maxe+1,i,j,0.0);
A[1,1]=1.0;
for(n=2,maxn,
for(i=0,maxe,for(l=0,maxe-i-1,
e2=i+l+1;coeff=-A[n-1,i+1]*(-1)^l*(n-0.5)^(-l-1)/e2;
A[n,e2+1]+=coeff));
A[n,1]+=sum(e=0,maxe,A[n-1,e+1]*0.5^e-A[n,e+1]*(-0.5)^e))}

ff();
F(x,maxe=100)={n=ceil(x);return(sum(e=0,maxe,(x-n+0.5)^e*A[n,e+1]))}
? F(1.3)
%4 = 0.73763573553250894796450401311904560284
? 1-log(1.3)
%5 = 0.73763573553250894796450401311904560280
? F(2)
%6 = 0.30685281944005469058276787854182343204
? F(3)
%7 = 0.048608388291131566907183039343407421419
? F(10)
%8 = 2.7701718377259589887581212167960043114 E-11
? F(20)
%9 = 2.4617828364115816365 E-29
? F(40)
%10 = 6.381005074572958815 E-38
? F(50)
%11 = 6.381005074572958815 E-38
?
[/CODE]
(compare these values to [url]https://en.wikipedia.org/wiki/Dickman_function[/url])
Don't know how far you want these values, F(50) is already a very small number (reached machine precision even for n=25). In A[n,] you can see the coefficients for the n-th polynom, with that you can get the rho function in [n-1,n].
rho(x)=F(x)=sum(e=0,maxe,(x-n+0.5)^e*A[n,e+1]) for n=ceil(x).
Then use any method, if x falls within the same [n-1,n] interval then you need to integrate only a polynom (ofcourse we are only approximating but if maxe isn't that small we can get very good approx.)
t shouldn't run from N*log(2) in method=2?

preda 2020-08-13 05:00

[QUOTE=R. Gerbicz;553414]t shouldn't run from N*log(2) in method=2?[/QUOTE]

Yes. Fixing that, and changing the log() from base 'e' to base 2, we get:
integral(1/x * rho((x - a)/b), x, N, infinity),
where a=log2(P), b=log2(B1). i.e. the integral has the nice property of being invariant relative to the base of the log().

Thank you for the interesting Pari code. Nevertheless, it still only offers a way for a numeric evaluation of the integral, not an approximate analytic formula. Well, that's fine, probably I'll proceed numerically then.

R. Gerbicz 2020-08-13 10:24

[QUOTE=preda;553499]
Nevertheless, it still only offers a way for a numeric evaluation of the integral, not an approximate analytic formula. Well, that's fine, probably I'll proceed numerically then.[/QUOTE]

Yes, a numeric evaluation, integration is quite easy.
But an analytic formula for a fixed [n-1,n] (or subinterval) is still possible, I've done exactly the same way to get the above rho(x) when done integrate(-rho(x-1)/x), for this you need to get the Taylor polynom for 1/x with center=a, and it is possible:
1/x=sum(e=0,inf,(-1)^e*a^(-e-1)*(x-a)^e) for abs(x-a)<abs(a)

with Pari-Gp and with maxe=100:
[CODE]
? G(x,a)=sum(e=0,100,(-1)^e*a^(-e-1)*(x-a)^e)
%1 = (x,a)->sum(e=0,100,(-1)^e*a^(-e-1)*(x-a)^e)
?
? G(7.129,7.593)
%2 = 0.14027212792818067050077149670360499368
? 1/7.129
%3 = 0.14027212792818067050077149670360499369
?
[/CODE]


All times are UTC. The time now is 18:22.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.