Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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 Jul 31 '13 at 11:41
    
@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 Jul 31 '13 at 12:07
1  
I posted this link recently. Might be of interest... –  cormullion Jul 31 '13 at 12:22
3  
The 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 Jul 31 '13 at 14:58
1  
Take a look at Leonid's great answer in this question, it will give you an idea how to improve your code. –  sebhofer Aug 3 '13 at 8:26
show 5 more comments

2 Answers

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.

share|improve this answer
3  
Although 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. Aug 1 '13 at 16:07
add comment
up vote 0 down vote accepted

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.