(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the notebook
starts with the line of 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[ 16977, 527]*)
(*NotebookOutlinePosition[ 17905, 560]*)
(* CellTagsIndexPosition[ 17861, 556]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell[TextData[{
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[" Pearls",
Evaluatable->False,
AspectRatioFixed->True]
}], "Title",
Evaluatable->False,
AspectRatioFixed->True],
Cell["Problems and Solutions", "Subsubtitle",
Evaluatable->False,
AspectRatioFixed->True],
Cell["by Don Piele", "Subsubtitle",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox["Welcome back to ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" Pearls, the column devoted to examining interesting solutions to an \
assortment of simple yet appealing problems. The challenge for this issue is \
Keith numbers. Paul Wellin suggested this problem, which he discovered on \
the home page of Mike Keith at http://users.aol.com. There Mike writes about \
his numbers.",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell["Keith numbers", "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
" \"Fermat numbers, Mersenne numbers, Cullen numbers, Carmichael numbers, \
Smith numbers, Stirlingnumbers, Kaprekar numbers, Catalan numbers...and the \
list goes on. It seems like anybody who's anybody has a class of numbers with \
certain mathematical properties named after him or her. Well, it took a \
while, but finally, in 1987, I published the first paper on what I hereby \
declare are now called Keith numbers. Although completely and utterly \
recreational (perhaps that explains some of their charm), they have proved \
to be popular enough to inspire several papers by other authors, most \
notably a whole series of papers that appeared in 1994 in Volume 26, Number \
3 of the ",
StyleBox["Journal of Recreational Mathematics",
FontVariations->{"Underline"->True}],
". \n\n\"A Keith number is an ",
StyleBox["n",
FontSlant->"Italic"],
"-digit integer ",
StyleBox["N",
FontSlant->"Italic"],
" with the following property: If a Fibonacci-like sequence (in which each \
term in the sequence is the sum of the ",
StyleBox["n",
FontSlant->"Italic"],
" previous terms) is formed, with the first ",
StyleBox["n",
FontSlant->"Italic"],
" terms being the decimal digits of the number ",
StyleBox["N",
FontSlant->"Italic"],
", then ",
StyleBox["N",
FontSlant->"Italic"],
" itself occurs as a term in the sequence. For example, 197 is a Keith \
number since it generates the sequence \n\n 1, 9, 7, 17, 33, \
57, 107, 197, ... "
}], "Text"],
Cell[TextData[
"\n\n\n\"Keith numbers are intriguing for several reasons. They are very \
hard to find, the only methods known being essentially exhaustive search. \
Some techniques are available to speed up the search, but there is no known \
technique for finding a Keith number 'quickly'. They are in some ways \
reminiscent of the primes in their erratic distribution among the integers. \
However, they are much rarer than the primes \[Dash] there are only 52 Keith \
numbers less than 15 digits long. \n\n \" For your amusement, here is the \
complete list of Keith numbers less than 15 digits in length: \n\n \
14 19 28 47 61 75 \n 197 742\n 1104 1537 2208 2580 3684 4788 \
7385 7647 7909 \n 31331 34285 34348 55604 62662 86936 93993 \n \
120184 129106 147640 156146 174680 183186 298320 355419 694280 926993 \n \
1084051 7913837 \n 11436171 33445755 441212607 \n \
129572008 251133297 \n (none with 10 digits) \n 24769286411 \
96189170155 \n 171570159070 202366307758 239143607789 296658839738 \n\
1934197506555 8756963649152 \n 43520999798747 \
74596893730427 97295849958669 "], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
\"In addition, three 15-digit Keith numbers are known (but not \
necessarily all of them).
\"Another feature of these numbers is that there are so many \
unanswered questions. Here are two of them:
1.Is the number of Keith numbers finite or infinite?
2.Define a cluster of Keith numbers as a set of two or
more Keith numbers (all
having the same number of digits) in which all the \
numbers are integer multiples of
the smallest number in the set. There are only three \
known clusters: (14, 28),
(1104, 2208), and the remarkable (for having three \
members) (31331, 62662,
93993). Question: is the number of Keith clusters
finite or infinite?
\"We conjecture that the answer to #1 is \"infinite\" and the
answer to #2 is \"finite\", but a proof of neither result is known. Since we
suspect that there are an infinite number of Keith numbers, the problem of
finding the next such number always remains a tantalizing one. \"\
\>", "Text"],
Cell[CellGroupData[{
Cell["Keith[n]", "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox["Our problem for this issue will be to write the function ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Keith[n]",
Evaluatable->False,
AspectRatioFixed->True,
FontWeight->"Bold"],
StyleBox[
" which returns the list of all n digit Keith numbers. Clearly, speed is \
important here since the order of your algorithm is at least ",
Evaluatable->False,
AspectRatioFixed->True],
Cell[BoxData[
\(TraditionalForm\`\(O(10\^n\)\)]],
StyleBox[
"). If you can do better than that, you will have discovered a truely \
unique pearl.\n",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Open ]]
}, Open ]],
Cell[" factorialPrimeDecomposition[n]", "Subsection"],
Cell["\<\
In the last issue, the following problem was posed. Create a \
function, factorialPrimeDecomposition[n], that will generate the prime \
decomposition for n! Speed counts, so time your solution for n=5000 and \
compare it with
the time for FactorInteger[n!]. A good solution will have more than a \
10-fold speed increase. Like many good
pearls, it can also be written as a one-liner. \
\>", "Text"],
Cell[TextData[{
"Solution\n",
StyleBox["by Michael Trott (mtrott@wolfram.com) ",
FontWeight->"Plain",
FontSlant->"Italic"]
}], "Subsection"],
Cell["\<\
The solution that appears at the end of this discussion was \
submitted by Michael Trott . It is the same algorithm that I used but
with some nice optimizations added. I have added an introductory example to
make it a bit easier for less experienced readers.
Let's walk through the algorithm by taking a small prime, say 7, and see
how to find the power to which this prime occurs in the prime factorization
of 100! This is called the multiplicity for 7. \
\>", "Text"],
Cell["\<\
The number of times that 7 divides 100 tells us something very \
important. \
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(Floor[100/7]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`14\)], "Output"]
}, Open ]],
Cell["\<\
It tells us there are 14 numbers between 1 and 100 that have 7 as a \
factor. The numbers are:\
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(7\ Range[14]\)], "Input"],
Cell[BoxData[
\(TraditionalForm
\`{7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98}\)], "Output"]
}, Open ]],
Cell[TextData[{
"The number of times that ",
Cell[BoxData[
\(TraditionalForm\`7\^\(2\ \)\ divides\)]],
"100 tells us how many numbers between 1 and 100 have ",
Cell[BoxData[
\(TraditionalForm\`7\^2\)]],
" as a factor."
}], "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(Floor[100/7\^2]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`2\)], "Output"]
}, Open ]],
Cell["And the numbers are:", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(\(7\^\(2\ \)\) Range[2]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`{49, 98}\)], "Output"]
}, Open ]],
Cell[TextData[{
"There are clearly none with ",
Cell[BoxData[
\(TraditionalForm\`\(7\^3\ \)\)]],
"as a factor, since ",
Cell[BoxData[
\(TraditionalForm\`\(7\^3\ \)\)]],
"is larger than 100. Thus, we can stop checking as soon as\n",
Cell[BoxData[
\(TraditionalForm\`\(7\^x > 100\ \)\)]],
"or x > Log[7,100]. This means we only need to check up to \
Floor[Log[7,100]]. Remember Log[b,x] is the Log of x to the base b."
}], "Text"],
Cell["\<\
So the total number of powers of 7 in the prime factorization of \
100! is 14+2 or 16. Let's check it out. \
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(\(FactorInteger[\(100!\)]\)[\([4]\)]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`{7, 16}\)], "Output"]
}, Open ]],
Cell[TextData[{
"Clearly, for 100! the only primes we need to consider are those less than \
or equal to 100. These are easy to list using\nPrimePi[k] (the number of \
primes ",
Cell[BoxData[
\(TraditionalForm\`\(\ \[LessEqual] \ \)\)]],
"k) and Prime[k] (the kth prime)."
}], "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(Prime[Range[PrimePi[100]]]\)], "Input"],
Cell[BoxData[
\(TraditionalForm
\`{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97}\)], "Output"]
}, Open ]],
Cell[TextData[{
"So the algorithm for factorialPrimeDecomposition[n] is as follows:\n \
",
StyleBox[
" For every prime p \[LessEqual] n, \n the \
multiplicity for prime p = Sum[ ",
FontWeight->"Bold"],
Cell[BoxData[
\(TraditionalForm\`Floor[n/p\^\(i\ \)]\)],
FontWeight->"Bold"],
StyleBox[", {i, 1, Log[p,n]}",
FontWeight->"Bold"],
"\n "
}], "Text"],
Cell["Now let's implement this algorithm in Mathematica. ", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(n = 100; \n
Table[{Prime[k],
Sum[Floor[n/Prime[k]^i], {i, Floor[Log[Prime[k], n]]}]}, {k,
PrimePi[n]}]\)], "Input"],
Cell[BoxData[
\({{2, 97}, {3, 48}, {5, 24}, {7, 16}, {11, 9}, {13, 7}, {17, 5}, {19,
5}, {23, 4}, {29, 3}, {31, 3}, {37, 2}, {41, 2}, {43, 2}, {47, 2}, {
53, 1}, {59, 1}, {61, 1}, {67, 1}, {71, 1}, {73, 1}, {79, 1}, {83,
1}, {89, 1}, {97, 1}}\)], "Output"]
}, Open ]],
Cell[TextData[
"Of course, this could also be done without the Table command by mapping a \
pure function onto the primes \[LessEqual] n."], "Text"],
Cell[BoxData[
\(\({#, \ Sum[Floor[n/#^i], \n\ {i, \ Floor[Log[#, \ n]]}]}&\)/@\ \
Prime[Range[PrimePi[n]]]\)], "Input"],
Cell["\<\
We can optimize this solution a bit by noting that once a prime \
divides into n less than twice, then its multiplicity is automatically
one. Notice, that for the prime 53, Floor[100/53]=1, and the multiplity for \
53
and all higher primes is 1. So we add this optimization feature to the
algorithm with an If statement.\
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(\({#,
If[\ n/# < 2, 1, \
Sum[Floor[n/#^i], \n\ {i, \ Floor[Log[#, \ n]]}]]}&\)/@\ \
Prime[Range[PrimePi[n]]] // Timing\)], "Input"],
Cell[BoxData[
\({0.116666666666787932`\
Second, {{2, 97}, {3, 48}, {5, 24}, {7, 16}, {11, 9}, {13, 7}, {17,
5}, {19, 5}, {23, 4}, {29, 3}, {31, 3}, {37, 2}, {41, 2}, {43, 2}, {
47, 2}, {53, 1}, {59, 1}, {61, 1}, {67, 1}, {71, 1}, {73, 1}, {79,
1}, {83, 1}, {89, 1}, {97, 1}}}\)], "Output"]
}, Open ]],
Cell[TextData[{
"Finally, compiling always helps speed up simple arithmetic \
computations.The first part of the compile function specifies the number \
type for n. The additional code at the end is there to specify a type for the \
function ",
StyleBox["Prime[Range[PrimePi[n]]]",
FontWeight->"Bold"],
" used in the code. It tells the system it is a one dimensional list of \
integers."
}], "Text"],
Cell[BoxData[
\(\( (*\ Solution\ submitted\ by\ Michael\ Trott*) \n
\(factorialPrimeDecomposition\ =
Compile[{{n, _Integer}}, \n
\({#, \ If[\ n/# < 2\ , \ 1, \
Sum[Floor[n/#^i], \n\ {i, \ Floor[Log[#, \ n]]}]]}&\)\ /@\ \
Prime[Range[PrimePi[n]]], \n
\ \ \ \ \ \ \ \ \ {{Prime[Range[PrimePi[n]]], \ _Integer, 1}}];
\)\)\)], "Input"],
Cell[TextData[{
"Now we are ready to check the time improvement over the built-in ",
StyleBox["FactorInteger",
FontWeight->"Bold"],
" function. The times below are for Power Mac 7100/66 and will vary with \
other systems."
}], "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(\(\ time1\ = \
\(Timing[\((facs1 = FactorInteger[\(5000!\)])\); ]\)[\([1]\)]\)\)],
"Input"],
Cell[BoxData[
\(TraditionalForm\`54.1666666666666696`\ Second\)], "Output"]
}, Open ]],
Cell["Let's run the fast algorithm 100 times.", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(time2 =
\(Timing[\n\ \ \ \ \ \ \ \
Do[\((facs2 = factorialPrimeDecomposition[5000])\), \n{100}]]\)[
\([1]\)]\)], "Input"],
Cell[BoxData[
\(17.8666666666649698`\ Second\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(speedup\ = 100\ \ time1/\ time2\)], "Input"],
Cell[BoxData[
\(303.171641791073609`\)], "Output"]
}, Open ]],
Cell["The speedup factor is somewhere between 300 and 400.", "Text"],
Cell["Just to make sure, let's check that the solutions agree.", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(facs1 == facs2\)], "Input"],
Cell[BoxData[
\(True\)], "Output"]
}, Open ]],
Cell["Show me the pearls!", "Subsection"],
Cell["\<\
Have you got a Mathematica Pearl burried deep in your files? I \
invite you to share them in this column and add them to our string. If
you have the goods, then \"Show me the pearls!\" The best way to get my
attention is to send your work via email to . Little pearls \
slip nicely into an email message and can zip to me from anywhere. You will
get immediate feedback. Accompany your work with an explanation of how it was \
created. That's it, and thanks.\
\>", "Text"],
Cell["About The Editor", "Subsubsection"],
Cell[TextData[{
"Don Piele has been interested in creating programming problems since he \
began the International Computer Problem Solving Contest in the pages of ",
StyleBox["Creative Computing",
FontVariations->{"Underline"->True}],
" in 1981. In 1992, he organized the USA Computing Olympiad, which \
selects the top four American high school computer programmers to represent \
the USA at the annual International Computing Olympiad. He also writes the \
column, ",
StyleBox["Cowculations",
FontSlant->"Italic"],
", devoted to computer algorithms using Mathematica in ",
StyleBox["Quantum Magazine",
FontVariations->{"Underline"->True}],
". The Web address for the USACO is: \[Not]http:// usaco.uwp.edu."
}], "Text"],
Cell["\<\
Don Piele
Mathematics Department
University of Wisconsin-Parkside
Kenosha, WI 53141
piele@cs.uwp.edu
\
\>", "Subsubsection"]
},
FrontEndVersion->"NeXT 3.0",
ScreenRectangle->{{0, 1053}, {0, 832}},
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{638, 440},
WindowMargins->{{157, Automatic}, {Automatic, 169}},
PrintingCopies->1,
PrintingPageRange->{1, Automatic},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1.25
]
(***********************************************************************
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[1731, 51, 254, 10, 122, "Title"],
Cell[1988, 63, 93, 2, 56, "Subsubtitle"],
Cell[2084, 67, 83, 2, 56, "Subsubtitle"],
Cell[2170, 71, 654, 18, 101, "Text"],
Cell[2827, 91, 83, 2, 54, "Subsection"],
Cell[2913, 95, 1540, 35, 332, "Text"],
Cell[4456, 132, 1259, 18, 542, "Text"],
Cell[5718, 152, 1090, 23, 416, "Text"],
Cell[CellGroupData[{
Cell[6833, 179, 78, 2, 54, "Subsection"],
Cell[6914, 183, 721, 22, 101, "Text"]
}, Open ]]
}, Open ]],
Cell[7662, 209, 53, 0, 54, "Subsection"],
Cell[7718, 211, 410, 8, 143, "Text"],
Cell[8131, 221, 152, 5, 73, "Subsection"],
Cell[8286, 228, 484, 9, 185, "Text"],
Cell[8773, 239, 101, 3, 38, "Text"],
Cell[CellGroupData[{
Cell[8899, 246, 45, 1, 31, "Input"],
Cell[8947, 249, 53, 1, 33, "Output"]
}, Open ]],
Cell[9015, 253, 118, 3, 38, "Text"],
Cell[CellGroupData[{
Cell[9158, 260, 45, 1, 31, "Input"],
Cell[9206, 263, 111, 2, 33, "Output"]
}, Open ]],
Cell[9332, 268, 253, 8, 59, "Text"],
Cell[CellGroupData[{
Cell[9610, 280, 48, 1, 34, "Input"],
Cell[9661, 283, 52, 1, 33, "Output"]
}, Open ]],
Cell[9728, 287, 36, 0, 38, "Text"],
Cell[CellGroupData[{
Cell[9789, 291, 56, 1, 34, "Input"],
Cell[9848, 294, 59, 1, 33, "Output"]
}, Open ]],
Cell[9922, 298, 461, 12, 101, "Text"],
Cell[10386, 312, 131, 3, 59, "Text"],
Cell[CellGroupData[{
Cell[10542, 319, 69, 1, 31, "Input"],
Cell[10614, 322, 58, 1, 33, "Output"]
}, Open ]],
Cell[10687, 326, 298, 7, 80, "Text"],
Cell[CellGroupData[{
Cell[11010, 337, 59, 1, 31, "Input"],
Cell[11072, 340, 159, 3, 52, "Output"]
}, Open ]],
Cell[11246, 346, 431, 13, 101, "Text"],
Cell[11680, 361, 67, 0, 38, "Text"],
Cell[CellGroupData[{
Cell[11772, 365, 157, 4, 82, "Input"],
Cell[11932, 371, 286, 4, 83, "Output"]
}, Open ]],
Cell[12233, 378, 149, 2, 59, "Text"],
Cell[12385, 382, 130, 2, 48, "Input"],
Cell[12518, 386, 348, 7, 122, "Text"],
Cell[CellGroupData[{
Cell[12891, 397, 188, 4, 65, "Input"],
Cell[13082, 403, 334, 5, 100, "Output"]
}, Open ]],
Cell[13431, 411, 410, 9, 101, "Text"],
Cell[13844, 422, 396, 8, 99, "Input"],
Cell[14243, 432, 245, 6, 59, "Text"],
Cell[CellGroupData[{
Cell[14513, 442, 122, 3, 31, "Input"],
Cell[14638, 447, 79, 1, 33, "Output"]
}, Open ]],
Cell[14732, 451, 55, 0, 38, "Text"],
Cell[CellGroupData[{
Cell[14812, 455, 167, 4, 65, "Input"],
Cell[14982, 461, 62, 1, 32, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[15081, 467, 65, 1, 31, "Input"],
Cell[15149, 470, 54, 1, 32, "Output"]
}, Open ]],
Cell[15218, 474, 68, 0, 38, "Text"],
Cell[15289, 476, 72, 0, 38, "Text"],
Cell[CellGroupData[{
Cell[15386, 480, 47, 1, 31, "Input"],
Cell[15436, 483, 38, 1, 32, "Output"]
}, Open ]],
Cell[15489, 487, 41, 0, 54, "Subsection"],
Cell[15533, 489, 504, 8, 164, "Text"],
Cell[16040, 499, 41, 0, 51, "Subsubsection"],
Cell[16084, 501, 752, 15, 164, "Text"],
Cell[16839, 518, 134, 7, 146, "Subsubsection"]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)