Cycle detection

From Rosetta Code
Revision as of 12:53, 28 May 2016 by Rdm (talk | contribs) (J: show a "more traditional" representation of the 1 plus x squared polynomial)
Cycle detection is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Detect a cycle in an iterated function using Brent's algorithm.


Detecting cycles in iterated function sequences is a sub-problem in many computer algorithms, such as factoring prime numbers. Some such algorithms are highly space efficient, such as Floyd's cycle-finding algorithm, also called the "tortoise and the hare algorithm". A more time efficient algorithm than "tortoise and hare" is Brent's Cycle algorithm. This task will implement Brent's algorithm.

See https://en.wikipedia.org/wiki/Cycle_detection for a discussion of the theory and discussions of other algorithms that are used to solve the problem.

When testing the cycle detecting function, you need two things:

1) An iterated function

2) A starting value

The iterated function used in this example is: f(x) = (x*x + 1) modulo 255.

The starting value used is 3.

With these as inputs, a sample program output would be:

3,10,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5

Cycle length = 6

Start index = 2

The output prints the first several items in the number series produced by the iterated function, then identifies how long the cycle is (6) followed by the zero-based index of the start of the first cycle (2). From this you can see that the cycle is:

101,2,5,26,167,95


C_sharp

This solution uses generics, so may find cycles of any type of data, not just integers.

<lang csharp>

// First file: Cycles.cs // Author: Paul Anton Chernoch

using System;

namespace DetectCycles { 

 /// <summary>
 /// Find the length and start of a cycle in a series of objects of any IEquatable type using Brent's cycle algorithm.
 /// </summary>
 public class Cycles<T> where T : IEquatable<T>
 {
   /// <summary>
   /// Find the cycle length and start position of a series using Brent's cycle algorithm.
   /// 
   ///  Given a recurrence relation X[n+1] = f(X[n]) where f() has
   ///  a finite range, you will eventually repeat a value that you have seen before.
   ///  Once this happens, all subsequent values will form a cycle that begins
   ///  with the first repeated value. The period of that cycle may be of any length.
   /// </summary>
   /// <returns>A tuple where:
   ///    Item1 is lambda (the length of the cycle) 
   ///    Item2 is mu, the zero-based index of the item that started the first cycle.</returns>
   /// <param name="x0">First item in the series.</param>
   /// <param name="yielder">Function delegate that generates the series by iterated execution.</param>
   public static Tuple<int,int> FindCycle(T x0, Func<T,T> yielder)
   {
     int power, lambda;
     T tortoise, hare;
     power = lambda = 1;
     tortoise = x0;
     hare = yielder(x0);

     // Find lambda, the cycle length
     while (!tortoise.Equals (hare)) {
       if (power == lambda) {
         tortoise = hare;
         power *= 2;
         lambda = 0;  
       }
       hare = yielder (hare);
       lambda += 1;
     }

     // Find mu, the zero-based index of the start of the cycle
     var mu = 0;
     tortoise = hare = x0;
     for (var times = 0; times < lambda; times++) 
       hare = yielder (hare);
     
     while (!tortoise.Equals (hare)) 
     {
       tortoise = yielder (tortoise);
       hare = yielder (hare);
       mu += 1;
     }

     return new Tuple<int,int> (lambda, mu);
   }
 }

}

// Second file: Program.cs

using System;

namespace DetectCycles { class MainClass { public static void Main (string[] args) { // A recurrence relation to use in testing Func<int,int> sequence = (int _x) => (_x * _x + 1) % 255;

// Display the first 41 numbers in the test series var x = 3; Console.Write(x); for (var times = 0; times < 40; times++) { x = sequence(x); Console.Write(String.Format(",{0}", x)); } Console.WriteLine();

// Test the FindCycle method var cycle = Cycles<int>.FindCycle(3, sequence); var clength = cycle.Item1; var cstart = cycle.Item2; Console.Write(String.Format("Cycle length = {0}\nStart index = {1}\n", clength, cstart)); } } }

</lang>

Elixir

Translation of: Ruby

<lang elixir>defmodule Cycle_detection do 

 def find_cycle(x0, f) do
   lambda = find_lambda(f, x0, f.(x0), 1, 1)
   hare = Enum.reduce(1..lambda, x0, fn _,hare -> f.(hare) end)
   mu = find_mu(f, x0, hare, 0)
   {lambda, mu}
 end
 
 # Find lambda, the cycle length
 defp find_lambda(_, tortoise, hare, _, lambda) when tortoise==hare, do: lambda
 defp find_lambda(f, tortoise, hare, power, lambda) do
   if power == lambda, do: find_lambda(f, hare, f.(hare), power*2, 1),
                     else: find_lambda(f, tortoise, f.(hare), power, lambda+1)
 end
 
 # Find mu, the zero-based index of the start of the cycle
 defp find_mu(_, tortoise, hare, mu) when tortoise==hare, do: mu
 defp find_mu(f, tortoise, hare, mu) do
   find_mu(f, f.(tortoise), f.(hare), mu+1)
 end

end

  1. A recurrence relation to use in testing

f = fn(x) -> rem(x * x + 1, 255) end

  1. Display the first 41 numbers in the test series

Enum.map_reduce(0..40, 3, fn _,acc -> {acc, f.(acc)} end) |> elem(0) |> Enum.join(",") |> IO.puts

  1. Test the find_cycle function

{clength, cstart} = Cycle_detection.find_cycle(3, f) IO.puts "Cycle length = #{clength}\nStart index = #{cstart}"</lang>

Output:
3,10,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5
Cycle length = 6
Start index = 2

J

Brute force implementation:

<lang J>cdetect=:1 :0 

 r=. ~.@(,u@{:)^:_ y
 n=. u{:r
 (,(#r)-])r i. n

)</lang>

In other words: iterate until we stop getting new values, find the repeated value, and see how it fits into the sequence.

Example use:

<lang J> 255&|@(1 0 1&p.) cdetect 3 2 6</lang>

(Note that 1 0 1 are the coefficients of the polynomial 1 + (0 * x) + (1 * x * x) or, if you prefer (1 * x ^ 0) + (0 * x ^ 1) + (1 * x ^ 2) - it's easier and probably more efficient to just hand the coefficients to p. than it is to write out some other expression which produces the same result.)

Java

Works with: Java version 8

<lang java>import java.util.function.*; import static java.util.stream.IntStream.*;

public class CycleDetection { 

   public static void main(String[] args) {
       brent(i -> (i * i + 1) % 255, 3);
   }

   static void brent(IntUnaryOperator f, int x0) {
       int cycleLength;
       int hare = x0;
       FOUND:
       for (int power = 1; ; power *= 2) {
           int tortoise = hare;
           for (int i = 1; i <= power; i++) {
               hare = f.applyAsInt(hare);
                if (tortoise == hare) {
                   cycleLength = i;
                   break FOUND;
               }
           }
       }

       hare = x0;
       for (int i = 0; i < cycleLength; i++)
           hare = f.applyAsInt(hare);

       int cycleStart = 0;
       for (int tortoise = x0; tortoise != hare; cycleStart++) {
           tortoise = f.applyAsInt(tortoise);
           hare = f.applyAsInt(hare);
       }

       printResult(x0, f, cycleLength, cycleStart);
   }

   static void printResult(int x0, IntUnaryOperator f, int len, int start) {
       System.out.printf("Cycle length: %d%nCycle: ", len);
       iterate(x0, f).skip(start).limit(len)
               .forEach(n -> System.out.printf("%s ", n));
   }

}</lang> 

Cycle length: 6
Cycle: 101 2 5 26 167 95

ooRexx

<lang ooRexx>/* REXX */ x=3 list=x Do i=1 By 1 

 x=f(x)
 p=wordpos(x,list)
 If p>0 Then Do
   list=list x
   Say 'list ='list '...'
   Say 'Start index = ' (wordpos(x,list)-1) '(zero based)'
   Say 'Cycle length =' (words(list)-p)
   Say 'Cycle        =' subword(list,p,(words(list)-p))
   Leave
   End
 list=list x
 End

Exit f: Return (arg(1)**2+1)//255 </lang>

Output:
list =3 10 101 2 5 26 167 95 101 ...
Start index =  2 (zero based)
Cycle length = 6
Cycle        = 101 2 5 26 167 95

Perl 6

Works with: Rakudo version 2016-01

Pretty much a line for line translation of the Python code on the Wikipedia page.

<lang perl6>sub cyclical-function (\x) { (x * x + 1) % 255 };

my ( $l, $s ) = brent( &cyclical-function, 3 );

put join ', ', (3, -> \x { cyclical-function(x) } ... *)[^20], '...'; say "Cycle length $l."; say "Cycle start index $s."; say (3, -> \x { cyclical-function(x) } ... *)[$s .. $s + $l - 1];

sub brent (&f, $x0) { 

   my $power = my $λ = 1;
   my $tortoise = $x0;
   my $hare = f($x0);
   while ($tortoise != $hare) {
       if $power == $λ {
           $tortoise = $hare;
           $power *= 2;
           $λ = 0;
       }
       $hare = f($hare);
       $λ += 1;
   }

   my $μ = 0;
   $tortoise = $hare = $x0;
   $hare = f($hare) for ^$λ;

   while ($tortoise != $hare) {
       $tortoise = f($tortoise);
       $hare = f($hare);
       $μ += 1;
   }
   return $λ, $μ;

}</lang>

Output:
3, 10, 101, 2, 5, 26, 167, 95, 101, 2, 5, 26, 167, 95, 101, 2, 5, 26, 167, 95, ...
Cycle length 6.
Cycle start index 2.
(101 2 5 26 167 95)

Python

Function from the Wikipedia article: <lang python>import itertools

def brent(f, x0): 

   # main phase: search successive powers of two
   power = lam = 1
   tortoise = x0
   hare = f(x0)  # f(x0) is the element/node next to x0.
   while tortoise != hare:
       if power == lam:  # time to start a new power of two?
           tortoise = hare
           power *= 2
           lam = 0
       hare = f(hare)
       lam += 1

   # Find the position of the first repetition of length lam
   mu = 0
   tortoise = hare = x0
   for i in range(lam):
   # range(lam) produces a list with the values 0, 1, ... , lam-1
       hare = f(hare)
   # The distance between the hare and tortoise is now lam.

   # Next, the hare and tortoise move at same speed until they agree
   while tortoise != hare:
       tortoise = f(tortoise)
       hare = f(hare)
       mu += 1

   return lam, mu

def iterate(f, x0): 

   while True:
       yield x0
       x0 = f(x0)

if __name__ == '__main__': 

   f = lambda x: (x * x + 1) % 255
   x0 = 3
   lam, mu = brent(f, x0)
   print("Cycle length: %d" % lam)
   print("Cycle start index: %d" % mu)
   print("Cycle: %s" % list(itertools.islice(iterate(f, x0), mu, mu+lam)))</lang>
Output:
Cycle length: 6
Cycle start index: 2
Cycle: [101, 2, 5, 26, 167, 95]

A modified version of the above where the first stage is restructured for clarity: <lang python>import itertools

def brent_length(f, x0): 

   # main phase: search successive powers of two
   hare = x0
   power = 1
   while True:
       tortoise = hare
       for i in range(1, power+1):
           hare = f(hare)
           if tortoise == hare:
               return i
       power *= 2

def brent(f, x0): 

   lam = brent_length(f, x0)

   # Find the position of the first repetition of length lam
   mu = 0
   hare = x0
   for i in range(lam):
   # range(lam) produces a list with the values 0, 1, ... , lam-1
       hare = f(hare)
   # The distance between the hare and tortoise is now lam.

   # Next, the hare and tortoise move at same speed until they agree
   tortoise = x0
   while tortoise != hare:
       tortoise = f(tortoise)
       hare = f(hare)
       mu += 1

   return lam, mu

def iterate(f, x0): 

   while True:
       yield x0
       x0 = f(x0)

if __name__ == '__main__': 

   f = lambda x: (x * x + 1) % 255
   x0 = 3
   lam, mu = brent(f, x0)
   print("Cycle length: %d" % lam)
   print("Cycle start index: %d" % mu)
   print("Cycle: %s" % list(itertools.islice(iterate(f, x0), mu, mu+lam)))</lang>
Output:
Cycle length: 6
Cycle start index: 2
Cycle: [101, 2, 5, 26, 167, 95]

REXX

Brent's algorithm

<lang rexx>/*REXX pgm detects a cycle in an iterated function [F] using Brent's algorithm*/ init=3; @=init /* [↓] show a line of function values.*/ 

                do  until length(@)>79; @=@ f(word(@,words(@))); end; say @ ...

call Brent init /*invoke Brent algorithm for func F. */ parse var result cycle idx /*obtain the two values returned from F*/ say 'cycle length = ' cycle /*display the cycle. */ say 'start index = ' idx " ◄─── zero index" /* " " index. */ say 'the sequence = ' subword(@, idx+1, cycle) /* " " sequence.*/ exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ Brent: procedure; parse arg x0 1 tort; pow=1; #=1 /*tort is set to X0*/ hare=f(x0) /*get 1st value for the func. */ 

           do  while tort\==hare
           if pow==#  then do;  tort=hare      /*set value of TORT to HARE.  */
                                pow=pow+pow    /*double the value of  POW.   */
                                #=0            /*reset  #  to zero  (lambda).*/
                           end
           hare=f(hare)
           #=#+1                               /*bump the lambda count value.*/
           end   /*while*/

hare=x0 

           do #;          hare=f(hare)         /*generate number of F values.*/
           end   /*j*/

tort=x0 /*find position of the 1st rep*/ 

           do mu=0  while tort\==hare          /*MU  is a  zero─based  index.*/
                          tort=f(tort)
                          hare=f(hare)
           end   /*mu*/

return # mu /*────────────────────────────────────────────────────────────────────────────*/ f: return (arg(1)**2 + 1) // 255 /*this defines/executes the function F.*/</lang> output   using the defaults: 

3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95 101 ...
cycle length =  6
start index  =  2   (zero index)
start index  =  2   ◄─── zero index
the sequence =  101 2 5 26 167 95

sequential search algorithm

<lang rexx>/*REXX pgm detects a cycle in an iterated function [F] using sequential search*/ x=3; @=x 

                       do  until cycle\==0;   x=f(x)  /*calculate another num*/
                       cycle=wordpos(x,@)             /*is this a repeat ?   */
                       @=@ x                          /*append number to list*/
                       end   /*until*/
  1. =words(@)-cycle /*compute cycle length.*/

say ' original list=' @ ... /*display the sequence.*/ say 'numbers in list=' words(@) /*display number of #s.*/ say ' cycle length =' # /*display the cycle. */ say ' start index =' cycle-1 " ◄─── zero based" /* " " index. */ say 'cycle sequence =' subword(@,cycle,#) /* " " sequence.*/ exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ f: return (arg(1)**2 + 1) // 255 /*this defines/executes the function F.*/</lang> output   is the same as the 1st REXX version   (Brent's algorithm).

hash table algorithm

This REXX version is a lot faster   (than the sequential search algorithm)   if the   cycle length   and/or   start index   is large. <lang rexx>/*REXX pgm detects a cycle in an iterated function [F] using a hash table. */ x=3; @=x;  !.=.;  !.x=1 

                        do n=1+words(@); x=f(x); @=@ x /*add number to list. */
                        if !.x\==.  then leave         /*A repeat? Then leave*/
                        !.x=n                          /*N: numbers in @ list*/
                        end   /*n*/

say ' original list=' @ ... /*maybe show the list.*/ say 'numbers in list=' n /*display num of nums.*/ say ' cycle length =' n-!.x /* " " cycle. */ say ' start index =' !.x-1 ' ◄─── zero based' /* " " index. */ say 'cycle sequence =' subword(@, !.x, n-!.x) /* " " sequence*/ exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ f: return (arg(1)**2 + 1) // 255 /*this defines/executes the function F.*/</lang> output   is the same as the 1st REXX version   (Brent's algorithm).

robust hash table algorithm

This REXX version implements a   robust   hash table algorithm, which means that when the divisor
(in the   F   function)   is fairly large, the cycle length can be very large, making the creation of the
iterated function sequence.   This becomes problematic when the mechanism to append a number to
the sequence becomes time consuming because the sequence contains many hundreds of thousands
of numbers.

This REXX version allows the divisor for the   F   function to be specified, which can be chosen to stress
test the hash table algorithm.   A divisor which is   two raised to the 49th power   was chosen;   it
generates a cyclic sequence that contains over 1.5 million numbers. <lang rexx>/*REXX pgm detects a cycle in an iterated function [F] using robust hashing.*/ parse arg power . /*obtain optional arg. */ if power== | power="," then power=8 /*use the default power*/ numeric digits 500 /*handle the big number*/ divisor=2**power-1 /*compute the divisor. */ numeric digits max(9, length(divisor) * 2 + 1) /*allow for square + 1.*/ say ' power =' power /*display the power. */ say ' divisor =' "{2**"power'}-1 = ' divisor /* " " divisor. */ say x=3; @=x; @@=; m=100 /*M: max nums to show.*/ !.=.;  !.x=1; do n=1+words(@); x=f(x); @@=@@ x 

                     if n//2000==0  then do;    @=@ @@;  @@=;  end   /*rejoin*/
                     if !.x\==.     then leave        /*this number a repeat?*/
                     !.x=n
                     end   /*n*/                      /*N:  size of  @  list.*/

@=space(@ @@) /*append residual nums.*/ if n<m then say 'original list=' @ ... /*maybe show the list. */ 

            say 'cycle length =' n-!.x                /*display the cycle.   */
            say 'start index  =' !.x-1     "  ◄─── zero based"  /*show index.*/

if n<m then say 'the sequence =' subword(@,!.x,n-!.x) /*maybe show the seq. */ exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ f: return (arg(1)**2 + 1) // divisor</lang> output   when the input (power of two) used is:   49

Note that the listing of the original list and the cyclic sequence aren't displayed as they are much too long. 

       power = 49
     divisor = {2**49}-1 =  562949953421311

cycle length = 1500602
start index  = 988379   ◄─── zero based

variant of the hash table algorithm

There is more information in the "hash table"
and f has no "side effect". <lang rexx>/*REXX pgm detects a cycle in an iterated function [F] */ x=3; list=x; p.=0; p.x=1 Do q=2 By 1 

 x=f(x)                           /* the next value                 */
 list=list x                      /* append it to the list          */
 If p.x>0 Then Do                 /* x was in the list              */
   cLen=q-p.x                     /* cycle length                   */
   Leave
   End
 p.x=q                            /* note position of x in the list */
 End

Say 'original list=' list ... Say 'cycle length =' cLen /*display the cycle len */ Say 'start index =' p.x-1 " (zero based)" /* " " index. */ Say 'the sequence =' subword(list,p.x,cLen) /* " " sequence. */ Exit /*-------------------------------------------------------------------*/ f: Return (arg(1)**2+1)// 255; /*define the function F*/</lang>

Ruby

Works with: ruby version 2.0

<lang Ruby># Author: Paul Anton Chernoch

  1. Purpose:
  2. Find the cycle length and start position of a numerical seried using Brent's cycle algorithm.
  4. Given a recurrence relation X[n+1] = f(X[n]) where f() has
  5. a finite range, you will eventually repeat a value that you have seen before.
  6. Once this happens, all subsequent values will form a cycle that begins
  7. with the first repeated value. The period of that cycle may be of any length.
  9. Parameters:
  10. x0 ...... First integer value in the sequence
  11. block ... Block that takes a single integer as input
  12. and returns a single integer as output.
  13. This yields a sequence of numbers that eventually repeats.
  14. Returns:
  15. Two values: lambda and mu
  16. lambda .. length of cycle
  17. mu ...... zero-based index of start of cycle

def findCycle(x0) 

 power = lambda = 1
 tortoise = x0
 hare = yield(x0)
 
 # Find lambda, the cycle length
 while tortoise != hare
   if power == lambda
     tortoise = hare
     power *= 2
     lambda = 0
   end
   hare = yield(hare)
   lambda += 1
 end
 
 # Find mu, the zero-based index of the start of the cycle
 hare = x0
 lambda.times { hare = yield(hare) }
 
 tortoise, mu = x0, 0
 while tortoise != hare
   tortoise = yield(tortoise)
   hare = yield(hare)
   mu += 1
 end
 
 return lambda, mu

end

  1. A recurrence relation to use in testing

def f(x) (x * x + 1) % 255 end

  1. Display the first 41 numbers in the test series

puts (1..40).reduce([3]){|acc,_| acc << f(acc.last)}.join(",")

  1. Test the findCycle function

clength, cstart = findCycle(3) { |x| f(x) } puts "Cycle length = #{clength}\nStart index = #{cstart}"</lang>

Output:
3,10,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5,26,167,95,101,2,5
Cycle length = 6
Start index = 2

zkl

Algorithm from the Wikipedia <lang zkl>fcn cycleDetection(f,x0){ // f(int), x0 is the integer starting value of the sequence 

 # main phase: search successive powers of two
 power:=lam:=1;
 tortoise,hare:=x0,f(x0);  # f(x0) is the element/node next to x0.
 while(tortoise!=hare){
    if(power==lam){  # time to start a new power of two?

tortoise,lam=hare,0; power*=2; 

    }
    hare=f(hare);
    lam+=1;
 }
 # Find the position of the first repetition of length λ
 mu:=0; tortoise=hare=x0;
 do(lam){ hare=f(hare) } # The distance between the hare and tortoise is now λ

 # Next, the hare and tortoise move at same speed till they agree
 while(tortoise!=hare){
    tortoise,hare=f(tortoise),f(hare);
    mu+=1;
 } 
 return(lam,mu);

}</lang> <lang zkl>cycleDetection(fcn(x){ (x*x + 1)%255 }, 3).println(" == cycle length, start index");</lang>

Output:
L(6,2) == cycle length, start index
Retrieved from "https://rosettacode.org/wiki/Cycle_detection?oldid=264740"