Re: a tricky limit

*To*: mathgroup at smc.vnet.net*Subject*: [mg15650] Re: [mg15327] a tricky limit*From*: Jurgen Tischer <jtischer at col2.telecom.com.co>*Date*: Sat, 30 Jan 1999 04:28:45 -0500 (EST)*Organization*: Universidad del Valle*References*: <199901080915.EAA03988@smc.vnet.net.>*Sender*: owner-wri-mathgroup at wolfram.com

Once more the tricky limit. I worked out a different solution to calculate more digits. It's a bit involved so here it goes: The problem was to solve Product[(1-x^p[n]/(n+1)),{n,2,inf}]/(1-x). where p[n]=n-n/Divisors[n][[2]]. Taking Log it changes to the sum -Log[1-x]+Sum[Log[1-x^p[n]/(n+1)],{n,2,inf}]. Now use the series for Log[1-x] and Log[1-x^p[n]/(n+1)] to get Sum[x^n/n,{n,inf}]-Sum[(x^p[k]/(k+1))^n/n,{n,inf},{k,2,inf}]. Now for x<1 everything is absolutely convergent and wi can reorder to get: x + x^2/2 -Sum[(x^p[k]/(k+1))^n / n,{k,2,inf},{n,2,inf} + Sum[x^(n+1)/(n+1)-x^p[n]/(n+1),{n,2,inf}]. Let's start with one of the first sums, say for n==2. A lower bound at x is given by In[1]:= Sum[(x^k/(1 + k))^2/2, {k, 2, Infinity}] Out[1]= (-4*x^2 - x^4 + 4*PolyLog[2, x^2])/(8*x^2) In[2]:= Limit[%, x -> 1] Out[2]= 1/24*(-15 + 2*Pi^2) and an upper bound by In[3]:= Sum[(x^(k/2)/(k + 1))^2/2, {k, 2, Infinity}] Out[3]= (-4*x - x^2 + 4*PolyLog[2, x])/(8*x) In[4]:= Limit[%, x -> 1] Out[4]= 1/24*(-15 + 2*Pi^2) So the limits are the same, and the same as Sum[1/k^2/2,{k,2,inf}], so by monotone convergence we may interchange limits. Inspection of higher n's (or knowledge of Riemann's Zeta function, see Help on Zeta), leads to the formula Limit[Sum[(x^p[k]/(k+1))^n/n ,{n,2,inf}],x->1] == Zeta[n]-1-2^-n and we need to sum Sum[Zeta[n]-1 -2^-n,{n,2,inf}] Mathematica doesn't want to solve it but in parts: In[5]:= Sum[-(1/(2^n*n)), {n, 2, Infinity}] Out[5]= 1/2*(1 - 2*Log[2]) In[6]:= Sum[(Zeta[n] - 1)/n, {n, 2, Infinity}] Out[6]= 1 - EulerGamma In[7]:= s1 = 1 + 1/2 - 1/2*(1 - 2*Log[2]) + 1-EulerGamma)//Simplify Out[7]= EulerGamma+Log[2] We are left with Limit[Sum[x^(x+1)/(n+1),{n,2,inf}]-Sum[x^p[n]/(n+1),{n,48}],x->1]. The idea now is to extract from the sum the terms of one same second divisor of n, similar to the sieve for finding primes. One problem is that there is no closed formula for the sum over all terms with the same divisor, one has to split the sum into sums of the type Sum[x^(n(k-1)/k),{k,a,inf,b}], and even for small p this splits into many sums. Now: In[8]:= Sum[x^(n + 1)/(n + 1) - x^((n*(k - 1))/k)/(n + 1), {n, a, Infinity, b}] Out[8]= (x^(1 + a)*LerchPhi[x^b, 1, (1 + a)/b])/b - (x^((a*(-1 + k))/k)* LerchPhi[(x^((-1 + k)/k))^b, 1, (1 + a)/b])/ b Inspection of cases for different values of a and b shows that the limit of the sum does not depend of a and is given by 1/b Log[1 - 1/k]. So the sum of all terms with second divisor k=Prime[m] is equal to c/bLog[1-1/k], where b is the period for k, which is equal (without proof, found by inspection) to the product \[CapitalPi] of all primes up to k (included), and c is the number of multiples of k in [k,k+P-1] with second divisor k. I calculated the first 8 c's and then searched the sequence at Sloane's, (again without proof) c is equal to the product of the primes-1 up to k-1. So we have for upper and lower bounds ub, lb: ub[n]: Sum[c[k]/b[k] Log[1-1/k],{k,n}] lb[n]: ub[n]-Log[1-1/(n+1)]+Sum[c[k]/b[k]](Log[1-1/(n+1)])== ub[n]-(1-Sum[c[k]/b[k]])Log[1-1/(n+1)] where the meaning of the symbols is {n,p,s,s2}== {n, c[n]/b[n], sum[c[k]/b[k],{k,n-1}], sum[c[k]/b[k] Log [1-1/k],{k,n-1}]}. After summing up 106100002 terms (it took 45 hours) I found the following bounds for the limit: ub=2.29217369524, lb=2.29217369419, this is an error of 1.06 10^-9. but I still don't know the 8th digit. At least my "best bet" 2.29217 coincides with this solution. Of course I was willing to spend another (estimated) 45 hours to find this last digit, but there came a new problem. Somewhere between 105 million and 106 million the speed went down so calculation time of 100 000 terms rose from 100 seconds to 5000 seconds. This would mean some 100 days of calculating. (It seems that the algorithm for Prime changes about that number.) So that's the end for me for the moment. Jurgen If you happen to have a faster computer, you could run it on the following: {n,p,s,s2}={106100002,1.20410332625684268556672671364644320606271`30.243*^-11, 0.97388264577427072385607700884973105920738`30.3338,-0.\ 440862266213089289557616316841531205845154`30.6699} li={}; While[True,{n,p,s,s2}=Nest[fu,{n,p,s,s2},100000]; AppendTo[li,{n,p,s,s2}];Print[{Exp[s1+s2],Exp[s1+s2+s Log[1-1/(n+1)]]}]] The Print is just so you know your calculation is still alive and if you can stop. Arnold Knopfmacher wrote: > > I wish to obtain a numerical estimate (say 8 decimal digits) of the > limit as x tends to 1 from below of the function > h[x]=(Product[(1-fm[x]/(m+1)),{m,2,Infinity}])/(1-x) where > fm[x]=x^(m-m/d) and d is the smallest divisor of m that is greater than > 1. The problem is that when I replace Infinity by say 1000 as the upper > limit of the product, the function blows up near 1. Visual inspection > of the graph of h[x] for 0<x<0.9 say, suggests that the limit should > have a value around 2.1. Can anyone help? > > Thanks > Arnold Knopfmacher > Dept of Computational and Applied Math Witwatersrand University > South Africa

**References**:**a tricky limit***From:*"Arnold Knopfmacher" <arnoldk@gauss.cam.wits.ac.za>

**Re: Block vs. Module**

**Re: Button to delete all graphics**

**Re: a tricky limit**

**Re: a tricky limit**