![]() |
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! |
[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? |
[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. |
[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.