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