5  Lost in Hyperspace

The best algorithm we found for this (which is far from the best algorithm that exists!) is the differential evolution method, a type of genetic algorithm included in Mathematica's NMinimize function. Of course, one must first get an efficient way to get the error for a given cubic. Several observations help get that quickly, the most important of which is that the maximum error occurs on the upper half of the unit circle.

The code that follows (for version 4.2 of Mathematica) produces 15 digits (but with no proof that they are correct). Note that the problem did not specify that p(z) was to have real coefficients. The Eindhoven team's solution (link given earlier) provides a proof. We did find another method that reduced the problem from a 4-dimensional optimization to a 1-dimensional optimization of the results of a 3-dimensional optimization, and that technique was good enough to get us 15 digits.
  
But now I see that experts (Laurie, Bornemann, Jansen & Boersma) use the method of Chebyshev approximation, which is efficient and provides a proof of correctness. Jansen and Boersma obtained 100 digits (and their method yields 500 with no problem), which are:

0.21433523459045963946152640018474939611340728778951670807393498597098767607405459816947988478883053305903475633.

Needs["NumericalMath`NMinimize`"] error = Compile[{a, b, c, d, θ}, Abs[a ᢺ ... ion, ,, PrecisionGoal16, ,, AccuracyGoal16, ,, MaxIterations100000}], ]}]

0.214335

We tried to get more digits. The best cubic we found is

      FormBox[RowBox[{RowBox[{RowBox[{-,  , 0.60334322047992600}], z^3}], +, RowBox[{0.625212, z^2}], +, RowBox[{1.0197618528346884, z}], +, 0.00554195}], TraditionalForm],

for which the maximum error is 0.21433 52345 90459 64277. It seemed likely that 16 or 17 of these digits are correct. It turns out that 16 are; 17 if one rounds.

An earlier dimenson-reduction approach of ours is interesting and did get us 15 digits. First one observes that, when one is close, the error curve has local maxima at 1.40, 2.26, and π, where the first two are approximate and the last is exact. Moreover, the error at π is a - b + c - d. So we can specify e, the error at π, and look at the family of cubics of the form a z^3 + b z^2 + c z + a - b + c + e. We can try to find the cubic in this 3-dimensional family for which the overall max-error is least. Then it is simply a matter of finding the value of e for which this abc-minimum value is least, and FindMinimum is up to the task on this 3-dimensional problem (and then on the 1-dimensional problem).

This approach can be viewed as a coordinate transformation of the 4-dimensional problem. Simply using FindMinimum on the raw 4-dimensional problem did not work well. Trying to lower the dimensions directly (fixing d, say, and optimizing a, b, c) does not work. A possible reason for this difference is that as we change e a little bit the overall shape of the curve does not change much; this is analogous to using Bézier polynomials as a design tool. This means that the location of the other two bumps does not change much, and the seed for that search is good. We store the locations found by the abc problem and use them for the next value of e, thus adding efficiency to the 3-dimensional search.

The solution by Boersma and Jansen of the Eindhoven team is noteworthy, as they found a lemma of Vidensky in a 1968 book by Smirnov and Lebedev [Functions of a Complex Variable, Constructive Theory, p. 450-452] that reduces the problem in quite a simple way to six simultaneous equations that can be solved by standard root-finding. Here is their code, which easily yields 500 digits.

w[a_, b_, c_, d_, z_] := 1/Gamma[z] - a - b z - c z^2 - d z^3 ; det[u_, z_] := Det[(   ... }], }}]}], }}], ,,  , , WorkingPrecision25, ,,  , AccuracyGoal15}], ]}]}]

0.214335234590459639480708


Created by Mathematica  (June 27, 2004)