home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / EXAMPLES.PAK / COLLATZ.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  4.3 KB  |  149 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3. (* RCS test line *)
  4.  
  5. (*:Name: NumberTheory`Collatz` *)
  6.  
  7.  
  8. (*:Title: The Collatz problem, Also Known as The 3x+1 Problem *)
  9.  
  10.  
  11. (*:Author: Ilan Vardi
  12.  
  13. *)
  14.  
  15. (*:Keywords: Collatz problem, 3x+1 problem
  16. *)
  17.  
  18. (*:Requirements: none. *)
  19.  
  20. (*:Warnings: The Collatz map is taken to be x -> x/2 if x is even and
  21.              x -> (3x+1)/2 if x is odd. This simplifies the analysis
  22.              of the problem. The previous package Collatz.m used the
  23.              map x -> 3x+1 if x is odd. 
  24.            
  25.              A table is computed when the package is loaded in, so it 
  26.              could take as long as 5 minutes to read in the package.
  27.  
  28.              $RecursionLimit is set to Infinity, since there is a high
  29.              degree of recursion in the algorithm.
  30.  
  31. *)
  32.  
  33. (*:Source: J.C. Lagarias, The 3x+1 problem and its generalizations,
  34.            American Math. Monthly #92 (1985), 3-23 
  35.  
  36.            I. Vardi, Computational Recreactions in Mathematica,
  37.            Addison-Wesley 1991, Chapter 7
  38.  
  39. *)
  40.  
  41. (*:Limitations: The program terminates only if an iterate reaches one of
  42.            the 4 known cycles. If the input reaches another cycle, then
  43.            the algorithm will not terminate.
  44.  
  45.            TotalStoppingTime seems to work up to 30,000 digit numbers
  46.            but runs out of memory after that.
  47.              
  48. *)
  49.  
  50. (* :Discussion: This package computes the iterates of the Collatz map 
  51.             x -> x/2 if x is even, x -> (3x+1)/2 if x is odd, until an
  52.             iterate reaches one of the 4 known cycles (the program runs
  53.             on positive and negative integers):
  54.  
  55.             {1,2}, {-1}, {-5, -7, -10}, 
  56.             {-17, -25, -37, -55, -82, -41, -61, -91, -136, -68, -34}
  57.  
  58.             An efficient algorithm is used to compute how many iterations
  59.             there are up to a cycle (the total stopping time). This
  60.             algorithm is discussed in detail in Computational Recreations 
  61.             in Mathematica, Chapter 7. 
  62.             
  63. *)
  64.  
  65.  
  66. BeginPackage["Examples`Collatz`"] 
  67.  
  68. Collatz::usage = "Collatz[n] returns the iterates of the Collatz map up
  69. to one of the known cycles."
  70.  
  71. TotalStoppingTime::usage = "TotalStoppingTime[n] gives the number of iterates
  72. of the Collatz map up to one of the known cycles. It gives the same result as
  73. Length[Collatz[n]]-1 but is much more efficient. TotalStoppingTime[n]  
  74. returning an answer gives a check of the 3x+1 conjecture for n."
  75.  
  76. Begin["`Private`"]  
  77.  
  78. (* V1.2   FoldList[f_, x_, list_]:= Accumulate[f, Prepend[list, x]] *)
  79.  
  80. CollatzT[n_]:= If[EvenQ[n], n/2, (3 n + 1)/2]
  81.  
  82. Collatz[1]:= {1} 
  83.  
  84. Collatz[-1]:= {-1} 
  85.  
  86. Collatz[-5]:= {-5}
  87.  
  88. Collatz[-17]:= {-17}
  89.  
  90. Collatz[n_Integer]:= Prepend[Collatz[CollatzT[n]], n] /; 
  91.                      Abs[n] < 2050
  92.  
  93.  
  94. Collatz[n_Integer]:= 
  95.  Block[{cr}, 
  96.         cr = Flatten[NestList[CollatzT, #, 9]& /@
  97.                      CollatzIterate[n]]; 
  98.        Join[cr, Rest[Collatz[Last[cr]]]]]
  99.  
  100. TotalStoppingTime[n_Integer]:= 
  101.   Length[Collatz[n]] - 1     /; Abs[n] < 2055
  102.  
  103. TotalStoppingTime[n_Integer] := 
  104.   10 + TotalStoppingTime[{n, 1} . 
  105.   CollatzTable10[[1 + Mod[n, 1024]]]] /;
  106.             2055 <= Abs[n] <= 2^1002
  107.  
  108. TotalStoppingTime[n_Integer]:= 
  109.  1000 + 
  110.  TotalStoppingTime[
  111.  NestList[Block[
  112.         {a = CollatzTable10[[1+Mod[#[[2]], 1024]]]},
  113.         {SemiProduct[#[[1]], a], Mod[{#[[2]], 1} . a, 
  114.                  2^1000]}]&,
  115.            {{1, 0}, Mod[n, 2^1000]},
  116.             100] [[-1, 1]] . {n, 1}]      /; 
  117.                         Abs[n] > 2^1002
  118.  
  119. SemiProduct[{a_, b_}, {c_, d_}]:= {a c, b c + d}
  120.  
  121.  
  122. CollatzIterate[n_Integer]:= {}     /; Abs[n] < 2055
  123.  
  124. CollatzIterate[n_Integer]:= 
  125.   Prepend[CollatzIterate[{n, 1} . 
  126.           CollatzTable10[[1 + Mod[n, 1024]]]], n] 
  127.  
  128.  
  129. CollatzTable[k_Integer]:= 
  130.   RotateRight[
  131.       Map[Function[x,  
  132.          {3^Apply[Plus, x] / 2^Length[x], 
  133.              Reverse[x] . 
  134.              FoldList[#1 #2 &, 1/2,
  135.                       1/2 3^Reverse[Rest[x]]]}],
  136.              Mod[Map[NestList[CollatzT,
  137.                          #, k-1] &,
  138.              Range[1, 2^k]], 2]]]
  139.  
  140. CollatzTable10 = CollatzTable[10] 
  141.  
  142. $RecursionLimit = Infinity;
  143.  
  144. End[]   (* Examples`Collatz`Private` *)
  145.  
  146. Protect[Collatz, TotalStoppingTime]
  147.  
  148. EndPackage[]   (* Examples`Collatz` *)    
  149.