(* 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;