From: "Peter L. Montgomery" Subject: Re: Dickman's function Date: Wed, 13 Sep 2000 07:34:19 GMT Newsgroups: sci.math.num-analysis,sci.crypt Summary: [missing] In article Francois Grieu writes: >I'm trying to find or devise simple C code to compute Dickman's >function. For non-negative reals a, this function gives the proportion >of integers N such that the highest prime factor of N is less that N^a. >It verifies: > /1 >F[a] = 1 for a>=1 F[a] = 1 - | F[t/(1-t)]/t dt for 0<=a<=1 > /a >Reference: Donald E. Knuth, The Art of Computer Programming, volume 2, >section 4.5.4, p367 (2nd ed) or p383 (3rd ed). >Online ref: >Things I tried so far are very imperfect, but here they are. It is handy >to define the auxiliary function f[b] = F[1/b] and we get > /b >f[b] = 1 for 0<=b<=1 f[b] = f[c] - | f[t-1]/t dt for b>=c>=1 > /c Andrew Granville spoke twice at a recent Algorithmic Number Theory workshop in Berkeley, CA. Go to http://www.msri.org . Under Publications, select Lectures on Streaming Video. Choose the Fall 2000 lectures. He suggests using u u f(u) = int f(t) dt (u >= 1) u-1 [The two sides are both zero at u = 1. Their derivatives agree when u >= 1, so they must be equal everywhere.] -- E = m c^2. Einstein = Man of the Century. Why the squaring? Peter-Lawrence.Montgomery@cwi.nl Home: San Rafael, California Microsoft Research and CWI ============================================================================== From: George Marsaglia Subject: Re: Dickman's function Date: Wed, 13 Sep 2000 06:58:40 -0400 Newsgroups: sci.math.num-analysis,sci.crypt Summary: [missing] Francois Grieu wrote: > I'm trying to find or devise simple C code to compute Dickman's > function. For non-negative reals a, this function gives the proportion > of integers N such that the highest prime factor of N is less that N^a. > It verifies: > > /1 > F[a] = 1 for a>=1 F[a] = 1 - | F[t/(1-t)]/t dt for 0<=a<=1 > /a > > Reference: Donald E. Knuth, The Art of Computer Programming, volume 2, > section 4.5.4, p367 (2nd ed) or p383 (3rd ed). > Online ref: > > Things I tried so far are very imperfect, but here they are. It is handy > to define the auxiliary function f[b] = F[1/b] and we get > > /b > f[b] = 1 for 0<=b<=1 f[b] = f[c] - | f[t-1]/t dt for b>=c>=1 > /c > > In the above, c could be any real but using c = floor[b] which gives a > simple recurence. For exemple on the segment 1<=b<=2 we get > > /b > f[b] = 1 - | 1/t dt thus f[b] = 1 - Log[b] for 1<=b<=2 > /1 > > More symbolic integration does not help much: on segment 2<=b<=3 we need > the PolyLog function, and farther the expression grows madly. > > Recursive numerical integration with a general method works for b up to > say 5, then gets terribly slow. A faster variant is to build sucessive > polynomial approximations at points j+1/2 for integers j>=0: > > fj[x] = f[x + (j+1/2)] = sum(Aij * x^i) + O(x^n) > 0<=i > Convergence is reasonably quick. However any rounding error or > truncation made in computing a segment affects the next ones as an > absolute error, and since the function decreases quickly to 0, we get > totaly innacurate results for b>12 when using 64 bits IEEE arithmetic. > > Anyone has a method working for b up to say 40 without using arbitrary > precision arithmetic or too many magic values ? > 19.25 7.2851e-28 | 19.5 2.3646e-28 | 19.75 7.6446e-29 | 20 2.4618e-29 I developed a method some years ago, ~ 1970. It was described in ``Numerical Solution of Some Classical Differential-Difference Equations,'' (1989) , Mathematics of Computation, 53, 191-201, but it requires extended precision arithmetic. Many reported results, based on numerical integration and not using extended precision, were found to be far off the true values. George Marsaglia ============================================================================== From: bach@jalapeno.cs.wisc.edu (Eric Bach) Subject: Re: Dickman's function Date: 13 Sep 2000 18:51:35 GMT Newsgroups: sci.math.num-analysis,sci.crypt In article <39BF5DDF.372D3F1B@stat.fsu.edu>, George Marsaglia wrote: > > >Francois Grieu wrote: > >> I'm trying to find or devise simple C code to compute Dickman's >> function. Patterson and Rumsey devised a method based on piecewise approximation by analytic functions. I don't think they ever published it but it is described in E. Bach and R. Peralta, Asymptotic Semismoothness Probabilities, Mathematics of Computation 65 (1996), 1701-1715. We give references to other methods too. --Eric Bach ============================================================================== From: djb@cr.yp.to (D. J. Bernstein) Subject: Re: Dickman's function Date: 13 Sep 2000 07:10:27 GMT Newsgroups: sci.math.num-analysis,sci.crypt Francois Grieu wrote: > I'm trying to find or devise simple C code to compute Dickman's > function. There's a rho() function in psibound: http://cr.yp.to/psibound.html ---Dan ============================================================================== From: djb@cr.yp.to (D. J. Bernstein) Subject: Re: Dickman's function Date: 14 Sep 2000 20:35:00 GMT Newsgroups: sci.math.num-analysis,sci.crypt Francois Grieu wrote: > Somehow this code > circumvents the error propagation problem encountered in my methods, Expand rho as a power series to the right of 1, to the right of 3/2, to the right of 2, etc. Use the integral to evaluate the constant term of each power series. Use the differential equation to evaluate the next 20 or 30 terms. This is numerically stable. ---Dan