Programmatic Mathematica XVII: The Collatz Conjecture

There has been a lot of interest recently in the Collatz Conjecture. A lot of video blogs are going into it, particularly Numberphile, a vlog present on YouTube. It might have something to do with the fact that this year is the 70th anniversary of the conjecture. It is a simple idea, easy enough for a child to understand. Yet, it has been difficult enough that no one has been able to either prove or disprove it to this day.

The Collatz Conjecture is the hunch, or guess, or idea, that performing a certain recursive operation on any positive integer leads to the inevitable result that repeated operations on all successors will lead to the number 1. After that, the sequence of {1, 4, 2, …} occurs in an infinite repetition.

This problem was first posed by Lothar Collatz in 1937. The reason it is only a conjecture is that no one has been able to prove it for all positive integers. It is only conjectured to work as such. Over the past seventy years, no one has been able to furnish a counterexample where the number 1 is not reached. So by now, we’re “pretty sure” Collatz is correct for all positive integers.

I thought of some Mathematica code to write for this. The algorithm would go something like:

  1. Precondition: n > 0; n \in Z
  2. If n is 1, return 1 and exit
  3. If n is even, return n/2
  4. If n is odd, return 3n + 1
  5. Go back to line 2.

Like Fermat’s Last Theorem, which has been proved once and for all in 1995 by Professor Andrew Wiles, and aided by Richard Taylor, the Collatz Conjecture is simple enough to describe to any lay person (as I just did), but its proof has eluded us.

The application of the above algorithm to Mathematica code involves some new syntax. Sow[n] acts as a kind of array for anyone who doesn’t want to declare and implement an array. I would suppose that the programmers of the Mathematica language didn’t see the need for an array for many implementations, such as sequences of numbers. If you want to generate a sequence, you want the numbers in order from some lower bound, up to some upper bound. If you want to list them, you want to do the same thing. It is not often that you want to access only one particular value inside the sequence. This is for those people who just want the whole sequence uninterrupted.

I guess what Sow[n] does is leave the members of the sequence lying around in some pre-defined region in computer memory. That memory is likely to be freed once the Reap[n] function is called, which lists all the members of the stored sequence in the order generated.

EvenQ[] and OddQ[] are employed to check if n if odd or even before executing the rest of the line. If false, control passes through the next line. The testing is inefficient here, since each statement is tested all the time. So, if we already know the number is even, OddQ[] is executed anyway.

ClearAll[Co];
Co[1] = 1;
Co[n_ /; EvenQ[n]] := (Sow[n]; Co[n/2])
Co[n_ /; OddQ[n]] := (Sow[n]; Co[3*n + 1])
Collatz[n_] := Reap[Co[n]]

But Reap[n] by itself gives a nested array (or more accurately, a “ragged” array) with the final “1” outside of the innermost nesting, where the other numbers are.

In[10]:= Collatz[7]
Out[10]= {1, {{7, 22, 11, 34, 17, 52, 26, 13,
    40, 20, 10, 5, 16, 8, 4, 2}}}

Nested arrays are un-necessary, but the remedy to this gets rid of the number “1” which is the number the Collatz function is supposed to always land on. So we then rely on the presence of the number “2”, the number arrived at before going to “1”, at the end of the sequence. Getting rid of the nested array relies on using Flatten[Reap[Co[n]]]. But when you do that, this happens:

In[11]:= Collatz[7]
Out[11]= {1, 7, 22, 11, 34, 17, 52, 26, 13,
    40, 20, 10, 5, 16, 8, 4, 2}

Flattening has the effect of placing the ending 1 at the beginning of the array. If we can live with this minor inconvenience, then we are able to test the Collatz Conjecture on wide ranges of positive integers. So, this is the code we ended up with:

ClearAll[Co];
Co[1] = 1;
Co[n_ /; EvenQ[n]] := (Sow[n]; Co[n/2])
Co[n_ /; OddQ[n]] := (Sow[n]; Co[3*n + 1])
Collatz[n_] := Flatten[Reap[Co[n]]]

The sequences generated by the Collatz conjecture have the well-documented property of having common endings. Using the Table[] command, we can observe the uncanny phenomena that most of these sequences end in “8, 4, 2” (or, to be more precise, “8, 4, 2, 1”). Here are the sequences generated for the numbers from 1 to 10:

In[38]:= Table[Collatz[i], {i, 10}]

Out[38]= {{1}, 
          {1, 2}, 
          {1, 3, 10, 5, 16, 8, 4, 2}, 
          {1, 4, 2}, 
          {1, 5, 16, 8, 4, 2}, 
          {1, 6, 3, 10, 5, 16, 8, 4, 2}, 
          {1, 7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2},
          {1, 8, 4, 2}, 
{1, 9, 28, 14, 7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2},    
          {1, 10, 5, 16, 8, 4, 2}}

Because even numbers are to be divided by 2, somewhere along the meanderings of the sequence, a power of 2 is encountered, and from there it’s a one-way trip to the number “1”.

Programmatic Mathematica XVI: Patterns in Highly Composite Numbers

This article was inspired by a vlog from Numberphile, on the discussion of “5040: an anti-prime number”, or some title like that.

A contributor to the OEIS named Jean-François Alcover came up with a short bit of Mathematica code that I modified slightly:

Reap[
   For[
      record = 0; n = 1, n <= 110880, n = If[n < 60, n + 1, n + 60], tau = DivisorSigma[0, n]; 
      If[tau > record, record = tau; Print[n, "\t\t", tau];
      Sow[tau]]]][[2, 1]]

This generates a list of a set of numbers with an unusually high amount of factors called “highly composite numbers” up to 110,880. The second column of the output are the number of factors.

1      1
2       2
4       3
6       4
12      6
24      8
36      9
48      10
60      12
120     16
180     18
240     20
360     24
720     30
840     32
1260        36
1680        40
2520        48
5040        60
7560        64
10080       72
15120       80
20160       84
25200       90
27720       96
45360       100
50400       108
55440       120
83160       128
110880      144

For a number like 110,880, there is no number before it that has more than 144 factors.

Highly composite numbers (HCNs) are loosely defined as a natural number which has more factors than any others that came before it. 12 is such a number, with 6 factors, as is 6 itself with 4. The number 5040 has 60 factors, and is also considered highly composite.

5040=24×32×5×7
This works out to 60, because with 24, for example, we get the factors 2, 4, 8, and 16. With 24×32, we get 2, 3, 4, 6, 8, 9, 16, 18, 36, 72, and 144, all which evenly divide 5040. The total number of factors including 1 and 5040 itself can be had from adding 1 to each exponent and multiplying: (4+1)(2+1)(1+1)(1+1)=5×3×2×2=60.

Initially, facotorization of HCNs was done in Maple using the “ifactor()” command. But there is a publication circulating the Internet referring to a table created by Ramanujan that has these factors. A partial list of these are summarized in a table below. The top row headers are the prime numbers that can be the prime factors, from 2 to 17. The first column is the number to factorize. The numbers in the same columns below these prime numbers are the exponents on the primes, such as: 10,080=25×32×51×71. The last column are the total number of factors on these HCNs. So, by adding 1 to each exponent in the row and multiplying, we find that 10,080 has 6×3×2×2=72 factors.

NUMBER PATTERNS OBSERVED

As a number of factors (underneath the “# facotrs” column), We get overlapping patterns starting from 60. One of them would be the sequence: 120, 240, 360, 480, 600, and 720. But the lack of an 840 breaks that pattern. But then we get 960, then 1080 is skipped, but then we get 1200.

For numbers of factors that are powers of 2, it seems to go right off the end of the table and beyond: 64, 128, 256, 512, 1024, 2048, 4096, 8192, … . Before 5040, the pattern is completed, since 2 has 2 factors, 6 has 4 factors, 24 has 8 factors, 120 has 16 factors, and 840 has 32 factors. The HCN with 8192 factors is 3,212,537,328,000. We have to go beyond that to see if there is a number with 16,384 factors.

Multiples of 12 make their appearance as numbers of factors: 12, 24, 36, 48, 60 (which are the numbers of factors of 5040), 72, 84, 96, 108, 120, but a lack of a 132 breaks that pattern. But then we see: 144, 288, 432, 576, 720, 864, 1008, 1152, and the pattern ends with the lack of a 1296.

We also observe short runs of numbers of factors in the sequence 100, 200, 400, 800, until we reach the end of this table. But the pattern continues with the number 2,095,133,040, which has 1600 factors. Then, 3200 is skipped.

There are also multiples of 200: 200, 400, 600, 800, but the lack of a 1000 breaks that pattern. But when seen as multiples of 400, we get: 400, 800, 1200, 1600, but then 2000 is skipped.

There are also peculiarities in the HCNs themselves. Going from 5040 to as high as 41,902,660,800, only 4 of the 60 HCNs were not multiples of 5040. The rest had the remainder 2520, which is one-half of 5040.

Also beginning from the HCN 720,720, we observe a run of numbers containing 3-digit repeats: 1081080, 1441440, 2162160, 2882880, 3603600, 4324320, 6486480, 7207200, 8648640, 10810800, and 14414400.

Number 2   3   5   7   11  13  17  # of
                                factors
-----------------------------------------------------------------------
5040    4   2   1   1               60  
7560    3   3   1   1               64  
10080   5   2   1   1               72  
15120   4   3   1   1               80  
20160   6   2   1   1               84  
25200   4   2   2   1               90  
27720   3   2   1   1   1           96  
45360   4   4   1   1               100 
50400   5   2   2   1               108 
55440   4   2   1   1   1           120 
83160   3   3   1   1   1           128 
110880  5   2   1   1   1           144 
166320  4   3   1   1   1           160 
221760  6   2   1   1   1           168 
332640  5   3   1   1   1           192 
498960  4   4   1   1   1           200 
554400  5   2   2   1   1           216 
665280  6   3   1   1   1           224 
720720  4   2   1   1   1   1       240 
1081080 3   3   1   1   1   1       256 
1441440 5   2   1   1   1   1       288 
2162160 4   3   1   1   1   1       320 
2882880 6   2   1   1   1   1       336 
3603600 4   2   2   1   1   1       360 
4324320 5   3   1   1   1   1       384 
6486480 4   4   1   1   1   1       400 
7207200 5   2   2   1   1   1       432 
8648640 6   3   1   1   1   1       448 
10810800    4   3   2   1   1   1       480 
14414400    6   2   2   1   1   1       504 
17297280    7   3   1   1   1   1       512 
21621600    5   3   2   1   1   1       576 
32432400    4   4   2   1   1   1       600 
61261200    4   2   2   1   1   1   1   720 
73513440    5   3   1   1   1   1   1   768 
110270160   4   4   1   1   1   1   1   800 
122522400   5   2   2   1   1   1   1   864 
147026880   6   3   1   1   1   1   1   896 
183783600   4   3   2   1   1   1   1   960 
245044800   6   2   2   1   1   1   1   1008    
294053760   7   3   1   1   1   1   1   1024    
367567200   5   3   2   1   1   1   1   1152    
551350800   4   4   2   1   1   1   1   1200    

After that run, we see a 4-digit overlapping repeat. The digits of the HCN 17297280 could be thought of as an overlap of 1728 and 1728 to make 1729728 as part of that number. The 3-digit run continues with: 21621600, 32432400, 61261200, and after that the pattern is broken.

Programmatic Mathematica XV: Lucas Numbers

The Lucas sequence follows the same rules for its generation as the Fibonacci sequence, except that the Lucas sequence begins with t1 = 2 and t2 =1.

Lucas numbers are found in the petal counts of flowing plants and pinecone spirals much the same as the Fibonnaci numbers. Also, like the Fibonacci numbers, successive pairs of Lucas numbers can be divided to make the Golden Ratio, \phi. The Mathematica version (10) which I am using has a way of  highlighting certain numbers that meet certain conditions. One of them is the Framed[] function, which draws a box around numbers. Framed[] can be placed into If[] statements so that an array of numbers can be fed into it (using a Table[] command).

For example, let’s frame all Lucas numbers that are prime:

In[1]:= If[PrimeQ[#], Framed[#], #] & /@ Table[L[n], {n, 0, 30}]

The If[] statement is best described as:

If[Condition[#], do_if_true[#], do_if_false[#]]

The crosshatch # is a positional parameter upon which some condition is placed by some function we are calling Condition[]. This boolean function returns True or False. In the statement we are using above, the function PrimeQ will return true if the number in the positional parameter is prime; false if 1 or composite.

The positional parameters require a source of numbers by which to make computations, so for this source, we shall look to a sequence of Lucas numbers generated by the Table command. The function which generates the numbers is a user-defined function L[n_]:

In[2]:= L[0] := 2
In[3]:= L[1] := 1
In[4]:= L[n_] := L[n-2] + L[n-1]

With that, I can generate an array with the Table[] command to get the first 31 Lucas numbers:

In[5]:= Table[L[n], {n, 0, 30}]
{2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843, 1364, \
2207, 3571, 5778, 9349, 15127, 24476, 39603, 64079, 103682, 167761, \
271443, 439204, 710647, 1149851, 1860498}

This list (or “table”) of numbers is passed through the If[] statement thusly:

In[6]:= If[PrimeQ[#], Framed[#], #] & /@ Table[L[n], {n, 0, 30}]

to produce the following output:

In[7]:=Array_Frame_Lucas_Prime

Note that this one was an actual screenshot, to get the effect of the boxes. So, these are the first 31 Lucas numbers, with boxes around the prime numbers. The Table[] command appears to feed the Lucas numbers into the positional parameters represented by #.

There was a sequence I created. Maybe it’s already famous; I have no idea. On the other hand, maybe no one cares. But I wanted to show that with any made-up sequence that is recursive in the same way Fibonacci and Lucas numbers were, that I could show, for example, that as the numbers grow, neighbouring numbers can get closer to the Golden Ratio. The Golden Ratio is \phi = \frac{1 + \sqrt{5}}{2}. I want to show that this is not really anything special that would be attributed to Fibonacci or François Lucas. It can be shown that, for any recursive sequence involving the next term being the sum of the previous two terms, sooner or later, you will always approach the Golden Ratio in the same way. It doesn’t matter what your starting numbers are. In Lucas’s sequence, the numbers don’t even have to begin in order. So let’s say I have:

K[0] := 2
K[1] := 5

K[n_] := K[n-2] + K[n-1]

So, just for kicks, I’ll show the first 31 terms:

Table[K[n], {n, 0, 30}]
{2, 5, 7, 12, 19, 31, 50, 81, 131, 212, 343, 555, 898, 1453, 2351, \
3804, 6155, 9959, 16114, 26073, 42187, 68260, 110447, 178707, 289154, \
467861, 757015, 1224876, 1981891, 3206767, 5188658}

Now, let’s output the Golden Ratio to 15 decimals as a reference:

N[GoldenRatio, 15]
1.61803398874989

Now, let’s take the ratio of the last two numbers in my 31-member sequence:

N[K[30]/K[29], 15]
1.61803398874942

You may say that the last two digits are off, but trying against the Fibonacci sequence, the ratio of the 30th and 31st numbers yields merely: 1.61803398874820, off by 3 digits.

For Lucas: 1.61803398875159, off by 4 digits — even worse.

So, my made-up sequence is more accurate for \phi than either Lucas or Fibonacci. I have tried other made-up sequences. Some are more, and some are less accurate. If it depends on the starting numbers, I think some combinations work better, and you won’t necessarily get greater accuracy by starting and ending with larger numbers.

Mathematica: Piecewise functions and shortcuts

I work with piecewise functions a lot in the courses I teach, and sometimes to do a quick reality check I would enter a piecewise function on Mathematica. Mathematica (the old version 5 I have) can do piecewise functions by declaring a function beginning with a replaceable parameter:

f[x_]:=

After this, I press <ESC>pw<ESC><CTRL+ENTER> and I get:

piece_startwith four placeholders for math expressions and their restrictions, which can be added by pressing CTRL+ENTER. I can add more by pressing CTRL+ENTER again. If I want an exponent, I type the base, then CTRL+6 then my exponent. Mathematica 5 makes it only partially clear what keyboard shortcuts to use for math expressions so that you don’t have to go to the pallette each time, but they are there.

Also the hints are not always there if you glide your mouse on the palette, so you will need to look for “keyboard shortcuts” in the documentation.

Here is a table of some shortcuts you might frequently use:

Character Keypress Comments
x3 Ctrl+6 Superscript or exponent
Ctrl+/ Fraction
square root Ctrl+2
x3 Ctrl+_ Subscript
Ctrl+Spacebar Moves cursor out of a formula by 1 level
Ctrl+Enter Adds another matrix row or creates a 2×2 if one does not exist
α ESC+a+ESC Greek alpha
β ESC+b+ESC Greek beta
π ESC+pi+ESC Greek pi
ESC+\infty+ESC One of many TEX-style ways of getting special characters
Δ ESC+D+ESC Greek capital delta

VIII: The Programmatic Side of Mathematica: Sampling With Replacement

Arrays. I can declare an array and return an arbitrary element from it:

S := {2,4,6,8,10}
S[[3]]
6

I can also do the same things for an anonymous array:

{2,4,6,8,10}[[3]]
6

In both cases, I request the third element of the array in double square brackets, and Mathematica returns with the number 6. I can request a random element as was shown in previour posts, but it would be better to name the array:

S[[Random[Integer,{1,Length[S]}]]]

The above statement returns a random element from an array without needing to know the length of the array. I could include it in my bag of tricks, since it would resemble sampling, especially if I sample a certain number of times.

Sampling with replacement means picking an element e from a set S without removing it. With sets, that amounts to just choosing elements from a set while allowing repetition. We don’t even need numbers. What about letters? Let’s re-define the set S thusly:

S:={a,b,c,d,e,f,g,h,i,j,k,l}

Now, let’s define a function called “pick[n]” which chooses some quantity of letters from S by passing the number of choices to it as “n”:

pick[n_]:=Table[S[[Random[Integer,{1,Length[S]}]]],{n}]

That’s a lot of nested brackets. But recall that the Table[] function generated a set of size n by executing the function in the first parameter n times, and returning the result as an array. But here, the function is passed into the Random[] function, which determines the array index of the element chosen from S. So, two calls to pick returned these for me:

pick[4]
{b,e,a,d}
pick[4]
{h,g,h,d}

So, as you can see, we have repetition.

VI: The Programmatic Side of Mathematica: The perfect shuffle of a deck of cards

Making a deck of cards requires making a the basic thirteen card values:

Join[Range[2,10], {J, Q, K, A}]
{2,3,4,5,6,7,8,9,10,J,Q,K,A}

Then, you need to distribute these values among the four suits:

Outer[List, {c, d, h, s}, %]
{{{c,2},{c,3},{c,4},{c,5},{c,6},{c,7},{c,8},{c,9},{c,10},{c,J},{c,Q},{c,K},{c,
   A}},{{d,2},{d,3},{d,4},{d,5},{d,6},{d,7},{d,8},{d,9},{d,10},{d,J},{d,
   Q},{d,K},{d,A}},{{h,2},{h,3},{h,4},{h,5},{h,6},{h,7},{h,8},{h,9},{h,
   10},{h,J},{h,Q},{h,K},{h,A}},{{s,2},{s,3},{s,4},{s,5},{s,6},{s,7},{s,
   8},{s,9},{s,10},{s,J},{s,Q},{s,K},{s,A}}}

This set gives four subsets representing each suit, then expresses each value-suit combination as a set. You want to keep the value-suit combinations, but if you want to shuffle the deck, you need to “flatten” the set somewhat to remove the set boundaries for the suits. That would be a call to the Flatten[] function with a “1” as a second parameter, indicating that we know that this set is a superset of four subsets, and we would like to remove the boundaries for the four subsets (the suits) and create a superset of the 52 individual cards:

Flatten[%,1]
{{c,2},{c,3},{c,4},{c,5},{c,6},{c,7},{c,8},{c,9},{c,10},{c,J},{c,Q},{c,K},{c,
  A},{d,2},{d,3},{d,4},{d,5},{d,6},{d,7},{d,8},{d,9},{d,10},{d,J},{d,Q},{d,
  K},{d,A},{h,2},{h,3},{h,4},{h,5},{h,6},{h,7},{h,8},{h,9},{h,10},{h,J},{h,
  Q},{h,K},{h,A},{s,2},{s,3},{s,4},{s,5},{s,6},{s,7},{s,8},{s,9},{s,10},{s,
  J},{s,Q},{s,K},{s,A}}

The cards are in numerical order of clubs, diamonds, hearts and spades. We would now like to shuffle the deck. What would a good dealer do to shuffle a deck? I would guess that he would first split the deck:

Partition[%,26]
{{{c,2},{c,3},{c,4},{c,5},{c,6},{c,7},{c,8},{c,9},{c,10},{c,J},{c,Q},{c,K},{c,
    A},{d,2},{d,3},{d,4},{d,5},{d,6},{d,7},{d,8},{d,9},{d,10},{d,J},{d,
    Q},{d,K},{d,A}},{{h,2},{h,3},{h,4},{h,5},{h,6},{h,7},{h,8},{h,9},{h,
    10},{h,J},{h,Q},{h,K},{h,A},{s,2},{s,3},{s,4},{s,5},{s,6},{s,7},{s,
    8},{s,9},{s,10},{s,J},{s,Q},{s,K},{s,A}}}

That creates that three-layered set of sets of sets again, but this time the middle layer is in two pieces, where the deck was split, rather than being split at each suit. Now, what? Our dealer would likely shuffle the deck. If he is skilled he could do what is known as a perfect shuffle, where the cards from each half-deck are interleaved with cards in the other half:

Transpose[%]
{{{c,2},{h,2}},{{c,3},{h,3}},{{c,4},{h,4}},{{c,5},{h,5}},{{c,6},{h,6}},{{c,
   7},{h,7}},{{c,8},{h,8}},{{c,9},{h,9}},{{c,10},{h,10}},{{c,J},{h,J}},{{c,
   Q},{h,Q}},{{c,K},{h,K}},{{c,A},{h,A}},{{d,2},{s,2}},{{d,3},{s,3}},{{d,
   4},{s,4}},{{d,5},{s,5}},{{d,6},{s,6}},{{d,7},{s,7}},{{d,8},{s,8}},{{d,
   9},{s,9}},{{d,10},{s,10}},{{d,J},{s,J}},{{d,Q},{s,Q}},{{d,K},{s,K}},{{d,
   A},{s,A}}}

All of these commands can be nested inside of one command so that we get:

Transpose[Partition[Flatten[Outer[List, {c, d, h, s},
    Join[Range[2,10], {J, Q, K, A}]], 1], 26]]

Oh yeah, and don’t forget to flatten the set again to remove the split in the deck:

Flatten[Transpose[Partition[Flatten[Outer[List, {c, d, h, s},
    Join[Range[2,10], {J, Q, K, A}]], 1], 26]], 1]
{{c,2},{h,2},{c,3},{h,3},{c,4},{h,4},{c,5},{h,5},{c,6},{h,6},{c,7},{h,7},{c,
    8},{h,8},{c,9},{h,9},{c,10},{h,10},{c,J},{h,J},{c,Q},{h,Q},{c,K},{h,K},{c,
    A},{h,A},{d,2},{s,2},{d,3},{s,3},{d,4},{s,4},{d,5},{s,5},{d,6},{s,6},{d,
    7},{s,7},{d,8},{s,8},{d,9},{s,9},{d,10},{s,10},{d,J},{s,J},{d,Q},{s,Q},{d,
    K},{s,K},{d,A},{s,A}}

I wouldn’t personally use these to play cards at this point, since the cards are still shuffled predictably and reproducibly. In the real world, if it so happened that a dealer shuffled the deck by interleaving the cards, this might still be the result from a new deck. But what you want is for them to be shuffled randomly and not predictably. You may even want to store the array in a variable so you can’t see the deck right away.

So, what is known as the “perfect shuffle” is only “perfect” in the sense that the interleaving of two halves of an ordered deck is done without errors. If all interleavings are perfect beginning with an ordered deck, I would guess that all subsequent deck splittings and interleavings are also repreoducible, and therefore predictable, assuming all shuffles are perfect shuffles. In the real world, random errors would be introduced in the interleavings (two cards slip by on one hand before one slips by on the other for example), which introduce the element of chance in card games.

The Programmatic Side of Mathematica IV: Functions Associated with Pascal’s Triangle

The nth triangular number can be written as a function:

triangular[n_] := Sum[i, {i, 1, n}]

What about a series of such numbers? You may sequentialize triangular[n] and make it return an array:

Table[triangular[n], {n, 1, 100}]

But it’s better to give an arbitrary length array:

triarr[n_] := Table[triangular[m], {m, 1, n}]

so that we can do things like obtain the first 25 triangular numbers:

triarr[25]

{1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210,231,253,276,\
300,325}

And while we’re at it, we can declare another function for tetrahedral numbers:

tetrahedral[m_] := Apply[Plus, triarr[m]]

which gives the nth tetrahedral number. And this gets us into the world of something called “n-Simplex” numbers. All of the n-Simplex numbers can be illustrated on Pascal’s triangle. Pascal’s triangle itself can be generated by masssaging a Mathematica function, although the shape of the triangle is not the familiar isosceles one:

pascal[p_] := TableForm[Table[Binomial[n, k], {n, 0, p}, {k, 0, n}]]

pascal[10]

But it is still kind of OK, since now we have the triangular numbers arranged in a column instead of a diagonal. In this case, the third column from the left. The tetrahedral numbers are in the fourth column, the 4-Simplex numbers (alternatively called the pentatopic or pentachron numbers) in the fifth column, the 5-Simplex numbers in the sixth column and so on. Like the tetrahedral numbers which arise from the sums of the first n triangular numbers, the nth k-simplex numbers generally are the sums of the first n (k-1) simplex numbers.

I imagine that the reason the numbering of “k-simplex” numbers are one behind their popular names (pentatopic seems to suggest “5” in the name, yet it is a 4-simplex set of numbers) is because the counting starts from 0 in column 1. While I haven’t read much on this, I am guessing that the “0-Simplex” numbers must be a repeating set of 1’s going forever, like we have in the first column, which I ought to be calling “column zero,” since it matches all Pascal’s triangle members with their combinations formulae. The “1-Simplex’ numbers are just the set of natural numbers in sequence from 1 to infinity. The “2-Simplex” numbers are the triangular numbers; and the “3-Simplex” numbers are the tetrahedral numbers. This brings us back to the pentatopic “4-Simplex” numbers. It seems fitting, anyway.

The Simplex system is supposed to also relate to geometry such that 0-Simplex means 0 dimensions or a point. 1-Simplex can be represented by a line, or points joined together in one dimension. 2-Simplex, our triangular numbers, fit in two dimensions; 3-simplex, the tetrahedral numbers, fit in 3-dimensions. The triangle is a shape that has the fewest sides that can fit in 2 dimensions, while the tetrahedron is a shape that has the fewest sides that can fit in 3 dimensions. After that, we are discussing shapes with “hyperfaces” in hyperdimensional space. The rest of this discussion can be found on Wikipedia.

You might also recall that each row in Pascal’s triangle represents the coefficients of the binomial expansion of (a + b)n.

binexpand[z_] := Expand[(a + b)^z]

binexpand[7]

whose output becomes: a^7 + 7a^6b + 21a^5b^2 + 35a^4b^3 + 35a^3b^4 + 21a^2b^5 + 7ab^6 + b^7, the binomial expansion of (a + b)^7.

 

The Programmatic Side of Mathematica II

Since the last time I posted, I found a better way to do 100 dice rolls:

This gives 100 dice rolls in a tabulated format in a much simpler syntax. You can also make a “die” out of complex numbers, but it is a little more of an issue to truncate:

… which gives output like:

Strictly speaking, these are all a kind of complex number called a pure imaginary number, but I guess I can’t split hairs.

You can create a list of divisors of integers like 28, like so:

Divisors[28]
{1, 2, 4, 7, 14, 28}

28 is a called a perfect number, because like 6, the number can be created by the sum of its proper factors. A “proper factor” of a number k does not include k itself, so if you add the numbers generated in the above list, you will obtain 2k, or 56.

So to get the correct result, you can make your own function that can re-generate the list with the last element truncated. The Drop[] function can do this:

I have declared a small function called ProperFactors, and the underscore after the k indicates that it is a placeholder for a parameter. When you run this, it returns nothing, but of course all you really did was declare a function. The Drop[] function drops elements from the list generated by the Divisors[] function, beginning at the last element (the reason there is a “-1” passed to Drop[]. In the list generated before, the last element was the highest factor of 28, which was 28 itself. The remaining elements in the list are the proper factors.

So, the call:

ProperFactors[28]

generates the list:

{1, 2, 4, 7, 14}

which are the proper factors of 28. Let’s programmatically check to see if the sum of the proper factors add to 28. Of course, checking by adding the numbers up yourself will tell you this: 1 + 2 + 4 + 7 + 14 = 28. But making a function do the work will be beneficial if the factors become very large. So, we can make use of the Apply[] function:

Apply[Plus, {1,2,3,4,5}]

This returns:

15

or the sum of 1 + 2 + 3 + 4  + 5. It takes a list of numbers (an actual array), and treats them as though each element were parameters passed into the Plus[] function. I could get the same thing with:

Plus[1, 2, 3, 4, 5]
15

The curly braces aren’t needed, since Plus[] expects a list of arbitrary length, but not an array. Since ProperFactors[] generates an array, I need to convert the array into a list to be operated on by Plus[]. This is the reason for the Apply[] function.

List[] is a function which generates an array:

List[1,2,3,4,5]
{1, 2, 3, 4, 5}

And simply invoking the literal array {1, 2, 3, 4, 5} does the same thing as List[]. And you can assign the array to a variable.

s := List[2, 3, 4, 5, 6]

does the same thing as

s:= {2, 3, 4, 5, 6}

so that

Apply[Plus, s]

converts the array s to a list of parameters passed to Plus. The result is:

20

Now, let’s go back to our perfect number checking function. Let’s call it PerfectCheck[]:

PerfectCheck[k_] := Apply[Plus, ProperFactors[k]] == k

What does “Apply[Plus, ProperFactors[k]]” do? Let’s invoke that by itself in Mathematica to see, using k = 28:

Apply[Plus, ProperFactors[28]]
28

so, it adds up the proper factors of 28. To make sure, let’s use a non-perfect composite number:

Apply[Plus, ProperFactors[100]]
117

To check this, go further back and list the proper factors of 100 themselves:

You will see when adding by hand that these add to 117.

Thus, a PerfectCheck[100] would be false since 117 <> 100.

PerfectCheck[100]
false
PerfectCheck[28]
true
PerfectCheck[6]
true
PerfectCheck[16]
false

Checking numbers in the range 1 to n can be done using Select[]:

perfect[n_] := Select[Range[n], PerfectCheck]

So that:

perfect[10000]

yields:

{6, 28, 496, 8128}

I was trying to find a fifth perfect number, but it’s difficult. I even cranked up the call to perfect[10000000] but no luck, with Mathematica chugging away for about 10 minutes. But I was able to find information through Google more quickly. According to The Math Forum, they say that the fifth one was still some ways off: 33,550,336.

A sixth number may exist, but as you can see the trend seems to indicate that they are pretty sparse on the number line. The Math forum hinted that these perfect numbers obey a number pattern according to the formula: (2n)(2n+1 – 1). This isn’t really reliable; that is to say that the value of integer used for n will not generate the nth perfect number necessarily. This seems to generate a perfect number so long as 2n+1 – 1 equals the value of a Mersenne Prime. A Mersenne Prime is a prime number that is one less than some power of 2, such as the number 7, which is 23 – 1. When finding the “highest known prime number” these days, mathematicians look for Mersenne Primes.

Using a command like

can fool you into thinking you have just printed out the first 20 perfect numbers, but you would be wrong. The list it gives out:

has numbers, such as 120, 496, and 2016, that are not perfect. You still need to check each one with our PerfectCheck[] function (either that, or fish out the numbers with Mersenne prime factors among the lot).

You would have to keep hunting until you find 8,589,869,056 as being the sixth perfect number, closer to the end of the list. Then the seventh one is two more over: 137,438,691,328.

But, of course we have a computer, so life gets happy once again for the slack among us. We can mash up the Table[] call above, and Select[] numbers that are perfect with PerfectCheck:

which generates the array:

and all of these are perfect numbers. I can go even more nutzo and increase the upper limit to 100:

… thereby generating the first 10 perfect numbers.

There are to date only 47 known Mersenne primes according to Wikipedia. Therefore, there are only 47  known perfect numbers, assuming all perfect numbers are even.

The programmatic side of Mathematica I

I have been playing with Mathematica 5, and have a big thick book on it from Wolfram and even a Schaum’s Outline on the program. All of these books touch on all of the functions of the program as though that was all there was to know.

But of course you can make your own user-defined functions, but material on that is harder to find. I managed to get my hands on a book published in the early 90s by Gaylord, Kamin, and Wellin, entitled Introduction to Programming with Mathematica. It does not appear that any specific version number for this program is mentioned in the text. But the publication date of 1993 would suggest that the intended version would be 2.1 or 2.2. I have version 5.2, but so far in the first few pages, I haven’t run into any trouble.

In the old days when an operating system was nothing more than a BASIC command interpreter (anyone old enough to remember the Commodore 64? or the TRS-80?), you are by default placed into an immediate mode where you can issue commands, test out command syntax or run programs written in the same language as understood by your interpreter (for the Commodore or TRS-80, that language would be BASIC). In Mathematica, that is what is happening. You are in a sort of “immediate mode” where you can issue commands or run pre-defined functions or programs, load your own functions, and so on.

The Mathematica kernel seems to allow you to enter similar commands as in this discussion without the need to press “shift-enter” to execute it as you do in the so-called “Notebook” interface. Just the “enter” key will suffice. I needed to open the kernel from the Start menu of Windows. Then I tried to test it out:

In[1]:= 39/7

        39
Out[1]= --
        7

OK, the interpreter works, and it works the way it seems to in the “Notebook” interface. You can also get help from a command by placing a question mark before it:

In[2]:= ?Plot
Plot[f, {x, xmin, xmax}] generates a plot of f as a function of x from xmin to xmax.
Plot[{f1, f2, ... }, {x, xmin, xmax}] plots several functions fi.

But the plotting of graphs  in kernel mode is a little … uh … different:

In[15]:= Plot[Sin[x], {x,0,2Pi}]

    #               ####                                                      
   1#            ####  #####                                                  
    #          ###         ###                                                
    #        ###             ###                                              
    #       ##                 ##                                             
    #     ###                   ##                                            
 0.5#    ##                      ##                                           
    #   ##                        ##                                          
    #  ##                          ###                                        
    # ##                             ##                                       
    ###                               ##                                      
    ### #  # # #  # # # #  # # #  # # ###  # # #  # # # #  # # #  # # # #  ###
  ############################################################################
    #          1           2          3  ##       4          5           6##  
    #                                     ##                             ##   
    #                                      ##                          ###    
    #                                       ##                        ##      
    #                                        ##                      ##       
-0.5#                                         ###                   ##        
    #                                           ##                ###         
    #                                            ###             ##           
    #                                              ###         ###            
    #                                                ####  #####              
  -1#                                                   ####

So, it’s probably better to use the Notebook mode for graphs…

Being a math program, you can create your own lists or generate new lists from Mathematica’s “Range” function:

Range[min,max,increment]

A “List” function can generate arbitrary lists:

In[16]:= List[1,2,3,-4,dog, cat, Sin, "Hello World!", {a, b}, Pi, {}]

Out[16]= {1, 2, 3, -4, dog, cat, Sin, Hello World!, {a, b}, Pi, {}}

In[17]:= Range[-12,15,7]

Out[17]= {-12, -5, 2, 9}

Notice the last output was for a Range command I inputted. The command asked for every 7th number between -12 and 15, including -12.

It doesn’t have loops in the pascal sense, but it has iterators, which are more compact. Here, I generate 50 random numbers:

In[18]:= Table[Random[],{50}]

Out[18]= {0.384617, 0.0481281, 0.278996, 0.161237, 0.880606, 0.587204, 0.996814, 0.0335705, 0.367147, 0.730021, 

>    0.195816, 0.707583, 0.316029, 0.030375, 0.971311, 0.607889, 0.172245, 0.461959, 0.511372, 0.001557, 0.830695, 

>    0.517233, 0.716287, 0.984279, 0.446078, 0.469105, 0.437292, 0.823042, 0.565472, 0.881901, 0.440478, 0.789472, 

>    0.198325, 0.151881, 0.244661, 0.0818886, 0.882296, 0.121506, 0.27335, 0.474, 0.71005, 0.659546, 0.761978, 

>    0.472443, 0.879355, 0.142313, 0.0456913, 0.488163, 0.433277, 0.673208}

{50} refers to the number of times the Random[] command must be repeated. Trouble is, the output looks a bit messy. You can clean this up by dividing the {50} into two dimensions, then tacking on a //TableForm modifier to the end of the whole command:

In[20]:= Table[Random[],{10},{5}]//TableForm
Out[20]//TableForm= 0.503921    0.225364   0.462814   0.985969   0.276704
                    0.613243    0.320156   0.961345   0.836671   0.519041
                    0.827408    0.256483   0.473772   0.283214   0.222593
                    0.39894     0.308651   0.763449   0.460018   0.900552
                    0.908033    0.286537   0.931285   0.096198   0.404112
                    0.0611721   0.468471   0.110229   0.127408   0.447929
                    0.148315    0.148884   0.290737   0.928888   0.320907
                    0.8924      0.816965   0.645673   0.098314   0.493461
                    0.508314    0.882225   0.638296   0.592909   0.600281
                    0.595688    0.707011   0.496711   0.196169   0.534516

This looks a bit better. {10},{5} specifies 10 rows and 5 columns. But of course you can have counters:

In[25]:= Table[i*j, {i, 1, 8}, {j, i, 8}]//TableForm
Out[25]//TableForm= 1    2    3    4    5    6    7    8
                    4    6    8    10   12   14   16
                    9    12   15   18   21   24
                    16   20   24   28   32
                    25   30   35   40
                    36   42   48
                    49   56
                    64

Here are the equivalent of two nested loops. You can think of “i” as the outer loop, and “j” as the inner loop, and I think the whole thing will make sense. And the loop syntax will accept a fourth parameter for loop increment:

In[35]:= Table[i*j, {i, 1, 8}, {j, i, 8, 2}]//TableForm
Out[35]//TableForm= 1    3    5    7
4    8    12   16
9    15   21
16   24   32
25   35
36   48
49
64

This one causes j to go to every second number. Note the change in the output. With a few more number operations, I can have 100 dice rolls in table format:

In[57]:= Table[IntegerPart[Random[]*6]+1, {10},{10}]//TableForm
Out[57]//TableForm= 6   3   6   6   4   4   5   2   5   1
                    5   5   3   2   3   1   2   1   5   6
                    2   2   6   1   2   5   6   2   4   2
                    2   1   6   1   3   2   3   5   6   2
                    1   5   2   2   6   4   3   1   5   5
                    3   6   1   3   1   6   1   3   4   4
                    5   4   4   3   4   6   3   2   4   3
                    6   1   6   4   4   2   6   1   3   3
                    5   5   6   6   6   1   2   3   3   2
                    6   2   5   6   6   1   5   2   2   5

IntegerPart[] takes the place of Trunc[] in many languages, returning the integer part of a number containing decimals. The chopping, or truncation of decimals is necessary to guarantee a dice roll between 1 and 6. Round[] is also available if you wish to experiment with it. It will round, but you will be rolling some zeroes as a result.