(*
 
Efficient Cycle Detection - for Collatz 3n+C

Seed7 version using Brent's Cycle Detection Algorithm

see Wikioedia article: http://en.wikipedia.org/wiki/Cycle_detection

ecd010.sd7

*)

$ include "seed7_05.s7i";
include "bigint.s7i";

#
# use an output structure compatible with other cycle detection programs
#
const type: cycle is new struct
    var bigInteger: the_attractor is 0_;  # smallest number in the cycle
    var bigInteger: the_initiator is 0_;  # seed number that lead to the cycle
    var bigInteger:       the_gcd is 0_;  # gcd of cycle numbers
    var bigInteger:    the_length is 0_;  # length of cycle (in odd numbers)
  end struct;

#
# record all attractors found (with their stats)
#
var array cycle: cycle_stats is 0 times cycle.value;

#
# but only track unique ones, so keep a set of attractors found and
# only add a cycle_stats record if attractor is unique
#
var set of bigInteger: att is (set of bigInteger).EMPTY_SET;

#
# the Collatz funtion (assumes n is odd): (3n+C)/2**LSB
#
const func bigInteger: collatz (in bigInteger: n, in bigInteger: C) is func
  result
    var bigInteger: cn is 0_;
  local
    var integer: LSB is 0;
  begin
    cn := n * 3_ + C;
    LSB := lowestSetBit(cn);
    cn >>:= LSB;               # remove factors of 2 in one fell swoop
  end func;

#
# once cycle is detected, make one sweep of the complete cycle
# to find the attractor and the cycle length
#
const func bigInteger: attractor (in bigInteger: n, in bigInteger: C, in bigInteger: seed) is func
  result
    var bigInteger:   attractor is 0_;
  local
    var integer:            over is 0;
    var bigInteger: cycle_length is 0_;
    var bigInteger:  entry_point is 0_;
    var bigInteger:    cycle_gcd is 0_;
    var bigInteger:            a is 0_;
    var cycle:   attractor_stats is cycle.value;
  begin
    attractor   := n;
    entry_point := n;                             # keep track of this so we know when cycle complete
    a := n;
    while over = 0 do
      if a < attractor then                       # is the new node smaller than current attractor?
        attractor := a;
      end if;
      a := collatz(a,C);
      cycle_length +:= 1_;
      if cycle_gcd = 0_ then                      # only do this once, on first cycle step
        cycle_gcd := gcd(a,entry_point);
      end if;
      if a = entry_point then
        over := 1;
      end if;
    end while;
    if attractor not in att then                  # record the stats (will be printed later)
      incl(att,attractor);                        # add it to the set (used to record unique attractors)
      attractor_stats.the_attractor := attractor;
      attractor_stats.the_initiator := seed;
      attractor_stats.the_gcd       := cycle_gcd;
      attractor_stats.the_length    := cycle_length;
      cycle_stats &:= [](attractor_stats);        # add attractor & stats to cycle_stats records
    end if;
  end func;

#
# Brent's modified 'tortoise & hare' strategy
#
const func bigInteger: brent (in bigInteger: n, in bigInteger: C) is func
  result
    var bigInteger: tortoise is 0_;
  local
    var integer:     power is 1;
    var integer:     lam   is 1;
    var bigInteger:  hare  is 0_;
    var integer:     i     is 1;
  begin
    tortoise := n;
    hare     := collatz(n,C);
    while tortoise <> hare do
      if power = lam then
        tortoise := hare;
        power *:= 2;
        lam := 0;
      end if;
      hare := collatz(hare,C);
      lam +:= 1;
    end while;
    tortoise := n;
    hare     := n;
    for i range 0 to lam-1 do
      hare := collatz(hare,C);
    end for;
    while tortoise <> hare do
      tortoise := collatz(tortoise,C);
      hare     := collatz(hare,C);
    end while;
  end func;

const proc: main is func
local
  var bigInteger:          n   is 0_;   # MUST be odd!
  var bigInteger:          e   is 0_;
  var bigInteger:          C   is 0_;
  var bigInteger:          o   is 0_;
  var cycle:               a   is cycle.value;
  var bigInteger:       seed   is 0_;
begin
  if (length(argv(PROGRAM)) < 4) then
    writeln();
    writeln("ABORTING -- requires 4 parameters");
    writeln();
    writeln(" usage: ecd010 N E C O");
    writeln("        N - start             of seed range (must be odd)");
    writeln("        E - end               of seed range");
    writeln("        C - C                 of the 3n+C Collatz function");
    writeln("        O - options           none as of this version (enter 0)");
    writeln();
    writeln("EXAMPLE: $ ./ecd010 1 40 16777213 0");
    writeln("                    | |  |        |");
    writeln("                    | |  |        no options");
    writeln("                    | |  use 3n+16777213");
    writeln("                    | up to 40");
    writeln("                    start @ 1");
    writeln();
    writeln("EXAMPLE OUTPUT: 16777213        16777213        16777213        1");                     
    writeln("                1       1       1       1");   
    writeln("                143     3       1       254252");   
    writeln("                29      5       1       254241");   
    writeln("                119     7       1       254214");   
    writeln("                |       |       |       |");   
    writeln("                |       |       |       length of cycle");   
    writeln("                |       |       gcd of cycle");   
    writeln("                |       seed that led to this cycle");   
    writeln("                attractor (smallest node in cycle)");   
    writeln();
    writeln("C is always added to the attractor list even if outside the requested seed range");
    writeln("because EVERY 3n+C has C as an attractor whose cycle stats are always C C C 1");
    writeln();
    exit(PROGRAM);
  end if;
  n := bigInteger parse (argv(PROGRAM)[1]);
  e := bigInteger parse (argv(PROGRAM)[2]);
  C := bigInteger parse (argv(PROGRAM)[3]);
  seed := n;
  #
  # enter C as the first attractor because 3n+C always has C as an attractor,
  # always has a gcd of C, and always has a length 1 (sv=[1])
  #
  incl(att,C);
  a.the_attractor := C;
  a.the_initiator := C;
  a.the_gcd       := C;
  a.the_length    := 1_;
  cycle_stats &:= [](a);
  #
  # loop though seed values from n to e, some C have attractors > C
  # but can be run in stages (which may give you duplicate attractors
  # with different seeds, but they would be the same (and would be suppressed
  # if both occur in the seed range)
  #
  while seed < e do
    if seed <> C then
      o := brent(seed,C);        # find a number in the cycle using Brent
      o := attractor(o,C,seed);  #  and then do the complete cycle (and collect stats)
    end if;
    seed +:= 2_;                 # seeds must always be odd
  end while;
  for a range cycle_stats do     # all attractors found, print them out
    write(a.the_attractor);
    write("\t");
    write(a.the_initiator);
    write("\t");
    write(a.the_gcd);
    write("\t");
    writeln(a.the_length);
  end for;
end func;