computing Dickman's function
- To: mathgroup at smc.vnet.net
- Subject: [mg21322] computing Dickman's function
- From: fgrieu at micronet.fr (Francois Grieu)
- Date: Fri, 24 Dec 1999 03:42:26 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
In a previous message [mg20909] I asked for ways to numericaly compute
Dickman's function, which can be defined as (view in non-prop. font)
/1
F[x] = 1-| F[t/(1-t)]/t dt when 0<=x<=1 F[x] = 1 when x>=1
/x
This gives the proportion F[x] of the integers N such that the highest
prime factor of N is less that N^x. [Reference: Donald E. Knuth,
The Art of Computer Programming, volume 2, section 4.5.4, page 383
in third edition].
Thanks to all who answered, especially to Daniel Lichtblau, which ideas
inspired the following work.
A notable simplification arise by defining f[y] = F[1/y] :
/y
f[y] = 1-| f[t-1]/t dt when y>=1 f[y] = 1 when 0<=y<=1
/1
We can compute f[y] on sucessive segments of length 1. For a start,
f[y] = 1-Log[y] when 1<=y<=2. It is possible to derive symbolic
formulas on one or two more intervals, but they grow much too long.
It is possible to recursively integrate numerically, but execution time
is beyond reason for y>8. It is better to devise some approximation on
an interval, and integrate to obtain an approximation on the next
interval. A method that works well is a power series expansion to some
degree, at points t = n-1/2. In functional programming style:
Clear[DickmanF, Dickman$f, Dickm$pr, Dickm$ex, Dickm$dv, Dickm$ma];
DickmanF::usage =
"DickmanF[x] is the proportion of integers N with factors below N^x";
Dickman$f::usage =
"Dickman$f[y] is the proportion of integers N with factors below N^(1/y)";
Dickm$pr = 40 (* desired ABSOLUTE (rather than relative) accuracy *);
Dickm$ex = Ceiling[Dickm$pr*Log[3, 10]]+1 (* order of expansions *);
Dickm$ma = Ceiling[Dickm$ex/Log[Dickm$ex]] (* domain limit *);
(* Dickm$dv[n] approximates Dickman$f[t + n - 1/2] for -1/2<=t<=1/2 *)
Dickm$dv[1] = N[1, Dickm$pr+1];
Dickm$dv[n_] := Dickm$dv[n] = Module[{u = Dickm$dv[n - 1], v},
v = Integrate[(u + O[t]^Ceiling[Dickm$ex-n])/(t + n - 1/2), t];
(Normal[u] /. t -> 1/2) + (Normal[v] /. t -> -1/2) - v];
Dickman$f[y_?NumberQ] := Which[y<=0, 1, y>Dickm$ma, N[1, y Log[10., y]]-1,
True, Normal[Dickm$dv[Ceiling[y]]] /. t -> y - Ceiling[y] + 1/2];
DickmanF[x_?NumberQ] = If [x<=0, 0, Dickman$f[1/x]];
I get Timing[DickmanF[1/20]] = {6.5 Second, 2.4617828288*10^-29}
Note that the power series expansions are cached, so the next
evaluations are reasonably fast, which is good news for Plot.
The above technique works better than any other I tried; it appears
to be right down to the last digit displayed, even when Dickm$pr is
increased. There are still some problems:
1) accuracy of coeficients in high powers of the series is overkill
2) a little too many coefficients are kept: less than Dickm$ex-n
terms are actually necessary
3) DickmanF gets slow and memory hungry when Dickm$pr rise
4) DickmanF is undetermined below some limit
5) precision is not adaptative:
N[DickmanF[1/20],30] does not work automagically
6) lack of symbolic handling, when we could have
Simplify[DickmanF[x], x>=1/3 && x<=1/2] equal to 1+Log[x]
D[DickmanF[x], x] equal to -DickmanF[x/(1-x)] when x>1
Regarding 3 AND 4, there is some hope, as a more direct computation is
conceivable: reportedly, if s if the positive solution to
E^s-1 == s/x then DickmanF[x] is close to
E^(EulerGamma - s/x + Integrate[(E^t-1)/t, {t,0,s}]) / Sqrt[2 Pi/x]
[Reference: N. de Bruijn, On the number of positive integers <= x and
free of prime factors >= y, Indagationes Mathematicae, vol. 13,
pp. 50-60, 1951 and vol. 28, pp. 236-247, 1966; could not get at it]
It turns out that de Bruijn's approximation can be worked out
directly using buit-in Mathematica 3 and 4 functions:
Dickman$dB[x_] = -E^(1+ExpIntegralEi[-x-ProductLog[-1,-x*E^-x]]+
ProductLog[-1,-x*E^-x]/x)*Sqrt[x/2/Pi]/(x+ProductLog[-1,-x*E^-x]);
But this approximation appears to improve slowly when x gets smaller
(-13.5% for x=1/8, -11.4% for x=1/16, -9.8% for x=1/32); could this
be improved ?
Francois Grieu