Collatz Tool Box --- any speed ups possible?

4

3

I have created a tool box for experimenting with the Collatz conjecture. Info here: Wolfram MathWorld and here: Wikipedia The patterns I have found allow us to make the conjecture jump through some hoops.

Final Edit countOrbit has been re-cast as a Do Until loop because the logic behaves similar to a one-directional linked list and I don't think it can be improved.

Note: the number of sub orbits $\leq$ (total steps) / 6

Edit for the down-voters: the current maximum that collatz projects have reached is around $10^{18}$. The example below has two numbers $10^{301030}$ and $10^{477119}$ that my functions can use quickly and accurately. That is many magnitudes greater that what is currently being done. So, why the down-votes?

This is one example: m is the count of numbers in an ascending sequence. uniqueRank[1,m] returns the first number of the first unique sub orbit of that length. omegaSubOrbit[x,m] returns the last number of that sub orbit. We display the digit counts of these numbers. Then we count the complete orbit down to 1. This takes a bit over 12 minutes on my AMD 1100T. The left-hand number is the count of multiplies (we don't count divisions by two) and the right-hand count is the number of sub orbits used to get the counts.

m = 1000000; x = uniqueRank[1, m];
y = omegaSubOrbit[x, m];
{IntegerLength[x], IntegerLength[y]}
countOrbit[1, x]

{301030, 477119}
{4805005, 1903828}
{13420758, 1903828} 7.8 minutes to count all multiplies and divides

My question: Since I'm using For and While in a few places, is there a way to speed up those functions? Especially, the countOrbit function?

ascendingQ[x_] := 3 == Mod[x, 4]
uniqueQ[x_] := 0 != Mod[2 x - 1, 3]
subOrbit[x_] := 
 Block[{y = x}, 
Flatten[Join[
Table[y \[DirectedEdge] (y = 1/2 (3 y + 1)), {n, 1, 
  IntegerExponent[x + 1, 2] - 1}], {y \[DirectedEdge] (3 y + 1)/
   2^IntegerExponent[3 y + 1, 2]}]]]
createSubOrbitGraph[n_] :=
 Block[{j, v = {}},
  For[j = 1, j <= n, j++,
   If[uniqueQ[2 j + 1], AppendTo[v, subOrbit[2 j + 1]]];
   ];
  Flatten[v]
  ]
uniqueRank[n_, m_] := Block[{a = If[1 != n && OddQ[n], n - 1, n]},
   (n + IntegerExponent[a, 2]) 2^m - 1]  
nonUniqueRank[n_, m_] := (2 n - 1)   2^m - 1
createOrbit[w_, x_] := Block[{u = subOrbit[x]},
  While[w < u[[-1, 2]], 
   u = Flatten[AppendTo[u, subOrbit[u[[-1, 2]]]]]];
   u]
createOrbitGraph[n_] :=
 Block[{j, u = {}},
  For[j = 1, j <= n, j++,
   If[uniqueQ[2 j + 1], u = Union[u, createOrbit[2 j + 1, 2 j + 1]]];
   ];
  Flatten[u]
  ]
omegaSubOrbit[x_, m_] := 
 Block[{z = (-1 + (3/2)^m (1 + x))}, 2^(1 - IntegerExponent[2 z, 2]) z]
omegaSubOrbit[x_] := 
 Block[{m = IntegerExponent[x + 1, 2], z = (-1 + (3/2)^m (1 + x))}, 
  2^(1 - IntegerExponent[2 z, 2]) z]
countOrbit[w_, x_] := 
 Block[{h = x, t, c, d = 0, s = 0, l = 0, m, n, a},
  While[True,
   m = IntegerExponent[h + 1, 2];
   t = -1 + (3/2)^m (1 + h);
   n = IntegerExponent[t, 2];
   l += m;
   d += n;
   s += 1;
   h = t/2^n;
   If[h <= w, Break[]];
   ];
  c = 2 l + d;
  a = If[1 == w, "Full ", "Extended Sub "];
  Print["---Counts ", a, "Orbit---",
   "\nNumber of Sub Orbits: ", s,
   "\nMultiplications by 3: ", l,
   "\nDivisions by 2:       ", l,
   "\nExtra Divisions by 2: ", d,
   "\nTotal for Orbit:      ", c]
  ]    

Some brief documentation:
ascendingQ[x] Returns True if odd x is ascending, False if descending
uniqueQ[x] Returns True if odd x is not embedded in any longer sequence
subOrbit[x] Returns graph data in the form {Head[DirectedEdge]Tail}
createSubOrbitGraph[n] Returns graph data of sub orbits up through n
uniqueRank[n, m] Returns the first number of the n-th unique occurance of a sub orbit of size m
nonUniqueRank[n, m] Returns the first number of the n-th occurance of a sub orbit of size m, which may not be unique.
createOrbit[w, x] Returns graph data for w=1 complete orbit or w=x extended sub orbit
Note: for a graph, both are equivalent. w=x produces fewer duplicates.
createOrbitGraph[n] Returns graph data of extended sub orbits up through n
omegaSubOrbit[x] Returns the last odd number of a sub orbit. (The one that required multiple divides.)
omegaSubOrbit[x, m] Same as above, but uses m to speed up the first step
countOrbit[w, x] Returns the count of multiply and division steps and the count of sub orbits processed, for w=1 complete orbit or w=x extended sub orbit

Edit Replaced countOrbit[w,x] with refined version that counts divisions and multiplication. Reduced the time from 15 minutes to 7.8 minutes for the above example.

Fred Kline

Posted 2013-07-31T11:07:19.367

Reputation: 2 226

1

many more basic collatz ideas here & hope to hear from others in [chat]

– vzn – 2015-06-12T17:22:43.997

The Collatz toolbox is useful but it could be improved if it included optimized code for merging Collatz sequences to generate tree graphs of the sequences (ending in 4 -> 2 -> 1), with a few bells and options for coloring even and odd $n$ nodes, or placements based on iteration time, or $n$ value... – David G. Stork – 2017-10-10T16:22:45.200

I didn't downvote this, but it may be that some "speed up my code" questions aren't always well-received by everyone. relevant Meta question...

– cormullion – 2013-07-31T11:41:04.270

@cormullion, Because I'm an old procedural programmer, I feel most comfortable using For and While and I don't usually see how to eliminate them. This post is two-fold---One, to let others grab this stuff and play with it and Two, to find out how to speed up the countOrbit function to maybe get it below 10 minutes duration on my pc. – Fred Kline – 2013-07-31T12:07:37.783

@cormullion, and I had to show all the code because if you tweak one function, most likely you would have to tweak another to keep everything working together. – Fred Kline – 2013-07-31T12:10:47.827

1

I posted this link recently. Might be of interest...

– cormullion – 2013-07-31T12:22:37.540

I admit I am not familiar with this problem and do not really understand your code, and this understanding will certainly be necessary for anyone wishing to contribute any serious improvements. Do you really systematically check every number up to $10^{477119}$ in 12 minutes? – Oleksandr R. – 2013-07-31T12:46:22.050

@OleksandrR.,no, I use my short-cut, which is what this is all about. Notice these numbers from the example: $4805005, 1903828$. The first is the count of multiplies required to get down to $1$ and the second is the number of sub orbits processed to get to that number. Each sub orbit contains its own length, so the first sub orbit adds $1000000$ to the count and moves on to the next sub orbit. So, the whole thing is done with a series of additions. The key is: IntegerExponent[x+1,2] applied to the first number of the sub orbit provides that count. – Fred Kline – 2013-07-31T12:59:07.870

@OleksandrR., also the $10^{477119}$ number is calculated in omegaSubOrbit[x] with x being the uniqueRank number. Both calculations are very quick--under 2 seconds. – Fred Kline – 2013-07-31T13:05:51.527

This Q and A is better for learning how to eliminate inefficient procedural habits. – cormullion – 2013-07-31T14:39:15.660

4The reason for my downvoting this is: on the one hand the "speed up my code" reason that cormuillion mentioned, and that you give too many functions without explaining what they do to someone who doesn't know about the Collatz problem.

I am happy to retract the vote if you make this into a more minimal question with an explanation for what the conjecture is, why subOrbit[n] skips the even steps and some indication this is done correctly and so on.. – gpap – 2013-07-31T14:58:35.710

1

Take a look at Leonid's great answer in this question, it will give you an idea how to improve your code.

– sebhofer – 2013-08-03T08:26:17.307

Answers

2

As stated in the OP, this is a "linked-list" type of algorithm that I don't feel can be speeded up using the usual methods. So, I figured that compiling the functions to use Integers Of Unusual Size would be the only way to speed up things. As noted in an answer on that thread, Mathematica compiles are limited to machine-precision so, I'm debating with myself about crafting a C program using GMP to test some timings. I will accept this answer to clear this thread from the unanswered queue.

Fred Kline

Posted 2013-07-31T11:07:19.367

Reputation: 2 226

2

Simple speeds ups start with simple coding. Get rid of the modulus. Get rid of all division. If dividing by 2, 4, etc. use bit shifts. If diving by a fraction such as 1/2 multiple by the decimal instead (4 * 0.5 is usually faster than 4 * 1 / 2 as it requires fewer operations). Get rid of multiplication by 3 and use a left shift and an add ( x * 3 is usually slower than x<<1 + x. The above all depend upon how good the compiler is at optimizing the code but it usually works best when the programmer doesn't assume the compiler is smarter than he is.

user8827

Posted 2013-07-31T11:07:19.367

Reputation: 21

3Although I can see your point, this is not necessarily good advice for a very high level, interpreted language such as Mathematica. Division does not take appreciably longer than any other operation for large arbitrary-precision numbers, and nor is bit shifting normally faster than multiplication. Fundamentally, I disagree with the premise--these approaches, while they make sense in the right context, are not particularly simple, while speed-ups arguably start with performance profiling. – Oleksandr R. – 2013-08-01T16:07:51.880