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