Re: Q: Dickman function

*To*: mathgroup at smc.vnet.net*Subject*: [mg21036] Re: Q: Dickman function*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Sun, 12 Dec 1999 23:51:06 -0500 (EST)*Organization*: Wolfram Research, Inc.*References*: <819je5$81p@smc.vnet.net> <8263fn$d77@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Francois Grieu wrote: > I'm trying to build the Dickman function, defined by > > /1 > F(x)=1 for x>=1 F(x)=1-| F(t/(1-t))/t dt for 0<=x<=1 > /x > > (ref: Knuth's TAOCP, vol 2, 4.5.4, p383 in third edition) > > After a lot of trial and error, I came up with a tentative > minimal definition in Mathematica: > > F[x_]:=F[x,Ceiling[1/x]-1] > F[x_,0]:=1 > F[x_,n_]:=F[1/n,n-1]-Integrate[F[t[n]/(1-t[n]),n-1]/t[n],{t[n],x,1/n}] > > This does works for x>=1/3 both numericaly and symbolicaly. > However both F[0.3] and F[1/4] give an imaginary component, which > is an 'obvious' nonsense; it's numerical value is vanishingly > small, but the expression of F[x,3] is so complex that even > FullSimplify won't cut it. And it's sooo long when x gets small. > > More trial and error leads to > > F[x_]:=F[x,Ceiling[1/x]-1] > F[x_,0]:=1 > F[x_,1]:=1+Log[x] > F[x_,2]:=1-Pi^2/12+Log[x]+Log[x]^2/2+PolyLog[2,x] > F[x_,n_]:=F[1/n,n-1]-NIntegrate[F[t[n]/(1-t[n]),n-1]/t[n],{t[n],x,1/n}] > > Extending the usable range to maybe x>=0.15, but I am afraid I won't > ever get F[0.1] that way. > > Any clues ? > > Francois Grieu I can show ways to compute it quickly in the range of interest. Computing it accurately is another matter entirely, but I'll give ideas for that as well. One reason it takes so long to compute is the endless recursion entailed by the definition. Often one can quell this by memo-izing results e.g. F[x_,n_] := F[x,n] = ... but I doubt that will help here because NIntegrate may choose different evaluation points at various different stages (I confess I did not try this, though). Anyway one can in essence eliminate the massive recursion by storing results in interpolating functions. A natural way to do this is to use the numeric differential equation solver NDSolve, because it returns results in exactly this form. I preferred to work with the interval [1,infinity) rather than (0,1] so I did the transformation g[x] = f[1/x]. The relevant ode is g'[x] = -(1/x^2)*f'[1/x] = -(1/x^2)*f[(1/x)/(1-1/x)]/(1/x) = -f[1/(x-1)]/x = -g[x-1]/x with initial condition g[1] = 1. Note that there may be numerical benefits to working with the original interval. I've no idea if that is so, though I am guessing it is not the case. But I thought I should mention these various possibilities before proceeding. The code below will create interpolation functions for each unit interval and then use them to evaluate at values of interest. This includes both values explicitly requested by the user and values needed to get from one interval to the next by NDSolve. Clear[g,ifun]; g[x_?NumberQ] /; x<=1 := 1 g[x_?NumberQ] /; 1<x<=2 := 1 - Log[x] g[x_?NumberQ] /; x>2 := ifun[Ceiling[x]][x] Clear[ifun] ifun[n_Integer] /; n>2 := ifun[n] = First[r /. NDSolve[{r[n-1]==g[n-1], r'[x]==-g[x-1]/x}, r, {x,n-1,n}]] The examples indicate this is fast but not accurate beyond 5 or so (by comparison to F[1/5]; they agree to two significant digits). One can improve this by playing with the precision controls. For example, try ifun[n_Integer] /; n>2 := ifun[n] = First[r /. NDSolve[{r[n-1]==g[n-1], r'[x]==-g[x-1]/x}, r, {x,n-1,n}, AccuracyGoal->20, PrecisionGoal->20, WorkingPrecision->30]] We get good agreement between g[7] and F[1/7] with this. They are also comparable in speed. The advantage of the NDSolve approach is that each subsequent unit interval will take roughly constant time (about 3 seconds on my machine), whereas the speed of F as given above deteriorates much worse than linearly. Hence it is hard to tell whether the code above continues to agree with F much beyond that point. Here are some results for later reference. In[40]:= Timing[g6a = g[6]] Out[40]= {14.94 Second, 0.000019649696353955308876418005453} In[41]:= Timing[g10a = g[10]] -11 Out[41]= {11. Second, 2.770171839308505982248911625 10 } A reasonable thing to try now is to raise the precision controls and see if new results agree with old, say at x=10. To this end, I did ifun[n_Integer] /; n>2 := ifun[n] = First[r /. NDSolve[{r[n-1]==g[n-1], r'[x]==-g[x-1]/x}, r, {x,n-1,n}, AccuracyGoal->30, PrecisionGoal->30, WorkingPrecision->60]] Unfortunately this hits some deficiencies in the internal NDSolve code. Rob Knapp, who knows alot more about numerically solving odes than I do, explained that in any case it is difficult to attain that level of global control of accuracy and precision, and numeric integration might be much better suited to handling this issue. So the next step was to recast this as an integration problem (trivial: it was such to begin with) and to use interpolating functions in that context (not as easy but I'll show something that works). Again we'll start with basic default precision controls. Also I worked again with g[x] = F[1/x], transforming the integral as appropriate. The idea is to find the greatest integer n less than x (the Floor[]), start at g[n], and integrate up to x. For each unit interval we create an interpolation function from 101 evenly spaced evaluation points. I used a higher than customary interpolation order to keep things smooth for NIntegrate use. A small amount of fancy work was needed to avoid infinite recursion at the endpoints, and quite likely this code can be simplified. Clear[g2,ifun2] g2[x_] /; x<=1 := 1 g2[x_] /; 1<x<=2 := 1 - Log[x] g2[x_] /; 2<x<=3 := 1 - Log[2] - NIntegrate[g2[t-1]/t,{t,2,x},Compiled->False] ifun2[n_Integer] /; n>=2 := ifun2[n] = Interpolation[ Table[{n+j/100,g2[n+j/100]}, {j,0,100}], InterpolationOrder->6] g2[x_] /; 3<x && Floor[x]==x := g2[x] = With[{n=Floor[x]}, ifun2[n-2][n-1] - NIntegrate[ifun2[n-2][t-1]/t,{t,n-1,n},Compiled->False]] g2[x_] /; 3<x := With[{n=Floor[x]}, ifun2[n-1][n] - NIntegrate[ifun2[n-1][t-1]/t,{t,n,x},Compiled->False]] In[8]:= Timing[g6b = g2[6]] Out[8]= {1.81 Second, 0.0000196497} In[9]:= Timing[g10b = g2[10]] -11 Out[9]= {2.46 Second, 2.77016 10 } (Why did I get a spelling warning for g10b but not g6b? Is it based on percentage of common letters? Who knows?). One good thing is that g10b agrees to several places with g10a. So I think we can say with some assurance that are in the right ballpark (also one can check that the rate of decay is consistent with that from earlier unit intervals, and on those we have agreement with the original definition for F[]). Various changes to NIntegrate option settings did little to change the result, although they played havoc with the speed. Lowering the interpolation order in the constructed interpolating functions gave a different result at x=10 although still of the same order of magnitude. So I am inclined to believe that what we have is correct to a few places, or at worst incorrect but of the right magnitude. On a different tack, I suppose one could run a Monte-Carlo simulation by reversing the original intent of the function. Specifically, we use the fact that for large integers x, F[a] approximates the probability that all divisors of x are less than x^a. So we pick, say, 10^12 random 40 digit integers. Factor them. Then see what fraction have all factors less than the tenth root of the number selected (this process might take a while....). More seriously, I did a scaled-down version of this to approximate F[1/3]: ee = Table[Random[Integer,10^14], {1000}]; Timing[ff = Map[{#,FactorInteger[#,FactorComplete->False]}&, ee];] gg = Map[MapAt[Map[First,#]&,#,2]&, ff]; (* kill off multiplicities *) hh = Select[gg, Max[#[[2]]]^3<#[[1]]&]; Length[hh] This takes about 5 seconds. In[235]:= Length[kk] Out[235]= 47 This means the fraction for this trial was 47/1000, which agrees quite well with g[3] = 0.048610. Subsequent trials netted results of 41, 43, and 51 which are not as good but still in the right ballpark. To estimate F[1/4] in a similar manner I did ee = Table[Random[Integer,10^16], {10000}]; ... hh = Select[gg, Max[#[[2]]]^4<#[[1]]&]; Length[hh] This took around a minute on my machine and gave a result of 40. So we have 40/10000 as an estimate for F[1/4] (actual value: 0.00491407), which is not unreasonable for this sort of simulation. Subsequent runs gave 41 and 52. Daniel Lichtblau Wolfram Research

**Re: Re: Abs[a] Sin[Abs[a]]**

**Minimum of f(x,y) ?**

**Re: Q: Dickman function**

**Re: Adding interpolated values to a list**