Mathematica Pearls

Problems and Solutions

by Don Piele

Welcome back to Mathematica 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.

Keith numbers

"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 Journal of Recreational Mathematics. 

"A Keith number is an n-digit integer N with the  following property: If a Fibonacci-like sequence (in which each term in the  sequence is the sum of the n previous terms) is formed, with the first n terms being  the decimal digits of the number N, then N itself occurs as a term in the sequence. For example, 197 is a Keith number since it generates the  sequence 

              1, 9, 7, 17, 33, 57, 107, 197, ...

"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 – there are only 52  Keith numbers less than 15 digits long. 

 " For your amusement, here is the complete list of Keith numbers less than 15 digits in length:  

         14 19 28 47 61 75 
         197 742
         1104 1537 2208 2580 3684 4788 7385 7647 7909 
         31331 34285 34348 55604 62662 86936 93993  
         120184 129106 147640 156146 174680 183186 298320 355419 694280 926993  
         1084051 7913837 
         11436171 33445755 441212607 
         129572008 251133297 
         (none with 10 digits) 
          24769286411 96189170155 
         171570159070 202366307758 239143607789  296658839738 
         1934197506555 8756963649152 
         43520999798747  74596893730427 97295849958669

"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. "

Keith[n]

Our problem for this issue will be to write the function Keith[n] which returns the list of all n digit Keith numbers. Clearly, speed  is important here since the order of your algorithm is at least O((10^n)). If you can do better than that, you will have discovered a truely unique pearl.

factorialPrimeDecomposition[n]

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.

Solution
by Michael Trott (mtrott@wolfram.com)

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.

The number of times that 7  divides 100 tells us something very important.

Floor[100/7]

It tells us there are 14 numbers between 1 and 100 that have 7 as a factor. The numbers are:

7 Range[14]

The number of times that (7^(2 )) divides100 tells us how many numbers between 1 and 100 have  (7^2) as a factor.

Floor[100/(7^2)]

And the numbers are:

(7^(2 )) Range[2]

There are clearly none with (7^3) as a factor, since (7^3) is larger than 100. Thus, we can stop checking as soon as
(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.

So the total number of powers of 7 in the prime factorization of 100! is 14+2 or 16. Let's check it out.

FactorInteger[100!][[4]]

Clearly, for 100! the only primes we need to consider are those less  than or equal to 100. These are easy to list using
PrimePi[k] (the number of primes   ≤  k) and Prime[k] (the kth prime).

Prime[Range[PrimePi[100]]]

So the algorithm for factorialPrimeDecomposition[n] is as  follows:
         For every prime p ≤ n, 
                       the multiplicity for prime p = Sum[ Floor[n/(p^(i ))], {i, 1, Log[p,n]}

Now let's implement this algorithm in Mathematica.

\(n = 100; 

    Table[{Prime[k], 
        Sum[Floor[n/Prime[k]^i], {i, Floor[Log[Prime[k], n]]}]}, {k, 
        PrimePi[n]}]\)

Of course, this could also be done without the Table command by mapping  a pure function onto the primes ≤ n.

\(\({#, \ Sum[Floor[n/#^i], 
\ {i, \ Floor[Log[#, \ n]]}]}&\)/@\ \ 
      Prime[Range[PrimePi[n]]]\)

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.

\(\({#, 
            If[\ n/# < 2, 1, \ 
              Sum[Floor[n/#^i], 
\ {i, \ Floor[Log[#, \ n]]}]]}&\)/@\ \ 
        Prime[Range[PrimePi[n]]] // Timing\)

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 Prime[Range[PrimePi[n]]] used in the code. It tells the system it is a one dimensional list of integers.

\(\( (*\ Solution\ submitted\ by\ Michael\ Trott*) 

    \(factorialPrimeDecomposition\  = 
      Compile[{{n, _Integer}}, 

        \({#, \ If[\ n/# < 2\ , \ 1, \ 
                Sum[Floor[n/#^i], 
\ {i, \ Floor[Log[#, \ n]]}]]}&\)\ /@\ \ 
          Prime[Range[PrimePi[n]]], 

        \ \ \ \ \ \ \ \ \ {{Prime[Range[PrimePi[n]]], \ _Integer, 1}}]; 
    \)\)\)

Now we are ready to check the time improvement over the built-in FactorInteger function. The times below are for Power Mac 7100/66 and will vary  with other systems.

time1  =  
      Timing[(facs1 = FactorInteger[5000!]); ][[1]]

Let's run the fast algorithm 100 times.

\(time2 = 
      \(Timing[
\ \ \ \ \ \ \ \ 
          Do[\((facs2 = factorialPrimeDecomposition[5000])\), 
{100}]]\)[
        \([1]\)]\)

speedup  = 100  time1/ time2

The speedup factor is somewhere between 300 and 400.

Just to make sure, let's check that the solutions agree.

facs1 == facs2

Show me the pearls!

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 <piele@cs.uwp.edu>. 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.

About The Editor

Don Piele has been interested in creating programming problems since  he began the International Computer Problem Solving Contest in the pages of Creative Computing 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, Cowculations, devoted to computer algorithms using Mathematica in Quantum Magazine. The Web address for the USACO is: ¬http:// usaco.uwp.edu.

Don Piele
Mathematics Department
University of Wisconsin-Parkside
Kenosha, WI 53141
piele@cs.uwp.edu