(*

Efficient Cycle Detection - for Collatz 3n+C

Seed7 version using Sedgewick's Cycle Detection Algorithm

scd002.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;

#
# table used by Sedgewick to detect cycles
# although fixed length, the length is set dynamically via command
# line parameter
#
const type: table is new struct
    var bigInteger: the_value is 0_;
    var integer:    the_index is -1;
end struct;

var array table: table_seq is 0 times table.value;

var integer:      position is 0;
var integer:             b is 1;
var integer:             g is 0;
var bigInteger:          C is 0_;
var integer:       counter is 0;
var integer:             M is 0;
var integer:           Max is 0;
var bigInteger:  seq_point is 0_;

#
# 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: SEED, in bigInteger: C) 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 attractor 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;

#
# BRandom is the function in which Sedgewick trys to find a cycle
# It's the same as collatz, only this modifies the global variable seq_point instead
# of returning a value as the function attracor expects
#
const proc: BRandom (inout bigInteger: seq_point, in bigInteger: C) is func
local
  var integer: LSB is 0;
begin
  seq_point := (seq_point * 3_) + C;
  LSB := lowestSetBit(seq_point);
  seq_point >>:= LSB;
end func;

#
# Bpurge used by Sedgewick to free space in the table
#
const proc: Bpurge () is func
local
  var integer: p is 0;
begin
  for p range 0 to M-1 do
    if ((table_seq[p+1].the_index mod (2*b)) > 0) then
      table_seq[p+1].the_index := -1;
    end if;
  end for;
  position := 0;
end func;

#
# Bsearch used by Sedgewick to search for...whatever it is it's looking for
#
const func integer: Bsearch () is func
  result
    var integer: found is -1;
  local
    var integer:     i is 0;
begin
  while (((i+1)<=M) and (found = -1)) do
    if ((i+1) > M) then
      writeln("Bsearch index error");
    end if;
    if (table_seq[i+1].the_value = seq_point) then
      found := i;
    else
      i +:= 1;
    end if;
  end while;
end func;

#
# Binsert adds the seq_point to the table
#
const proc: Binsert (in integer: i) is func
begin
  if ((position+1) > M) then
    writeln("Binsert1 index error");
  end if;
  while (table_seq[position+1].the_index <> -1) do
    position +:= 1;
    if ((position+1) > M) then
      write(" Binsert2 index error ");
      write(position);
      write(" ");
      write(i);
    end if;
  end while;
  table_seq[position+1].the_value := seq_point;
  table_seq[position+1].the_index := i;
  position +:= 1;
end func;

#
# Sedgewick's cycle detetion algorithm
# uses limited memory as opposed to Brent's memoryless algorithm
# hopefully, this one's faster
#
const proc: sedgewick (in integer: g, in bigInteger: SEED, in bigInteger: C, inout bigInteger: seq_point) is func
  local
    var integer:         i is 0;
    var integer:         m is 0;
    var integer:     found is -1;
    var integer: was_found is 0;
    var bigInteger:      a is 0_;
  begin
    while (found = -1) and (counter < Max) do
      #
      # Seed7 defaults to array indexing by 1 (grr...)
      # can be re-cast to index by 0, but I've already worked around it
      #
      if (((i mod b) = 0) and (m = M+1)) then
        Bpurge();
        b *:= 2;
        m := m div 2;
        i := 0;
      end if;
      #
      # Su's original program did not check for array bounds violations
      #
      if (position >= M) then
        Bpurge();
      end if;
      if (((i) mod b) = 0) then
        Binsert(i);
        m +:= 1;
      end if;
      BRandom(seq_point,C);
      i +:= 1;
      if (((i) mod (g*b)) < b) then
        found := Bsearch();
      end if;
      counter +:= 1;
    end while;
    if (counter = Max) then
      writeln("detection iteration counter exceeded");
    else
      a := attractor(seq_point,SEED,C);
    end if;
  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 bigInteger:         SEED is 0_;
  var cycle:               a   is cycle.value;
  var integer:             f   is 0;
  var table:               t   is table.value;
  var integer:             ii  is 0;
begin
  if (length(argv(PROGRAM)) < 7) then
    writeln();
    writeln("ABORTING -- requires 7 parameters");
    writeln();
    writeln(" usage: scd002 N E C O M G X");
    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("        M - table size        needs to be tuned for optimum performance (try > cycle_length/2)");
    writeln("        G - G                 Sedgewick scan parameter (needs tuning, try > M/2)");
    writeln("        X - iteration limit   when to give up (must be larger than cycle_length)");
    writeln();
    writeln("EXAMPLE: $ ./scd002 1 40 16777213 0 1000000 409600 999999");
    writeln("                    | |  |        | |       |      |");
    writeln("                    | |  |        | |       |      iteration limit");
    writeln("                    | |  |        | |       scan rate");
    writeln("                    | |  |        | use a table size of 1 million");
    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]);
  o   := bigInteger parse (argv(PROGRAM)[4]);
  M   := integer    parse (argv(PROGRAM)[5]);
  g   := integer    parse (argv(PROGRAM)[6]);
  Max := integer    parse (argv(PROGRAM)[7]);

  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 program 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)
  #

  table_seq := M times table.value;

  while SEED < e do
    if SEED <> C then
      seq_point := SEED;
      sedgewick(g,SEED,C,seq_point);        # find a number in the cycle using Sedgewick
    end if;
    SEED +:= 2_;                            # seeds must always be odd
    for ii range 1 to M do                  # clear table for next seed
      table_seq[ii].the_value := 0_;
      table_seq[ii].the_index := -1;
    end for;
    position := 0;                          # reset global variables for next seed
    b := 1;
    counter := 0;
  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;