Computing Knopfmacher's limit, or: My first foray into computational mathematics, reprise

Daniel Lichtblau
danl@wolfram.com
Wolfram Research, Inc.
100 Trade Center Dr.
Champaign IL 61820

ACA 2009
Montreal, Canada
Session: Applications of Math Software to Mathematical Research

ABSTRACT

I will discuss a problem I encountered over a decade ago, and worked on via internet with someone I (alas) never met. It involves a mix of number theory, real analysis, hard-core computation, and some slightly perplexing results.

In brief, we begin with a function expressed as a certain infinite product; Arnold Knopfmacher encountered it in an attempt to approximate the number of irreducible factors of univariate polynomials over Galois fields and raised the question of how to obtain a certain limit to this function. We derive and execute an effective algorithm for the task at hand. We'll also indicate why the most "obvious" approach does not work well in practice, or at all in theory.

The problem

As posed by Arnold Knopfmacher to the Usenet group comp.soft-sys.math.mathematica in January 1999

Let d(k) denote the smallest prime factor of k.

Define m(k)=k - k/d(k)

Define

p(x)=(((∏_(k=2)^∞)(1-((x^(m(k)))/(k+1))))/(1-x))

We wish to compute numerically, to at least eight decimal places, the following limit.

(lim_(x→(SuperMinus[1]))) p(x)

That is to say, we compute the limit as x approaches 1 from the left (lesser) side.

Why do we care?

There are similar formulas in a paper from 1995 by Knopfmacher and Warlimont, analyzing  probabilities related to numbers of irreducible factors of distinct degrees  in univariate polynomials over Galois fields

It's an interesting computation

My history with this problem

I worked on it off and on for several days. Then someone else reading the forum contributed a similar result, but much more precise. His name was Jürgen Tischer, a math department faculty member of Universidad del Valle, Columbia. We corresponded a bit over a period of months, and I wrote up the results. I lost contact with him a year or two later. I had wanted to invite him to this conference. After getting nowhere with an internet search for a current contact address, I learned he had passed away in January of 2008. I felt it fitting to talk a bit about this problem, since it used ideas of his and also was one of my first forays into computational mathematics.

Candidate results

There were four responses, with two (mine and Tischer's) giving roughly the same results. The proposed values were

1.3397

2 (exactly)

2.292

So which is correct?

More important: How do we even know the limit exists?

Easy to show...

A lim sup and lim inf  both are readily computed.

A lim inf  is given by (e^γ)(the exponential of the Euler gamma constant). This is around 1.871, so...1.3397 will exit stage left.

Definitions in Mathematica

d[k_] := Divisors[k][[2]]
m[k_]:=k-(k/(d[k]))
p[x_]:=(((∏_(k=2)^∞)(1-((x^(m[k]))/(k+1))))/(1-x))

General idea: Start with

(((∏_(k=2)^∞)(1-((x^(⌈(k/2)⌉))/(k+1))))/(1-x))<p[x]<(((∏_(k=2)^∞)(1-((x^(k-1))/(k+1))))/(1-x))

Now take logs to get summations. Expand logs at 1 as power series, obtaining double summations. Switch order of summation (requires justification), and we find that log of a lim sup is γ + log(2). Finding a lim inf  is similar though a bit more work.

What has changed in the past decade?

Ten years ago this computing took manual intervention. I had to do things like split sums, and do further contortions to take limits. Today some can be done directly. Here is one such that arose in the process.

Limit[(∑_(k=3)^M)(∑_(j=2)^∞)(1/(j (k^j))),M->∞]

A start at approximating the actual limit

Truncate the series for the logarithm

Evaluate using exact or  high precision arithmetic at x=1

Exponentiate

Log[p[x]]=-Log[1-x]+(∑_(2|k))Log[1-((x^(k/2))/(k+1))]+(∑_(2∤k,3|k))Log[1-((x^(2 k/3))/(k+1))]+...=-(∑_(k=1)^∞)((x^k)/k)+(∑_(2|k))((x^(k/2))/(k+1))+(∑_(2|k))(∑_(j=2)^∞)((x^(j k/2))/(j (((k+1))^j)))+(∑_(2∤k,3|k))((x^(2 k/3))/(k+1))+(∑_(2∤k,3|k))(∑_(j=2)^∞)((x^((2/3)j k))/(j (((k+1))^j))) +...

(The series that get truncated are the summations in j). From this tactic I was able to get 2.292 (so, as you may have guessed, exact 2 was not the correct result either).

Troublesome aspects

There are serious problems with this approach.

Difficult to get good precision

(Related, but more serious) It is quite difficult to bound the error. Indeed, it is not easy to show we have convergence.

Tischer's idea

Figure out exact forms for some of the infinite sums, so as to avoid truncation. In parts we cannot compute exactly, show that error is much better than what we have from above approach.

Start by writing log(p(x)) (after a bit of algebra) as

x+((x^2)/2)+(∑_(k=2)^∞)(((x^(k+1))/(k+1))-((x^(m[k]))/(k+1))) -(∑_(k=2)^∞)(∑_(j=2)^∞)((1/j))(((((x^(m[k]))/(k+1))))^j)

Proposition: This approaches γ+log(2)+(∑_(k=2)^∞)(((x^(k+1))/(k+1))-((x^(m(k)))/(k+1))) as x->1 from below.

Sketch of proof

Clearly we only need focus on the double summation. Switch summation order and use

(∑_(k=2)^∞)(((((x^k)/(k+1))))^j)≤(∑_(k=2)^∞)(((((x^(m[k]))/(k+1))))^j)≤(∑_(k=2)^∞)(((((x^(k/2))/(k+1))))^j)

Middle is squeezed to

(∑_(k=2)^∞)((((1/(k+1))))^j)

So we can find:

(∑_(j=2)^∞)((1/j))(∑_(k=2)^∞)((((1/(k+1))))^j)

Several steps require justification! We interchanged a summation order, then a sum with a limit...

That remaining summation

We now need to estimate the remaining part. We split by smallest divisors.

(∑_(k=2)^∞)(((x^(k+1))/(k+1))-((x^(m[k]))/(k+1)))=(∑_(d[k]=2)^∞)(((x^(k+1))/(k+1))-((x^(((k/2))))/(k+1))) +(∑_(d[k]=3)^∞)(((x^(k+1))/(k+1))-((x^((((2 k)/3))))/(k+1)))+(∑_(d[k]=5)^∞)(((x^(k+1))/(k+1))-((x^((((4 k)/5))))/(k+1))) + ...

The reordering is fine: for 0<x<1 each term is negative so we can do this.

Now we need to compute

(∑_(d[k]=Prime[j])^∞)(((x^(k+1))/(k+1))-((x^(((((Prime[j]-1) k)/(Prime[j])))))/(k+1)))

Remaining summation...

We need some functions.

q[j_]:=(∏_(k=1)^j)Prime[k]
r[j_]:=(∏_(k=1)^(j-1))(Prime[k]-1)
frac[j_]:=((r[j])/(q[j]))

The terms k for which prime(j) is the smallest divisor larger than 1 fall into finitely many congruence sets. For example, when the prime in question is 5, the applicable values for k are {5,25,35,55,65,85,95,...}. This may be partitioned as {5,35,65,95,...} and {25,55,85,...}. In each case, the step size is 30.

Remaining summation...

In general we have the following lemmas.

The step size of congruence classes for prime(j) is q(j) as defined above

The number of congruence classes is r(j)

When we partition in this way, the limit for each subsum depends only on the prime and is independent of congruence class

Upshot:

(lim_(x->(SuperMinus[1])))(∑_(d[k]=Prime[j])^∞)(((x^(k+1))/(k+1))-((x^(((((Prime[j]-1) k)/(Prime[j])))))/(k+1))) =
 frac[j]Log[1-(1/(Prime[j]))]

Remaining summation...

(lim_(x->(SuperMinus[1])))(∑_(k=2)^∞)(((x^(k+1))/(k+1))-((x^(m[k]))/(k+1)))=(∑_(j=1)^M)frac[j] Log[1-(1/(Prime[j]))]+(lim_(x->(SuperMinus[1])))(∑_(d[k]≥Prime[M+1])^∞)(((x^(k+1))/(k+1))-((x^(((((d[k]-1) k)/(d[k])))))/(k+1)))

We can readily bound that tail sum.

tailsumbnd= Sum[((x^(k+1))/(k+1))-((x^(((Prime[M]-1) k)/(Prime[M])))/(k+1)),{k,0,∞}]

Limit[tailsumbnd,x->1,Direction->1]

Computing our estimate

We can now put all this together.

estimate[n_]:=(∑_(j=1)^n)frac[j] Log[1-(1/(Prime[j]))]
error[n_]:=(∑_(j=n+1)^∞)frac[j] Log[1-(1/(Prime[j]))]

numestimate[n_,prec_]:= Module[{sum=-1/2*Log[2],frac=1/2,p1=Prime[1],p2},
Do[
p2=N[Prime[j],prec];
frac *= (((p1-1)/p2));
sum += frac*Log[1-(1/p2)];
p1=p2,
{j,2,n}];
{frac*p2,sum}]

Timing[{mfact,est} = numestimate[2*10^8+1,25]]

Computing...

Upper bound on error. Set n=2 10^8 .

mfact(∑_(j=n+1)^∞)(1/(Prime[j]))|Log[1-(1/(Prime[j]))]|<mfact(1+(1/(Prime[n])))(∑_(j=n+1)^∞)(1/((j^2) ((Log[j])^2)))<57/2000 (∫_n^∞)(1/(((j log[j]))^2))dj

errormax=(57/2000) Integrate[1/(j Log[j])^2,{j,2*10^8,Infinity}]
N[errormax]

Computing...

Finally we get our estimate and error bound.

p1 = Exp[EulerGamma+Log[2]+est]

errbound = N[Exp[errormax]-1]*p1

Further items of interest

We can investigate the error term of the "naive" summation approach by looking at the series of the log of the product.

logp[x_,n_]:=Normal[
Series[-Log[1-x]+(∑_(k=2)^(2*n))Log[1-((x^(m[k]))/(k+1))], {x,0,n}]]

We can use this to get the signs of the terms. I show them in a run-length form.

They seem to alternate, with sporadic runs of three positive terms. Strange...

Further items...

But stranger is the magnitudes of these coefficients. They are not even bounded by  O((1/n)).

Can show:

They are bounded by O(((log [log[ n]])/n))

This bound is tight (we can show that infinitely many coefficients will approach it closely).

The ones that approach closely have interesting factorization patterns (which is why they approach closely).

Figuring this out was more math than computation (our jobs are not going to the machines just yet).

Upshot: a naive summation will clearly give very poor convergence.

Some open problems

Understand the sign patterns of error approximants.

Find a more efficient way to compute the estimate to high precision.

Find an exact closed form for the limit.

Automate more of the symbolic analysis: some still requires manual intervention.

Determine whether the error bound/estimate is tight. If not, improve it (this would be a "cheap" way of getting more digits).

References

A. Knopfmacher (1999a). Usenet news group post to comp.soft-sys.math.mathematica.
http://library.wolfram.com/mathgroup/archive/1999/Jan/msg00023.html

A. Knopfmacher and R. Warlimont (1995). Distinct degree factorizations for  polynomials over a finite field. Transactions of the American Mathematical  Society 347(6):2235-2243.

DL (1999). The evaluation of Knopfmacher's curious limit. Manuscript.

J. Tischer (1999). Usenet news group post to comp.soft-sys.math.mathematica.
http://library.wolfram.com/mathgroup/archive/1999/Jan/msg00355.html