(************** Content-type: application/mathematica **************
Mathematica-Compatible Notebook
This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.
To get the notebook into a Mathematica-compatible application, do
one of the following:
* Save the data starting with the line of stars above into a file
with a name ending in .nb, then open the file inside the
application;
* Copy the data starting with the line of stars above to the
clipboard, then use the Paste menu command inside the application.
Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode. Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).
NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.
For more information on notebooks and Mathematica-compatible
applications, contact Wolfram Research:
web: http://www.wolfram.com
email: info@wolfram.com
phone: +1-217-398-0700 (U.S.)
Notebook reader applications are available free of charge from
Wolfram Research.
*******************************************************************)
(*CacheID: 232*)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 35829, 1220]*)
(*NotebookOutlinePosition[ 36538, 1245]*)
(* CellTagsIndexPosition[ 36494, 1241]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell["\<\
An Improved Algorithm for Diophantine Equations in One \
Variable\
\>", "Title"],
Cell["Adam Strzebonski", "Author"],
Cell["\<\
Wolfram Reseach Inc, 100 Trade Centre Dr, Champaign, IL 61820, \
USA\
\>", "Address"],
Cell[TextData[{
"We present a new algorithm for computing integer roots of univariate \
polynomials. For a polynomial ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ](t)\)]],
" we can find the set of its integer roots in a time polynomial in the size \
of the\nsparse encoding of",
Cell[BoxData[
\(TraditionalForm\`f\)]],
". An algorithm to do this was given by F. Cucker, P. Koiran, S. Smale in \
1999. Their paper introduces a polynomial time method for computing sign of \
",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" at integral points, and it finds integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" by isolating the real roots of all polynomials in the sparse derivative \
sequence of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", up to unit-length intervals. We propose to isolate the roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" using a sparse variant of Fourier`s theorem, and show that our algorithm \
requires a smaller number of polynomial sign evaluations. We present an \
empirical comparison of our implementations of the original and of the \
improved algorithm, and of an algorithm using modular root\nfinding and \
Hensel lifting, suitable for dense polynomials."
}], "Abstract"],
Cell[CellGroupData[{
Cell["Introduction", "Section"],
Cell["Let", "Text"],
Cell[BoxData[
\(f = \(a\_1\) t\^e\_1 + \[Ellipsis] + \(a\_k\)
t\^e\_k\)], "NumberedEquation"],
Cell[TextData[{
"where ",
Cell[BoxData[
\(TraditionalForm\`e\_1 > \[Ellipsis] > e\_k\)]],
", and ",
Cell[BoxData[
\(TraditionalForm\`a\_i \[Element] \[DoubleStruckCapitalZ]\\{0}\)]],
", and let us denote"
}], "Text"],
Cell[BoxData[
\(minexp \((f)\) := e\_k\)], "NumberedEquation"],
Cell[TextData[{
"The sparse derivative sequence of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" is the sequence ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
"defined by"
}], "Text"],
Cell[BoxData[{
\(f\_1 := f/t\^\(minexp \((f)\)\)\), "\[IndentingNewLine]",
\(f\_\(i + 1\) := f\_i'/t\^\(minexp \((f\_i')\)\)\)}], "NumberedEquation"],
Cell[TextData[{
"The algorithm defined in [CKS] (F.Cucker, P.Koiran, S.Smale, ``A \
Polynomial Time Algorithm for Diophantine Equations in One Variable'', \
J.Symbolic Comp., 27 (1999), 21-29), finds integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" by isolating the real roots of polynomials of the sparse derivative \
sequence of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", taken in the reverse order, up to intervals ",
Cell[BoxData[
\(TraditionalForm\`\((u, \ u + 1)\)\)]],
" or ",
Cell[BoxData[
\(TraditionalForm\`\([u, u]\)\)]],
", with ",
Cell[BoxData[
\(TraditionalForm\`u \[Element] \[DoubleStruckCapitalZ]\)]],
". The method is based on the fact that, by Rolle's theorem, a polynomial \
can have at most one root between two subsequent roots of its derivative. In \
the worst case the algorithm may require isolating ",
Cell[BoxData[
\(TraditionalForm\`k\^2\)]],
" roots."
}], "Text"],
Cell[TextData[{
"Our algorithm isolates the roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" directly, using a ``sparse variant'' of Fourier's theorem proven in the \
next section. It requires isolating at most ",
Cell[BoxData[
\(TraditionalForm\`2\ k - 2\)]],
" roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" or polynomials in its sparse derivative sequence. "
}], "Text"],
Cell[TextData[{
"We will also present an algorithm finding roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" modulo a prime and using Hensel lifting to obtain the integer roots. \
Since ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" may have ",
Cell[BoxData[
\(TraditionalForm\`e\_1\)]],
" roots modulo a prime ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
", even if ",
Cell[BoxData[
\(TraditionalForm\`k = 2\)]],
", for instance if ",
Cell[BoxData[
\(TraditionalForm\`f = t\^p - t\)]],
", this algorithm is exponential in the size of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
". However, our experiments suggest that for dense polynomials this method \
is the fastest."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Root Isolation Theorem", "Section"],
Cell[TextData[{
"Let ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], \(f\_k\) :
I\[LongRightArrow]\[DoubleStruckCapitalR]\)]],
" be a sequence of differentiable functions defined in an open interval ",
Cell[BoxData[
\(TraditionalForm\`I \[Subset] \[DoubleStruckCapitalR]\)]],
", and such\nthat for all ",
Cell[BoxData[
\(TraditionalForm\`1 \[LessEqual] i \[LessEqual] k - 1\)]],
", and for all ",
Cell[BoxData[
\(TraditionalForm\`x \[Element] I\)]]
}], "Text"],
Cell[BoxData[
\(sign \((\(f\_\(i + 1\)\) \((x)\))\) =
sign \((f\_i' \((x)\))\)\)], "NumberedEquation"],
Cell["and", "Text"],
Cell[BoxData[
\(\(f\_k\) \((x)\) \[NotEqual] 0\)], "NumberedEquation"],
Cell[TextData[{
"For ",
Cell[BoxData[
\(TraditionalForm\`x \[Element] I\)]],
", and ",
Cell[BoxData[
\(TraditionalForm\`1 \[LessEqual] p \[LessEqual] q \[LessEqual] k\)]],
" let ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(p, q\)\)(x)\)]],
" denote the number of sign changes in the sequence ",
Cell[BoxData[
\(TraditionalForm\`\(f\_p\)(x), \[Ellipsis], \(f\_q\)(x)\)]],
" with terms equal to zero removed, and let ",
Cell[BoxData[
\(TraditionalForm\`sgc(x) := \(sgc\_\(1, k\)\)(x)\)]],
"."
}], "Text"],
Cell[TextData[{
StyleBox["Theorem: ",
FontWeight->"Bold"],
StyleBox["For any ",
FontVariations->{"CompatibilityType"->0}],
Cell[BoxData[
\(TraditionalForm\`a \[Element] I\)]]
}], "Text"],
Cell[BoxData[{
\(lim\+\(x \[Rule] \(a\^+\)\)\ sgc \((x)\) =
sgc \((a)\)\), "\[IndentingNewLine]",
\(lim\+\(x \[Rule] \(a\^-\)\)\ sgc \((x)\) =
sgc \((a)\) + r + 2 s\)}], "NumberedEquation"],
Cell[TextData[{
"where ",
Cell[BoxData[
\(TraditionalForm\`r \[Element] \(\[DoubleStruckCapitalZ]\_+\)\)]],
", ",
Cell[BoxData[
\(TraditionalForm\`s \[Element] \(\[DoubleStruckCapitalZ]\_+\) \[Union] \
{0}\)]],
", and"
}], "Text"],
Cell[BoxData[{
\(\(f\_1\) \((a)\) = \(\[Ellipsis] = \(\(f\_r\) \((a)\) =
0\)\)\), "\[IndentingNewLine]",
\(\(f\_\(r + 1\)\) \((a)\) \[NotEqual] 0\)}], "NumberedEquation"],
Cell[TextData[{
"Moreover, ",
Cell[BoxData[
\(TraditionalForm\`s = 0\)]],
", unless there is ",
Cell[BoxData[
\(TraditionalForm\`t > r + 1\)]],
", such that ",
Cell[BoxData[
\(TraditionalForm\`\(f\_t\)(a) = 0\)]],
"."
}], "Text"],
Cell[TextData[{
StyleBox["Corollary:",
FontWeight->"Bold"],
" Let ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
" be the sparse derivative sequence of a polynomial ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", let ",
Cell[BoxData[
\(TraditionalForm\`0 \[LessEqual] a < b\)]],
", and let ",
Cell[BoxData[
\(TraditionalForm\`1 \[LessEqual] p \[LessEqual] q \[LessEqual] k\)]],
", be such that ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(p, k\)\)(a) = \(sgc\_\(p, k\)\)(b)\)]],
". Then"
}], "Text"],
Cell[TextData[{
"(1) ",
Cell[BoxData[
\(TraditionalForm\`f\_p\)]],
" has a constant nonzero sign on ",
Cell[BoxData[
\(TraditionalForm\`\((a, b]\)\)]]
}], "Text"],
Cell[TextData[{
"(2)",
" ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(p, k\)\)(c) = \(sgc\_\(p, k\)\)(a)\)]],
" for any ",
Cell[BoxData[
\(TraditionalForm\`a \[LessEqual] c \[LessEqual] b\)]]
}], "Text"],
Cell[TextData[{
"(3) ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(1, p\)\)(a) - \(sgc\_\(1, p\)\)(b) =
r + 2 s\)]],
", where ",
Cell[BoxData[
\(TraditionalForm\`r\)]],
" is the number of roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" in ",
Cell[BoxData[
\(TraditionalForm\`\((a, b]\)\)]],
", counted with multiplicities, and ",
Cell[BoxData[
\(TraditionalForm\`s\)]],
" is a nonnegative integer. In particular, if ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(1, p\)\)(a) = \(sgc\_\(1, p\)\)(b)\)]],
" then ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" has no roots in ",
Cell[BoxData[
\(TraditionalForm\`\((a, b]\)\)]],
", and if ",
Cell[BoxData[
\(TraditionalForm\`\(sgc\_\(1, p\)\)(a) - \(sgc\_\(1, p\)\)(b) = 1\)]],
" then ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" has exactly one root in ",
Cell[BoxData[
\(TraditionalForm\`\((a, b]\)\)]],
", and ",
Cell[BoxData[
\(TraditionalForm\`sign(f(a)) \[NotEqual] sign(f(b))\)]],
"."
}], "Text"],
Cell[TextData[{
StyleBox["Corollary:",
FontWeight->"Bold"],
" A univariate polynomial with ",
Cell[BoxData[
\(TraditionalForm\`k > 0\)]],
" nonzero terms has at most ",
Cell[BoxData[
\(TraditionalForm\`k - 1\)]],
"positive roots and at most ",
Cell[BoxData[
\(TraditionalForm\`k - 1\)]],
" negative roots, and hence at most ",
Cell[BoxData[
\(TraditionalForm\`2\ k - 1\)]],
" real roots. For completeness, let us note that the polynomial "
}], "Text"],
Cell[BoxData[
\(f = \[Product]\+\(i = 0\)\%\(k - 1\)\((t\^2 -
i\^2)\)\)], "DisplayFormula"],
Cell[TextData[{
"has exactly ",
Cell[BoxData[
\(TraditionalForm\`k\)]],
" nonzero terms and ",
Cell[BoxData[
\(TraditionalForm\`2 k - 1\)]],
" real roots."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["The Algorithm", "Section"],
Cell[TextData[{
"We present an algorithm computing the positive integer roots of a \
polynomial ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
". To find the negative integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" we compute the positive integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f(\(-t\))\)]],
". To compute the sign of a polynomial we use the algorithm described in \
[CKS]. As a positive integer root bound for ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", we take the smaller of the bound used in [CKS] (absolute value of the \
trailing coefficient of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
"), and the positive root bound described in J. R. Johnson, ``Algorithms \
for Polynomial Real Root Isolation'', in B. Caviness, J. Johnson (eds.), \
Quantifier Elimination and Cylindrical Algebraic Decomposition, \
Springer-Verlag 1998, 269-299. We also use the following two\nalgorithms."
}], "Text"],
Cell[TextData[{
StyleBox["BinarySearch(f,a,b): ", "Program"],
"for ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`a, b \[Element] \[DoubleStruckCapitalZ]\)]],
", such that ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" has exactly one root in ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
", and ",
Cell[BoxData[
\(TraditionalForm\`\(f(a)\) \(f(b)\) < 0\)]],
", finds out whether the root is an integer, and if yes returns it. The \
algorithm uses binary search\nwith integer subdivision points, and stops when \
it finds an integer root, or isolates the root to an interval of length 1."
}], "Text"],
Cell[TextData[{
StyleBox["ExhaustiveSearch(f,a,b,r):", "Program"],
" for ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`a, b \[Element] \[DoubleStruckCapitalZ]\)]],
", finds the integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" in ",
Cell[BoxData[
\(TraditionalForm\`\([a, b]\)\)]],
", assuming that there is at most ",
Cell[BoxData[
\(TraditionalForm\`r\)]],
" of them. The algorithm computes the sign of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" at subsequent integers in ",
Cell[BoxData[
\(TraditionalForm\`\([a, b]\)\)]],
", and stops when it finds ",
Cell[BoxData[
\(TraditionalForm\`r\)]],
" roots, or when it has checked all integers in ",
Cell[BoxData[
\(TraditionalForm\`\([a, b]\)\)]],
"."
}], "Text"],
Cell[TextData[{
StyleBox["Algorithm 1:",
FontWeight->"Bold"],
" Given ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
" returns the list of positive integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
"."
}], "Text"],
Cell[TextData[{
"(1) Compute the sparse derivative sequence ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
" of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", a positive integer root bound ",
Cell[BoxData[
\(TraditionalForm\`rb\)]],
" for ",
Cell[BoxData[
\(TraditionalForm\`f\_1\)]],
", and set ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
" to an empty list. While ",
Cell[BoxData[
\(TraditionalForm\`\(f\_1\)(rb) = 0\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`rb > 0\)]],
" add ",
Cell[BoxData[
\(TraditionalForm\`rb\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
" and decrement ",
Cell[BoxData[
\(TraditionalForm\`rb\)]],
". Put interval ",
Cell[BoxData[
\(TraditionalForm\`\((0, rb)\)\)]],
" on ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
"."
}], "Text"],
Cell[TextData[{
"(2) If ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
" is empty, return ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
", else take an interval ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
" off ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
"."
}], "Text"],
Cell[TextData[{
"(3) If ",
Cell[BoxData[
\(TraditionalForm\`sgc(a) = sgc(b)\)]],
" there are no roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" in ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
". Go to step 2."
}], "Text"],
Cell[TextData[{
"(4) If ",
Cell[BoxData[
\(TraditionalForm\`sgc(a) - sgc(b) = 1\)]],
", compute ",
Cell[BoxData[
\(TraditionalForm\`BinarySearch(f\_1, a, b)\)]],
". If it finds an integer root add it to ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
". Go to step 2."
}], "Text"],
Cell[TextData[{
"(5) If ",
Cell[BoxData[
\(TraditionalForm\`b - a \[LessEqual] k\)]],
", compute ",
Cell[BoxData[
\(TraditionalForm\`ExhaustiveSearch(f\_1, a, b, sgc(a) - sgc(b))\)]],
". Add the roots found to ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
". Go to step 2."
}], "Text"],
Cell[TextData[{
"(6) Choose integers ",
Cell[BoxData[
\(TraditionalForm\`a < c = d < b\)]],
". While ",
Cell[BoxData[
\(TraditionalForm\`\(f\_1\)(c) = 0\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`c > a\)]],
" add ",
Cell[BoxData[
\(TraditionalForm\`c\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
" and decrement ",
Cell[BoxData[
\(TraditionalForm\`c\)]],
". While ",
Cell[BoxData[
\(TraditionalForm\`\(f\_1\)(d) = 0\)]],
"\nand ",
Cell[BoxData[
\(TraditionalForm\`d < b\)]],
" add ",
Cell[BoxData[
\(TraditionalForm\`d\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`rootlist\)]],
" and increment ",
Cell[BoxData[
\(TraditionalForm\`d\)]],
". Put intervals ",
Cell[BoxData[
\(TraditionalForm\`\((a, c)\)\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`\((d, b)\)\)]],
" on ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
" and go to step 2."
}], "Text"],
Cell[TextData[{
StyleBox["Remark:",
FontWeight->"Bold"],
" There are several efficiency improvents that can be made to the algorithm \
described above."
}], "Text"],
Cell[TextData[{
"(a) ",
Cell[BoxData[
\(TraditionalForm\`sgc(0)\)]],
" is equal to the number of sign changes in the coefficient list of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
". We can compute it before computing the sparse derivative sequence, and \
if it is 0 or 1, we know that there are no roots, or we can find the only \
root using BinarySearch. Similarly,if ",
Cell[BoxData[
\(TraditionalForm\`rb\)]],
" is small (",
Cell[BoxData[
\(TraditionalForm\`rb < k\)]],
" in our implementation), we use ExhaustiveSearch on ",
Cell[BoxData[
\(TraditionalForm\`\([1, rb]\)\)]],
", without computing the sparse derivative sequence."
}], "Text"],
Cell[TextData[{
"(b) We can construct the loop somewhat differently, so that we check \
whether ",
Cell[BoxData[
\(TraditionalForm\`sgc(a) - sgc(b)\)]],
" is 0 or 1, or whether ",
Cell[BoxData[
\(TraditionalForm\`b - a\)]],
" is small before adding ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
" to ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
". This way ",
Cell[BoxData[
\(TraditionalForm\`intervalstack\)]],
" contains only intervals for which ",
Cell[BoxData[
\(TraditionalForm\`sgc(a) - sgc(b) > 1\)]],
", and hence its height never exceeds ",
Cell[BoxData[
\(TraditionalForm\`\((k - 1)\)/2\)]],
"."
}], "Text"],
Cell[TextData[{
"(c) Keep the computed signs of ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
" at interval's endpoints together with the interval, so that we don't need \
to recompute them."
}], "Text"],
Cell[TextData[{
"(d) For every interval ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
" which we need to subdivide, compute the smallest ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
" satisfying the conditions of Corollary. Then for any point ",
Cell[BoxData[
\(TraditionalForm\`c\)]],
" subdividing ",
Cell[BoxData[
\(TraditionalForm\`\((a, b)\)\)]],
" we only need to compute the signs of ",
Cell[BoxData[
\(TraditionalForm\`\(f\_1\)(c), \[Ellipsis], \(f\_\(p - 1\)\)(c)\)]],
"."
}], "Text"],
Cell[TextData[{
"(e) In point 6 of the algorithm we choose ",
Cell[BoxData[
\(TraditionalForm\`c = floor \((\((a + b)\)/2)\)\)]],
", except that since the changes in value of ",
Cell[BoxData[
\(TraditionalForm\`sgc\)]],
" tend to happen close to zero, for intervals with ",
Cell[BoxData[
\(TraditionalForm\`a = 0\)]],
"we choose ",
Cell[BoxData[
\(TraditionalForm\`c\)]],
" to be closer to zero than ",
Cell[BoxData[
\(TraditionalForm\`b/2\)]],
"."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Complexity Analysis", "Section"],
Cell[TextData[{
"Let us compute the number of polynomial sign evaluations required by our \
algorithm (let us call it FIR, for Fourier Integer Roots), and by the \
algorithm described in [CKS] (let us call it RIR for Rolle Integer Roots). \
Let ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
", let ",
Cell[BoxData[
\(TraditionalForm\`k\)]],
" be the number of nonzero terms in ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", and let ",
Cell[BoxData[
\(TraditionalForm\`M\)]],
" be a bound for integer roots of polynomials in the sparse derivative \
sequence ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
" of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
". For simplicity we take a common root bound for all ",
Cell[BoxData[
\(TraditionalForm\`f\_1, \[Ellipsis], f\_k\)]],
". Since the worst case estimate for size of coefficients of polynomials ",
Cell[BoxData[
\(TraditionalForm\`f\_i\)]],
" grows with ",
Cell[BoxData[
\(TraditionalForm\`i\)]],
", taking a common estimate is ``to the advantage'' of RIR."
}], "Text"],
Cell[TextData[{
"FIR isolates at most ",
Cell[BoxData[
\(TraditionalForm\`2 k - 2\)]],
" changes in value of ",
Cell[BoxData[
\(TraditionalForm\`sgc\)]],
", from an interval of length of at most ",
Cell[BoxData[
\(TraditionalForm\`M\)]],
" to an interval of length\nof at least 1, evaluating signs of at most ",
Cell[BoxData[
\(TraditionalForm\`k\)]],
" polynomials for each interval subdivision. Hence, the total number of \
polynomial sign evaluations is bounded by ",
Cell[BoxData[
\(TraditionalForm\`O(\(k\^2\) \(log(M)\))\)]],
"."
}], "Text"],
Cell[TextData[{
"To isolate roots of ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 1\)\)]],
" RIR needs to evaluate ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 1\)\)]],
" at at most ",
Cell[BoxData[
\(TraditionalForm\`\((m - 1)\)\^2\)]],
" endpoints of intervals isolating roots\nof ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 2\), \[Ellipsis], f\_k\)]],
", and then it needs to isolate at most ",
Cell[BoxData[
\(TraditionalForm\`2 m - 1\)]],
" roots of ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 1\)\)]],
" from an interval of length of at most ",
Cell[BoxData[
\(TraditionalForm\`M\)]],
" to intervals of length 1, evaluating 1 polynomial sign for each \
subdivision. Hence, the number of polynomial sign evaluations necessary for \
isolating roots of ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 1\)\)]],
", assuming the roots of ",
Cell[BoxData[
\(TraditionalForm\`f\_\(k - m + 2\), \[Ellipsis], f\_k\)]],
" have already been isolated, is bounded by ",
Cell[BoxData[
\(TraditionalForm\`O(m\^2 + m\ \(log(M)\))\)]],
". Therefore the total number of polynomial sign evaluations necessary for \
finding\nthe integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" using RIR is bounded by ",
Cell[BoxData[
\(TraditionalForm\`O(k\^3 + \(k\^2\) \(log(M)\))\)]],
"."
}], "Text"],
Cell[TextData[{
"Clearly when ",
Cell[BoxData[
\(TraditionalForm\`k \[GreaterGreater] log(M)\)]],
" FIR has a lower complexity than RIR. However, the complexity computed \
above is a worst case estimate. In practice we quite often get large \
intervals with ",
Cell[BoxData[
\(TraditionalForm\`sgc(a) - sgc(b) = 1\)]],
", in which case in the further subdivision of the interval we need only to \
compute signs of ",
Cell[BoxData[
\(TraditionalForm\`f\_1\)]],
". Also, the use of Remark helps to lower the number of polynomial signs \
that need to be computed with each interval subdivision."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["The Dense Polynomial Case", "Section"],
Cell[TextData[{
"Another possible approach to finding integer roots of a polynomial ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" is to find the roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" modulo a prime ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
", and then use Hensel lifting to find the integer roots. Such an algorithm \
is exponential in the size of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" because the number of roots modulo ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
" may be as large as the degree of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
". Therefore the algorithm is not suitable for sparse polynomials. However, \
it works quite well when ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
" is dense."
}], "Text"],
Cell[TextData[{
StyleBox["Algorithm 2:",
FontWeight->"Bold"],
" Given ",
Cell[BoxData[
\(TraditionalForm\`f \[Element] \[DoubleStruckCapitalZ][t]\\{0}\)]],
" returns the list of integer roots of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
"."
}], "Text"],
Cell[TextData[{
"(1) Compute the square free part ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" of ",
Cell[BoxData[
\(TraditionalForm\`f\)]],
", and an integer root bound ",
Cell[BoxData[
\(TraditionalForm\`rb\)]],
" for ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
"."
}], "Text"],
Cell[TextData[{
"(2) Select a prime ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
" such that ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" is square free modulo ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
"."
}], "Text"],
Cell[TextData[{
"(3) Compute roots ",
Cell[BoxData[
\(TraditionalForm\`r\_1, \[Ellipsis], r\_s\)]],
" of ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" modulo ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
", by computing ",
Cell[BoxData[
\(TraditionalForm\`h = gcd(g, t\^p - t)\)]],
", and factoring it using the algorithm\ndescribed in [Z] R. Zippel, \
``Effective Polynomial Computation'', Kluwer Academic Publishers, \
Boston/Dodrecht/London 1993, section 18.3."
}], "Text"],
Cell[TextData[{
"(4) Lift ",
Cell[BoxData[
\(TraditionalForm\`r\_1, \[Ellipsis], r\_s\)]],
" to roots ",
Cell[BoxData[
\(TraditionalForm\`r\_1', \[Ellipsis], r\_s'\)]],
" of ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" mod ",
Cell[BoxData[
\(TraditionalForm\`p\^\(2\^m\)\)]],
", with ",
Cell[BoxData[
\(TraditionalForm\`p\^\(2\^m\) \[GreaterEqual] 2 rb + 1\)]],
", using the quadratic iteration Hensel lifting described in [Z], section \
16.2."
}], "Text"],
Cell[TextData[{
"(5) Return those of ",
Cell[BoxData[
\(TraditionalForm\`r\_1' - p\^\(2\^m\), \[Ellipsis],
r\_s - \(p\^\(2\^m\)\)', r\_1', \[Ellipsis], r\_s'\)]],
" which are in the interval ",
Cell[BoxData[
\(TraditionalForm\`\([\(-rb\), rb]\)\)]],
" and are roots of ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
"."
}], "Text"],
Cell[TextData[{
"The integer root bound is the greater of the positive and negative integer \
root bounds computed as in Algorithm 1. In our implementation whenever the \
root bound is less than 10 we find the integer roots using a direct search. \
In step 2 we try an increasing sequence of primes, until we find one for \
which ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" is square free modulo ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
". Primes for which ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" is not square free modulo",
Cell[BoxData[
\(TraditionalForm\`p\)]],
"divide the discriminant of ",
Cell[BoxData[
\(TraditionalForm\`g\)]],
" therefore there is only a finite number of them."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Experimental Results", "Section"],
Cell[TextData[{
"In this section we present an empirical comparison of our implementations \
of Algorithm 1 (FIR, for Fourier Integer Roots), Algorithm 2 (HIR, for Hensel \
Integer Roots), and of the algorithm described in [CKS] (RIR, for Rolle \
Integer Roots). Both FIR and RIR use the refinement proposed in section 4 of \
[CKS]. All three algorithms have been\nimplemented as a part of ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" kernel. The computations have been done on a 700 MHz Pentium III laptop \
computer with 256 MB of RAM. Time is given in seconds, columns marked st\\# \
show the number of polynomial sign tests used in FIR and RIR. The results \
given are averages for sets of 10 randomly generated polynomials."
}], "Text"],
Cell["\<\
We have used the following three types of randomly generated \
polynomials.\
\>", "Text"],
Cell[CellGroupData[{
Cell[TextData[{
"Random sparse polynomials of degree ",
Cell[BoxData[
\(TraditionalForm\`n\)]],
" with ",
Cell[BoxData[
\(TraditionalForm\`k\)]],
" integer roots"
}], "Subsection"],
Cell[BoxData[
\(SP \((n, k,
p)\) := \((\(a\_1\) \(t\^\(n - k\)\) \(\(+\[Sum]\)\+\(i =
2\)\%\(p - 1\)\) \(a\_i\) t\^e\_i +
a\_p)\) \(\[Product]\+\(j = 1\)\%k\((t -
r\_j)\)\)\)], "DisplayFormula"],
Cell[TextData[{
"where ",
Cell[BoxData[
\(TraditionalForm\`a\_i\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`r\_j\)]],
" are randomly generated 100-bit integers, and ",
Cell[BoxData[
\(TraditionalForm\`r\_j\)]],
" are all distinct, and ",
Cell[BoxData[
\(TraditionalForm\`0 < e\_i < n - k\)]],
" are randomly generated distinct exponents. "
}], "Text"],
Cell["\<\
For all random sparse polynomials we tried FIR was the fastest \
algorithm, except for polynomials with the highest density (29%and 43%), \
where HIR was faster. The speed advantage of FIR over RIR grows with the \
number of terms. For very high degree polynomials we see that the number of \
polynomial sign tests performed both by FIR and RIR decreases. This is a \
result of using the refinement proposed in section 4 of [CKS].\
\>", "Text"],
Cell[BoxData[GridBox[{
{"Polynomial", "Deg", "Terms", "Roots", \(Time\[IndentingNewLine]
HIR\), \(Time\[IndentingNewLine]
FIR\), \(Time\[IndentingNewLine]
RIR\), \(st #\[IndentingNewLine]
FIR\), \(st #\[IndentingNewLine]
RIR\)},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "19.9", "1", "1.1", "0.005",
"0.005", "110", "317"},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "29.2", "2", "1.11", "0.011",
"0.085", "274", "3616"},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "49.2", "4", "1.14", "0.055",
"0.96", "799", "27299"},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "89.3", "8", "1.19", "0.22",
"7.21", "2198", "122415"},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "157.1", "16", "1.31", "0.88",
"65.6", "5196", "575958"},
{\(SP \((10\^3, 1, 10)\)\), \(10\^3\), "294.9", "32", "1.64",
"4.72", \(\(>\)\(300\)\), "12825", "?"},
{\(SP \((10\^2, 4, 10)\)\), \(10\^2\), "43.1", "4", "0.028", "0.034",
"0.68", "722", "22672"},
{\(SP \((10\^4, 4, 10)\)\), \(10\^4\), "50", "4", "40.4", "0.33",
"0.38", "456", "2874"},
{\(SP \((10\^5, 4, 10)\)\), \(10\^5\), "50", "4", \(\(>\)\(300\)\),
"5.61", "5.65", "418", "1031"},
{\(SP \((10\^6, 4, 10)\)\), \(10\^6\), "50", "4", \(\(>\)\(300\)\),
"76.5", "76.5", "420", "1031"}
}]], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData[{
"Random dense polynomials of degree ",
Cell[BoxData[
\(TraditionalForm\`n\)]],
" with ",
Cell[BoxData[
\(TraditionalForm\`k\)]],
" integer roots"
}], "Subsection"],
Cell[BoxData[
\(DP \((n,
k)\) := \((\[Sum]\+\(i = 0\)\%\(n - k\)\(a\_i\)
t\^i)\) \(\[Product]\+\(j = 1\)\%k\((t -
r\_j)\)\)\)], "DisplayFormula"],
Cell[TextData[{
"where ",
Cell[BoxData[
\(TraditionalForm\`a\_i\)]],
" and ",
Cell[BoxData[
\(TraditionalForm\`r\_j\)]],
" are randomly generated 100-bit integers, and ",
Cell[BoxData[
\(TraditionalForm\`r\_j\)]],
" are all distinct. "
}], "Text"],
Cell["\<\
When the polynomials are dense, HIR is clearly the fastest method. \
FIR is slower, but still much faster than RIR.\
\>", "Text"],
Cell[BoxData[GridBox[{
{"Polynomial", "Deg", "Roots", \(Time\[IndentingNewLine]
HIR\), \(Time\[IndentingNewLine]
FIR\), \(Time\[IndentingNewLine]
RIR\), \(st #\[IndentingNewLine]
FIR\), \(st #\[IndentingNewLine]
RIR\)},
{\(DP \((128, 128)\)\), "128", "128", "1.09",
"28.3", \(\(>\)\(300\)\), "22372", "?"},
{\(DP \((128, 32)\)\), "128", "32", "0.23", "1.33", "130", "6548",
"734241"},
{\(DP \((128, 8)\)\), "128", "8", "0.086", "0.26", "12.2", "2337",
"218165"},
{\(DP \((128, 2)\)\), "128", "2", "0.052", "0.13", "2.1", "879",
"58697"},
{\(DP \((128, 0)\)\), "128", "0", "0.008", "0.015", "0.11", "82",
"545"}
}]], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData[{
"Sparse polynomials with ",
Cell[BoxData[
\(TraditionalForm\`k + 1\)]],
" terms and the maximal number ",
Cell[BoxData[
\(TraditionalForm\`2 k + 1\)]],
"of integer roots"
}], "Subsection"],
Cell[BoxData[
\(PP \((k, p)\) :=
t \(\[Product]\+\(j = 1\)\%k\((t\^p - r\_j\^p)\)\)\)], "DisplayFormula"],
Cell[TextData[{
"where ",
Cell[BoxData[
\(TraditionalForm\`p\)]],
"is an even integer and ",
Cell[BoxData[
\(TraditionalForm\`r\_j\)]],
" are randomly generated 10-bit integers. "
}], "Text"],
Cell["\<\
For random sparse polynomials with the maximal number of roots the \
fastest algorithm is FIR. HIR gives the same timings as FIR for smaller \
degree polynomials, but for higher degrees it becomes much slower. RIR is \
slower than FIR by a factor which decreases with the degree.\
\>", "Text"],
Cell[BoxData[GridBox[{
{"Polynomial", "Deg", "Terms", "Roots", \(Time\[IndentingNewLine]
HIR\), \(Time\[IndentingNewLine]
FIR\), \(Time\[IndentingNewLine]
RIR\), \(st #\[IndentingNewLine]
FIR\), \(st #\[IndentingNewLine]
RIR\)},
{\(PP \((5, 2)\)\), \(10\^1\), "6", "10", "0.003", "0.003", "0.007",
"139", "337"},
{\(PP \((5, 20)\)\), \(10\^2\), "6", "10", "0.006", "0.007", "0.012",
"128", "314"},
{\(PP \((5, 200)\)\), \(10\^3\), "6", "10", "0.018", "0.018", "0.034",
"126", "268"},
{\(PP \((5, 2000)\)\), \(10\^4\), "6", "10", "16.3", "0.46", "0.51",
"129", "142"},
{\(PP \((5, 20000)\)\), \(10\^5\), "6", "10", \(\(>\)\(300\)\),
"14.5", "16.3", "126", "142"}
}]], "Text"]
}, Closed]]
}, Closed]]
}, Open ]]
},
FrontEndVersion->"4.1 for X",
ScreenRectangle->{{0, 1024}, {0, 768}},
Evaluator->"m50",
WindowSize->{802, 645},
WindowMargins->{{92, Automatic}, {Automatic, 17}},
Magnification->1.25,
StyleDefinitions -> "ArticleClassic.nb"
]
(*******************************************************************
Cached data follows. If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of the file. The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[1727, 52, 90, 3, 215, "Title"],
Cell[1820, 57, 34, 0, 48, "Author"],
Cell[1857, 59, 95, 3, 22, "Address"],
Cell[1955, 64, 1304, 29, 197, "Abstract"],
Cell[CellGroupData[{
Cell[3284, 97, 31, 0, 66, "Section"],
Cell[3318, 99, 19, 0, 31, "Text"],
Cell[3340, 101, 108, 2, 31, "NumberedEquation"],
Cell[3451, 105, 238, 8, 31, "Text"],
Cell[3692, 115, 66, 1, 31, "NumberedEquation"],
Cell[3761, 118, 223, 8, 31, "Text"],
Cell[3987, 128, 158, 2, 61, "NumberedEquation"],
Cell[4148, 132, 971, 25, 107, "Text"],
Cell[5122, 159, 418, 12, 50, "Text"],
Cell[5543, 173, 756, 25, 69, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[6336, 203, 41, 0, 33, "Section"],
Cell[6380, 205, 516, 14, 50, "Text"],
Cell[6899, 221, 113, 2, 31, "NumberedEquation"],
Cell[7015, 225, 19, 0, 31, "Text"],
Cell[7037, 227, 74, 1, 30, "NumberedEquation"],
Cell[7114, 230, 555, 17, 56, "Text"],
Cell[7672, 249, 206, 7, 31, "Text"],
Cell[7881, 258, 213, 4, 70, "NumberedEquation"],
Cell[8097, 264, 252, 9, 31, "Text"],
Cell[8352, 275, 190, 3, 49, "NumberedEquation"],
Cell[8545, 280, 261, 11, 31, "Text"],
Cell[8809, 293, 573, 19, 54, "Text"],
Cell[9385, 314, 182, 7, 32, "Text"],
Cell[9570, 323, 226, 8, 35, "Text"],
Cell[9799, 333, 1084, 39, 81, "Text"],
Cell[10886, 374, 499, 16, 50, "Text"],
Cell[11388, 392, 108, 2, 66, "DisplayFormula"],
Cell[11499, 396, 190, 8, 31, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[11726, 409, 32, 0, 33, "Section"],
Cell[11761, 411, 1019, 23, 145, "Text"],
Cell[12783, 436, 740, 20, 69, "Text"],
Cell[13526, 458, 902, 30, 69, "Text"],
Cell[14431, 490, 290, 10, 31, "Text"],
Cell[14724, 502, 932, 38, 69, "Text"],
Cell[15659, 542, 343, 14, 31, "Text"],
Cell[16005, 558, 266, 11, 31, "Text"],
Cell[16274, 571, 315, 11, 31, "Text"],
Cell[16592, 584, 321, 11, 31, "Text"],
Cell[16916, 597, 1018, 44, 50, "Text"],
Cell[17937, 643, 172, 5, 31, "Text"],
Cell[18112, 650, 697, 19, 88, "Text"],
Cell[18812, 671, 706, 24, 69, "Text"],
Cell[19521, 697, 234, 6, 50, "Text"],
Cell[19758, 705, 549, 17, 51, "Text"],
Cell[20310, 724, 513, 17, 50, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[20860, 746, 38, 0, 33, "Section"],
Cell[20901, 748, 1169, 33, 126, "Text"],
Cell[22073, 783, 600, 18, 69, "Text"],
Cell[22676, 803, 1430, 41, 126, "Text"],
Cell[24109, 846, 637, 15, 88, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[24783, 866, 44, 0, 33, "Section"],
Cell[24830, 868, 807, 25, 88, "Text"],
Cell[25640, 895, 282, 10, 31, "Text"],
Cell[25925, 907, 320, 14, 31, "Text"],
Cell[26248, 923, 248, 11, 31, "Text"],
Cell[26499, 936, 514, 16, 69, "Text"],
Cell[27016, 954, 509, 18, 52, "Text"],
Cell[27528, 974, 367, 12, 33, "Text"],
Cell[27898, 988, 759, 21, 88, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[28694, 1014, 39, 0, 33, "Section"],
Cell[28736, 1016, 757, 12, 126, "Text"],
Cell[29496, 1030, 99, 3, 31, "Text"],
Cell[CellGroupData[{
Cell[29620, 1037, 202, 8, 51, "Subsection"],
Cell[29825, 1047, 258, 5, 68, "DisplayFormula"],
Cell[30086, 1054, 396, 14, 51, "Text"],
Cell[30485, 1070, 454, 7, 88, "Text"],
Cell[30942, 1079, 1458, 27, 310, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[32437, 1111, 201, 8, 30, "Subsection"],
Cell[32641, 1121, 190, 4, 68, "DisplayFormula"],
Cell[32834, 1127, 279, 11, 32, "Text"],
Cell[33116, 1140, 139, 3, 31, "Text"],
Cell[33258, 1145, 786, 17, 164, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[34081, 1167, 227, 8, 30, "Subsection"],
Cell[34311, 1177, 116, 2, 68, "DisplayFormula"],
Cell[34430, 1181, 213, 8, 32, "Text"],
Cell[34646, 1191, 303, 5, 69, "Text"],
Cell[34952, 1198, 837, 17, 180, "Text"]
}, Closed]]
}, Closed]]
}, Open ]]
}
]
*)
(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)