Cycle detection

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

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.

Task

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

11l

Translation of: D

<lang 11l>F print_result(x0, f, len, start)

  print(‘Cycle length = ’len)
  print(‘Start index = ’start)
  V i = x0
  L 1..start
     i = f(i)
  V cycle = [0] * len
  L 0.<len
     cycle[L.index] = i
     i = f(i)
  print(‘Cycle: ’, end' ‘’)
  print(cycle)

F brent(f, x0)

  Int cycle_length
  V hare = x0
  V power = 1
  L
     V tortoise = hare
     L(i) 1..power
        hare = f(hare)
        I tortoise == hare
           cycle_length = i
           ^L.break
     power *= 2

  hare = x0
  L 1..cycle_length
     hare = f(hare)

  V cycle_start = 0
  V tortoise = x0
  L tortoise != hare
     tortoise = f(tortoise)
     hare = f(hare)
     cycle_start++

  print_result(x0, f, cycle_length, cycle_start)

brent(i -> (i * i + 1) % 255, 3)</lang>

Output:

Cycle length = 6
Start index = 2
Cycle: [101, 2, 5, 26, 167, 95]

8086 Assembly

<lang asm> cpu 8086 org 100h section .text mov ax,3 ; Print the first 20 values in the sequence mov bx,f xor cx,cx mov dx,20 call seq mov ax,3 ; Use Brent's algorithm to find the cycle mov bx,f call brent mov bp,ax ; BP = start of cycle mov di,dx ; DI = length of cycle mov dx,nl ; Print a newline call prstr mov dx,len ; Print "Length: " call prstr mov ax,di ; Print the length call prnum mov dx,ix ; Print "Index: " call prstr mov ax,bp ; Print the index call prnum mov dx,nl ; Print another newline call prstr mov ax,3 ; Print the cycle mov bx,f mov cx,bp mov dx,di jmp seq ;;; Brent's algorithm ;;; Input: AX = x0, BX = address of function ;;; Output: AX = mu, DX = lambda, CX, SI, BP destroyed ;;; The routine in BX must preserve BX, CX, DX, SI, BP brent: mov bp,ax ; BP = x0 mov cx,ax ; CX = tortoise xor dx,dx ; DX = lambda mov si,1 ; SI = power .lam: call bx ; Apply function (AX = hare) inc dx ; Lambda += 1 cmp ax,cx ; Done yet? je .mu cmp dx,si ; Time to start a new power of two? jne .lam mov cx,ax ; Tortoise = hare shl si,1 ; power *= 2 xor dx,dx ; lambda = 0 jmp .lam .mu: mov ax,bp ; Find position of first repetition mov cx,dx ; CX = lambda .apply: call bx loop .apply mov cx,bp ; CX = tortoise xor si,si ; SI = mu .lmu: cmp ax,cx ; Done yet? je .done call bx ; hare = f(hare) xchg ax,cx ; tortoise = f(tortoise) call bx xchg ax,cx inc si jmp .lmu .done: mov ax,si ret ;;; Function to use ;;; AX = (AX*AX+1) mod 255 f: mul al ; AX = AL*AL (we know anything mod 255 is <256) inc ax ; + 1 div byte [fdsor] ; This is faster than freeing up a byte register xchg al,ah ; Remainder into AL xor ah,ah ; Pad to 16 bits ret ;;; -- Utility routines -- ;;; Print decimal value of AL. Destroys AX, DX. prnum: mov dx,nbuf ; Number output buffer xchg bx,dx .dgt: aam ; Extract digit add al,'0' ; Make ASCII digit dec bx mov [bx],al ; Store in buffer xchg ah,al ; Continue with rest of digits test al,al ; As long as there are digits jnz .dgt xchg bx,dx ; Put buffer pointer back in DX prstr: mov ah,9 ; Print using MS-DOS call. int 21h ret ;;; Print DX values in sequence x, f(x), f(f(x))... starting at CX. ;;; AX = x0, BX = function. AX, CX, DX, SI destroyed. seq: test cx,cx ; CX = 0? jz .print .skip: call bx ; Otherwise skip CX numbers loop .skip .print: mov cx,dx ; Amount of numbers to print .loop: mov si,ax ; Keep a copy of AX (print routine trashes it) call prnum mov ax,si ; Restore x call bx ; Find the next value loop .loop ret section .data fdsor: db 255 ; This 255 is used for the modulus calculation db '...' ; Buffer for byte-sized numeric output nbuf: db ' $' ix: db 'Index: $' len: db 'Length: $' nl: db 13,10,'$'</lang>

Output:

3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95
Length: 6 Index: 2
101 2 5 26 167 95 

Ada

This implementation is split across three files. The first is the specification of a generic package: <lang ada> generic

  type Element_Type is private;

package Brent is

  type Brent_Function is access function (X : Element_Type) return Element_Type;
  procedure Brent(F : Brent_Function; X0 : Element_Type; Lambda : out Integer; Mu : out Integer);

end Brent; </lang> The second is the body of the generic package: <lang ada> package body Brent is

  procedure Brent (F : Brent_Function; X0 : Element_Type; Lambda : out Integer; Mu : out Integer) is
     Power : Integer := 1;
     Tortoise : Element_Type := X0;
     Hare : Element_Type := F(X0);
  begin
     Lambda := 1;
     Mu := 0;
     while Tortoise /= Hare loop
        if Power = Lambda then
           Tortoise := Hare;
           Power := Power * 2;
           Lambda := 0;
        end if;
        Hare := F(Hare);
        Lambda := Lambda + 1;
     end loop;
     Tortoise := X0;
     Hare := X0;
     for I in 0..(Lambda-1) loop
        Hare := F(Hare);
     end loop;
     while Hare /= Tortoise loop
        Tortoise := F(Tortoise);
        Hare := F(Hare);
        Mu := Mu + 1;
     end loop;
  end Brent;

end Brent; </lang> By using a generic package, this implementation can be made to work for any unary function, not just integers. As a demonstration two instances of the test function are created and two instances of the generic package. These are produced inside the main procedure: <lang ada> with Brent; with Text_Io; use Text_Io; procedure Main is

  package Integer_Brent is new Brent(Element_Type => Integer);
  use Integer_Brent;
  function F (X : Integer) return Integer is
    ((X * X + 1) mod 255);
  type Mod255 is mod 255;
  package Mod255_Brent is new Brent(Element_Type => Mod255);
  function F255 (X : Mod255) return Mod255 is
     (X * X + 1);
  lambda : Integer;
  Mu : Integer;
  X : Integer := 3;

begin

  for I in 1..41 loop
     Put(Integer'Image(X));
     if I < 41 then
        Put(",");
     end if;
     X := F(X);
  end loop;
  New_Line;
  Integer_Brent.Brent(F'Access, 3, Lambda, Mu);
  Put_Line("Cycle Length: " & Integer'Image(Lambda));
  Put_Line("Start Index : " & Integer'Image(Mu));

  Mod255_Brent.Brent(F255'Access, 3, Lambda, Mu);
  Put_Line("Cycle Length: " & Integer'Image(Lambda));
  Put_Line("Start Index : " & Integer'Image(Mu));

  Put("Cycle       : ");
  X := 3;
  for I in 0..(Mu + Lambda) loop
     if Mu <= I and I < (Lambda + Mu) then
        Put(Integer'Image(X));
     end if;
     X := F(X);
  end loop;

end Main; </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
Cycle Length:  6
Start Index :  2
Cycle       :  101 2 5 26 167 95

APL

Works with: Dyalog APL

<lang apl>brent←{

   f←⍺⍺
   lam←⊃{
       l p t h←⍵
       p=l: 1 (p×2) h (f h) ⋄ (l+1) p t (f h)
   }⍣{=/2↓⍺} ⊢ 1 1 ⍵ (f ⍵)
   mu←⊃{
       (⊃⍵+1),f¨1↓⍵
   }⍣{=/1↓⍺} ⊢ 0 ⍵ (f⍣lam⊢⍵)
   mu lam

}

task←{

   seq←{f←⍺⍺ ⋄ (⊃⍺)↓{⍵,f⊃⌽⍵}⍣(1-⍨+/⍺)⊢⍵}
   ⎕←0 20 ⍺⍺ seq ⍵                    ⍝ First 20 elements
   ⎕←(↑'Index' 'Length'),⍺⍺ brent ⍵   ⍝ Index and length of cycle
   ⎕←(⍺⍺ brent ⍺⍺ seq⊢)⍵              ⍝ Cycle

}

(255 | 1 + ⊢×⊢) task 3</lang>

Output:

3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95
Index  2
Length 6
101 2 5 26 167 95

BCPL

<lang BCPL>get "libhdr"

// Brent's algorithm. 'fn' is a function pointer, // 'lam' and 'mu' are output pointers. let brent(fn, x0, lam, mu) be $( let power, tort, hare = 1, x0, fn(x0)

   !lam := 1
   
   until tort = hare do
   $(  if power = !lam then
       $(  tort := hare
           power := power * 2
           !lam := 0
       $)
       hare := fn(hare)
       !lam := !lam + 1
   $)
   
   tort := x0
   hare := x0
   for i = 1 to !lam do
       hare := fn(hare)
   
   !mu := 0
   until tort = hare do
   $(  tort := fn(tort)
       hare := fn(hare)
       !mu := !mu + 1
   $)

$)

// Print fn^m(x0) to fn^n(x0) let printRange(fn, x0, m, n) be $( for i = 0 to n-1 do

   $(  if m<=i then writef("%N ", x0)
       x0 := fn(x0)
   $)
   wrch('*N')

$)

let start() be $( let lam, mu = 0, 0

   // the function to find a cycle in
   let f(x) = (x*x + 1) rem 255
   
   // print the first 20 values
   printRange(f, 3, 0, 20)
   
   // find the cycle
   brent(f, 3, @lam, @mu)
   writef("Length: %N*NIndex: %N*N", lam, mu)
   
   // print the cycle
   printRange(f, 3, mu, mu+lam)

$)</lang>

Output:

3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95
Length: 6
Index: 2
101 2 5 26 167 95

C

Translation of: Modula-2

<lang c>#include <stdio.h>

1. include <stdlib.h>

typedef int(*I2I)(int); typedef struct {

   int a, b;

} Pair;

Pair brent(I2I f, int x0) {

   int power = 1, lam = 1, tortoise = x0, hare, mu, i;
   Pair result;

   hare = (*f)(x0);
   while (tortoise != hare) {
       if (power == lam) {
           tortoise = hare;
           power = power * 2;
           lam = 0;
       }
       hare = (*f)(hare);
       lam++;
   }

   hare = x0;
   i = 0;
   while (i < lam) {
       hare = (*f)(hare);
       i++;
   }

   tortoise = x0;
   mu = 0;
   while (tortoise != hare) {
       tortoise = (*f)(tortoise);
       hare = (*f)(hare);
       mu++;
   }

   result.a = lam;
   result.b = mu;
   return result;

}

int lambda(int x) {

   return (x*x + 1) % 255;

}

int main() {

   int x0 = 3, x = 3, i;
   Pair result;

   printf("[3");
   for (i = 1; i <= 40; ++i) {
       x = lambda(x);
       printf(", %d", x);
   }
   printf("]\n");

   result = brent(lambda, x0);
   printf("Cycle length = %d\nStart index = %d\nCycle = [", result.a, result.b);

   x0 = 3;
   x = x0;
   for (i = 1; i <= result.b; ++i) {
       x = lambda(x);
   }
   for (i = 1; i <= result.a; ++i) {
       if (i > 1) {
           printf(", ");
       }

       printf("%d", x);
       x = lambda(x);
   }

   printf("]\n");

   return 0;

}</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
Cycle = [101, 2, 5, 26, 167, 95]

C#

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>

C++

<lang cpp>struct ListNode {

     int val;
     ListNode *next;
     ListNode(int x) : val(x), next(NULL) {}
};

ListNode* Solution::detectCycle(ListNode* A) {

   ListNode* slow = A;
   ListNode* fast = A;
   ListNode* cycleNode = 0;
   while (slow && fast && fast->next)
   {
       slow = slow->next;
       fast = fast->next->next;
       if (slow == fast)
       {
           cycleNode = slow;
           break;
       }
   }
   if (cycleNode == 0)
   {
       return 0;
   }
   std::set<ListNode*> setPerimeter;
   setPerimeter.insert(cycleNode);
   for (ListNode* pNode = cycleNode->next; pNode != cycleNode; pNode = pNode->next)
   setPerimeter.insert(pNode);
   for (ListNode* pNode = A; true; pNode = pNode->next)
   {
       std::set<ListNode*>::iterator iter = setPerimeter.find(pNode);
       if (iter != setPerimeter.end()) 
       {
           return pNode;
       }
   }

}</lang>

Cowgol

<lang cowgol>include "cowgol.coh";

typedef N is uint8; interface BrentFn(x: N): (r: N);

sub Brent(f: BrentFn, x0: N): (lam: N, mu: N) is

   var power: N := 1;
   var tortoise := x0;
   var hare := f(x0);
   while tortoise != hare loop
       if power == lam then
           tortoise := hare;
           power := power << 1;
           lam := 0;
       end if;
       hare := f(hare);
       lam := lam + 1;
   end loop;
   
   tortoise := x0;
   hare := x0;
   var i: N := 0;
   while i < lam loop
       i := i + 1;
       hare := f(hare);
   end loop;
   
   mu := 0;
   while tortoise != hare loop
       tortoise := f(tortoise);
       hare := f(hare);
       mu := mu + 1;
   end loop;

end sub;

sub PrintRange(f: BrentFn, x0: N, start: N, end_: N) is

   var i: N := 0;
   while i < end_ loop
       if i >= start then
           print_i32(x0 as uint32);
           print_char(' ');
       end if;
       i := i + 1;
       x0 := f(x0);
   end loop;
   print_nl();

end sub;

sub x2_plus1_mod255 implements BrentFn is

   r := ((x as uint16 * x as uint16 + 1) % 255) as uint8;

end sub;

PrintRange(x2_plus1_mod255, 3, 0, 20); var length: N; var start: N; (length, start) := Brent(x2_plus1_mod255, 3); print("Cycle length: "); print_i32(length as uint32); print_nl(); print("Cycle start: "); print_i32(start as uint32); print_nl(); PrintRange(x2_plus1_mod255, 3, start, length+start);</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: 2
101 2 5 26 167 95

D

Translation of: Java

<lang D>import std.range; import std.stdio;

void main() {

   brent(i => (i * i + 1) % 255, 3);

}

void brent(int function(int) 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(hare);
            if (tortoise == hare) {
               cycleLength = i;
               break FOUND;
           }
       }
   }

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

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

   printResult(x0, f, cycleLength, cycleStart);

}

void printResult(int x0, int function(int) f, int len, int start) {

   writeln("Cycle length: ", len);
   writefln("Cycle: %(%s %)", iterate(x0, f).drop(start).take(len));

}

auto iterate(int start, int function(int) f) {

   return only(start).chain(generate!(() => start=f(start)));

}</lang>

Output:

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

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

Stream.iterate(3, &f.(&1)) |> Enum.take(41) |> 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

Factor

This is a strict translation of the Python code from the Wikipedia article. Perhaps a more idiomatic version could be added in the future, although there is value in showing how Factor's lexical variables differ from most languages. A variable binding with ! at the end is mutable, and subsequent uses of that name followed by ! change the value of the variable to the value at the top of the data stack. <lang factor>USING: formatting kernel locals make math prettyprint ;

cyclical-function ( n -- m ) dup * 1 + 255 mod ;

brent ( x0 quot -- λ μ )

   1 dup :> ( power! λ! )
   x0 :> tortoise!
   x0 quot call :> hare!
   [ tortoise hare = not ] [
       power λ = [
           hare tortoise!
           power 2 * power!
           0 λ! 
       ] when
       hare quot call hare!
       λ 1 + λ!
   ] while

   0 :> μ!
   x0 dup tortoise! hare!
   λ [ hare quot call hare! ] times
   [ tortoise hare = not ] [
       tortoise quot call tortoise!
       hare quot call hare!
       μ 1 + μ!
   ] while
   λ μ ; inline

3 [ 20 [ dup , cyclical-function ] times ] { } make nip . 3 [ cyclical-function ] brent "Cycle length: %d\nCycle start: %d\n" printf</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: 2

FOCAL

<lang focal>01.10 S X0=3;T %3 01.20 S X=X0;F I=1,20;T X;D 2 01.30 D 3;T !" START",M,!,"LENGTH",L,! 01.40 S X=X0;F I=1,M;D 2 01.50 F I=1,L;T X;D 2 01.60 T ! 01.70 Q

02.01 C -- X = X*X+1 MOD 255 02.10 S X=X*X+1 02.20 S X=X-255*FITR(X/255)

03.01 C -- BRENT'S ALGORITHM ON ROUTINE 2 03.10 S X=X0;S P=1;S L=1;S T=X0;D 2 03.20 I (T-X)3.3,3.6,3.3 03.30 I (P-L)3.5,3.4,3.5 03.40 S T=X;S P=P*2;S L=0 03.50 D 2;S L=L+1;G 3.2 03.60 S T=X0;S X=X0;F I=1,L;D 2 03.70 S H=X;S M=0 03.80 I (T-H)3.9,3.99,3.9 03.90 S X=T;D 2;S T=X;S X=H;D 2;S H=X;S M=M+1;G 3.8 03.99 R</lang>

Output:

=   3=  10= 101=   2=   5=  26= 167=  95= 101=   2=   5=  26= 167=  95= 101=   2=   5=  26= 167=  95
 START=   2
LENGTH=   6
= 101=   2=   5=  26= 167=  95

Forth

Works with gforth. <lang Forth>

cycle-length { x0 'f -- lambda } \ Brent's algorithm stage 1

   1 1 x0 dup 'f execute
   begin 2dup <> while
       2over = if
           2swap nip 2* 0
           2swap nip dup
       then
       'f execute rot 1+ -rot
   repeat 2drop nip ;

iterations ( x f n -- x )

   >r swap r> 0 ?do over execute loop nip ;

cycle-start { x0 'f lambda -- mu } \ Brent's algorithm stage 2

   0 x0 dup 'f lambda iterations
   begin 2dup <> while
       swap 'f execute  swap 'f execute  rot 1+ -rot
   repeat 2drop ;

find-cycle ( x0 'f -- mu lambda ) \ Brent's algorithm

   2dup cycle-length dup >r cycle-start r> ;

\ --- usage ---

.cycle { start len x0 'f -- }

   x0
   start 1- 0 do  'f execute  loop
   len 0 do 'f execute dup .  loop
   drop ;

f(x) dup * 1+ 255 mod ;

3 ' f(x) find-cycle ." The cycle starts at offset " over . ." and has a length of " dup . cr ." The cycle is " 3 ' f(x) .cycle cr bye </lang>

Output:

The cycle starts at offset 2 and has a length of 6 
The cycle is 101 2 5 26 167 95 

FreeBASIC

Brent's algorithm

Translation of: Python

<lang freebasic>' version 11-01-2017 ' compile with: fbc -s console

' define the function f(x)=(x*x +1) mod 255 Function f(x As Integer) As Integer

   Return (x * x +1) Mod 255

End Function

Sub brent(x0 As Integer, ByRef lam As Integer, ByRef mu As Integer)

   Dim As Integer i, power = 1
   lam = 1

   ' main phase: search successive powers of two
   Dim As Integer tortoise = f(x0)   ' f(x0) is the element/node next to x0.
   Dim As Integer hare = f(f(x0))

   While tortoise <> hare
       If power = lam Then
           tortoise = hare
           power *= 2
           lam = 0
       End If
       hare = f(hare)
       lam += 1
   Wend

   ' Find the position of the first repetition of length ?
   mu = 0
   tortoise = x0
   hare = x0
   For i = 0 To lam -1
       ' range(lam) produces a list with the values 0, 1, ... , lam-1
       hare = f(hare)
   Next
   ' The distance between the hare and tortoise is now ?.

   ' Next, the hare and tortoise move at same speed until they agree
   While tortoise <> hare
       tortoise = f(tortoise)
       hare = f(hare)
       mu += 1
   Wend

End Sub

' ------=< MAIN >=------

Dim As Integer i, j, lam, mu, x0 = 3

brent(x0, lam, mu) Print " Brent's algorithm" Print " Cycle starts at position: "; mu; " (starting position = 0)" Print " The length of the Cycle = "; lam Print

j = f(x0) Print " Cycle: "; For i = 1 To lam + mu -1

   If i >= mu Then
       Print j;
       If i <> (lam + mu -1) Then Print ", "; Else Print ""
   End If
   j = f(j)

Next Print

' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>

Output:

 Brent's algorithm
 Cycle starts at position:  2 (starting position = 0)
 The length of the Cycle =  6

 Cycle:  101,  2,  5,  26,  167,  95

Tortoise and hare. Floyd's algorithm

Translation of: Wikipedia

<lang freebasic>' version 11-01-2017 ' compile with: fbc -s console

' define the function f(x)=(x*x +1) mod 255 Function f(x As Integer) As Integer

   Return (x * x +1) Mod 255

End Function

Sub floyd(x0 As Integer, ByRef lam As Integer, ByRef mu As Integer)

   ' Main phase of algorithm: finding a repetition x_i = x_2i.
   ' The hare moves twice as quickly as the tortoise and
   ' the distance between them increases by 1 at each step.
   ' Eventually they will both be inside the cycle and then,
   ' at some point, the distance between them will be
   ' divisible by the period ?.
   Dim As Integer tortoise = f(x0)   ' f(x0) is the element/node next to x0.
   Dim As Integer hare = f(f(x0))

   While tortoise <> hare
       tortoise = f(tortoise)
       hare = f(f(hare))
   Wend

   ' At this point the tortoise position, ?, which is also equal
   ' to the distance between hare and tortoise, is divisible by
   ' the period ?. So hare moving in circle one step at a time,
   ' and tortoise (reset to x0) moving towards the circle, will
   ' intersect at the beginning of the circle. Because the
   ' distance between them is constant at 2?, a multiple of ?,
   ' they will agree as soon as the tortoise reaches index µ.

   ' Find the position µ of first repetition.
   mu = 0
   tortoise = x0
   While tortoise <> hare
       tortoise = f(tortoise)
       hare = f(hare)         ' Hare and tortoise move at same speed
       mu += 1
   Wend

   ' Find the length of the shortest cycle starting from x_µ
   ' The hare moves one step at a time while tortoise is still.
   ' lam is incremented until ? is found.
   lam = 1
   hare = f(tortoise)
   While tortoise <> hare
       hare = f(hare)
       lam += 1
   Wend

End Sub

' ------=< MAIN >=------

Dim As Integer i, j, lam, mu, x0 = 3

floyd(x0, lam, mu) Print " Tortoise and hare. Floyd's algorithm" Print " Cycle starts at position: "; mu; " (starting position = 0)" Print " The length of the Cycle = "; lam Print

j = f(x0) Print " Cycle: "; For i = 1 To lam + mu -1

   If i >= mu Then
       Print j;
       If i <> (lam + mu -1) Then Print ", "; Else Print ""
   End If
   j = f(j)

Next Print

' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>

Go

Translation of: Python

Translation of: Wikipedia

Run it on the go playground, or on go play space. <lang Go> package main

import "fmt"

func brent(f func(i int) int, x0 int) (int, int) { var λ, µ, power, tortoise, hare int

// Main phase: search successive powers of two. power = 1 λ = 1 tortoise = x0 hare = f(x0) // f(x0) is the element/node next to x0.

for tortoise != hare { if power == λ { // Time to start a new power of two. tortoise = hare power *= 2 λ = 0 } hare = f(hare) λ++ }

// Find the position of the first repetition of length λ. µ = 0 tortoise, hare = x0, x0 for i := 0; i < λ; i++ { // produces a list with the values 0,1,...,λ-1. hare = f(hare) // The distance between hare and tortoise is now λ. }

// The tortoise and the hare move at the same speed until they agree. for tortoise != hare { tortoise = f(tortoise) hare = f(hare) µ++ }

return λ, µ }

func f(i int) int { return (i*i + 1) % 255 }

func main() { x0 := 3 λ, µ := brent(f, x0) fmt.Println("Cycle length:", λ) fmt.Println("Cycle start index:", µ) a := []int{} for i := 0; i <= λ+1; i++ { a = append(a, x0) x0 = f(x0) } fmt.Println("Cycle:", a[µ:µ+λ]) }</lang>

Output:

Cycle length: 6 
Cycle start index: 2 
Cycle: [101 2 5 26 167 95]

Groovy

Translation of: Java

<lang groovy>import java.util.function.IntUnaryOperator

class CycleDetection {

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

   static void brent(IntUnaryOperator f, int x0) {
       int cycleLength = -1
       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) {
       printf("Cycle length: %d%nCycle: ", len)

       int n = x0
       for (int i = 0; i < start + len; i++) {
           n = f.applyAsInt(n)
           if (i >= start) {
               printf("%s ", n)
           }
       }
       println()
   }

}</lang>

Output:

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

Haskell

Most of solutions given in other languages are not total. For function which does not have any cycles under iteration (i.e. f(x) = 1 + x) they will never terminate.

Haskellers, being able to handle conceptually infinite structures, traditionally consider totality as an important issue. The following solution is total for total inputs (modulo totality of iterated function) and allows nontermination only if input is explicitly infinite.

<lang haskell>import Data.List (findIndex)

findCycle :: Eq a => [a] -> Maybe ([a], Int, Int) findCycle lst =

 do l <- findCycleLength lst
    mu <- findIndex (uncurry (==)) $ zip lst (drop l lst)
    let c = take l $ drop mu lst
    return (c, l, mu)

findCycleLength :: Eq a => [a] -> Maybe Int findCycleLength [] = Nothing findCycleLength (x:xs) =

 let loop _ _ _ [] = Nothing
     loop pow lam x (y:ys)
       | x == y     = Just lam
       | pow == lam = loop (2*pow) 1 y ys
       | otherwise  = loop pow (1+lam) x ys
 in loop 1 1 x xs</lang>

Examples

λ> findCycle (cycle [1,2,3])
Just ([1,2,3],3,0)

λ> findCycle ([1..100] ++ cycle [1..12])
Just ([1,2,3,4,5,6,7,8,9,10,11,12],12,100)

λ> findCycle [1..1000]
Nothing

λ> findCycle (iterate (\x -> (x^2 + 1) `mod` 255) 3)
Just ([101,2,5,26,167,95],6,2)

λ> findCycle (take 100 $ iterate (\x -> x+1) 3)
Nothing

λ> findCycle (take 100000 $ iterate (\x -> x+1) 3)
Nothing

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

Julia

Works with: Julia version 0.6

Following the Wikipedia article:

<lang julia>using IterTools

function floyd(f, x0)

   local tort = f(x0)
   local hare = f(tort)
   while tort != hare
       tort = f(tort)
       hare = f(f(hare))
   end

   local μ = 0
   tort = x0
   while tort != hare
       tort = f(tort)
       hare = f(hare)
       μ += 1
   end

   λ = 1
   hare = f(tort)
   while tort != hare
       hare = f(hare)
       λ += 1
   end

   return λ, μ

end

f(x) = (x * x + 1) % 255

λ, μ = floyd(f, 3) cycle = iterate(f, 3) |>

   x -> Iterators.drop(x, μ) |>
   x -> Iterators.take(x, λ) |>
   collect

println("Cycle length: ", λ, "\nCycle start index: ", μ, "\nCycle: ", join(cycle, ", "))</lang>

Output:

Cycle length: 6
Cycle start index: 2
Cycle: 101, 2, 5, 26, 167, 95

Kotlin

<lang scala>// version 1.1.2

typealias IntToInt = (Int) -> Int

fun brent(f: IntToInt, x0: Int): Pair<Int, Int> {

   // main phase: search successive powers of two
   var power = 1
   var lam = 1
   var tortoise = x0
   var 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++
   }

   // Find the position of the first repetition of length 'lam'
   var mu = 0
   tortoise = x0
   hare = x0
   for (i in 0 until lam) 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++
   }
   return Pair(lam, mu)

}

fun main(args: Array<String>) {

   val f = { x: Int -> (x * x + 1) % 255 }
   // generate first 41 terms of the sequence starting from 3
   val x0 = 3
   var x = x0
   val seq = List(41) { if (it > 0) x = f(x) ; x }
   println(seq)
   val (lam, mu) = brent(f, x0)
   val cycle = seq.slice(mu until mu + lam)
   println("Cycle length = $lam")
   println("Start index  = $mu")
   println("Cycle        = $cycle")

}</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
Cycle        = [101, 2, 5, 26, 167, 95]

Lua

Fairly direct translation of the Wikipedia code, except that the sequence is stored in a table and passed back as a third return value. <lang Lua>function brent (f, x0)

   local tortoise, hare, mu = x0, f(x0), 0
   local cycleTab, power, lam = {tortoise, hare}, 1, 1
   while tortoise ~= hare do
       if power == lam then
           tortoise = hare
           power = power * 2
           lam = 0
       end
       hare = f(hare)
       table.insert(cycleTab, hare)
       lam = lam + 1
   end
   tortoise, hare = x0, x0
   for i = 1, lam do hare = f(hare) end
   while tortoise ~= hare do
       tortoise = f(tortoise)
       hare = f(hare)
       mu = mu + 1
   end
   return lam, mu, cycleTab

end

local f = function (x) return (x * x + 1) % 255 end local x0 = 3 local cycleLength, startIndex, sequence = brent(f, x0) print("Sequence:", table.concat(sequence, " ")) print("Cycle length:", cycleLength) print("Start Index:", startIndex)</lang>

Output:

Sequence:       3 10 101 2 5 26 167 95 101 2 5 26 167 95
Cycle length:   6
Start Index:    2

Modula-2

Translation of: Kotlin

<lang modula2>MODULE CycleDetection; FROM FormatString IMPORT FormatString; FROM Terminal IMPORT WriteString,WriteLn,ReadChar;

TYPE IntToInt = PROCEDURE(INTEGER) : INTEGER; TYPE Pair =

   RECORD
       a,b : INTEGER;
   END;

PROCEDURE Brent(f : IntToInt; x0 : INTEGER) : Pair; VAR power,lam,tortoise,hare,mu,i : INTEGER; BEGIN

   (* main phase: search successive powers of two *)
   power := 1;
   lam := 1;
   tortoise := x0;
   hare := f(x0);  (* f(x0) is the element/node next to x0. *)
   WHILE tortoise # hare DO
       (* time to start a new power of two? *)
       IF power = lam THEN
           tortoise := hare;
           power := power * 2;
           lam := 0
       END;
       hare := f(hare);
       INC(lam)
   END;

   (* Find the position of the first repetition of length 'lam' *)
   mu := 0;
   tortoise := x0;
   hare := x0;
   i := 0;
   WHILE i < lam DO
       hare := f(hare);
       INC(i)
   END;

   (* 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 DO
       tortoise := f(tortoise);
       hare := f(hare);
       INC(mu)
   END;

   RETURN Pair{lam,mu}

END Brent;

PROCEDURE Lambda(x : INTEGER) : INTEGER; BEGIN

   RETURN (x * x + 1) MOD 255

END Lambda;

VAR

   buf : ARRAY[0..63] OF CHAR;
   x0,x,i : INTEGER;
   result : Pair;

BEGIN

   x0 := 3;
   x := x0;

   WriteString("[3");
   FOR i:=1 TO 40 DO
       x := Lambda(x);
       FormatString(", %i", buf, x);
       WriteString(buf)
   END;
   WriteString("]");
   WriteLn;

   result := Brent(Lambda, x0);
   FormatString("Cycle length = %i\nStart index  = %i\nCycle        = [", buf, result.a, result.b);
   WriteString(buf);

   x0 := 3;
   x := x0;
   FOR i:=1 TO result.b DO
       x := Lambda(x)
   END;
   FOR i:=1 TO result.a DO
       IF i > 1 THEN
           WriteString(", ")
       END;

       FormatString("%i", buf, x);
       WriteString(buf);
       x := Lambda(x)
   END;

   WriteString("]");
   WriteLn;

   ReadChar

END CycleDetection.</lang>

Nim

Translation of Wikipedia Python program. <lang Nim>import strutils, sugar

func brent(f: int -> int; x0: int): (int, int) =

 # Main phase: search successive powers of two.
 var
   power, λ = 1
   tortoise = x0
   hare = f(x0)

 while tortoise != hare:
   if power == λ:
     # Time to start a new power of two.
     tortoise = hare
     power *= 2
     λ = 0
   hare = f(hare)
   inc λ

 # Find the position of the first repetition of length λ.
 tortoise = x0
 hare = x0
 for i in 0..<λ:
   hare = f(hare)
 # The distance between the hare and tortoise is now λ.

 # Next, the hare and tortoise move at same speed until they agree.
 var μ = 0
 while tortoise != hare:
   tortoise = f(tortoise)
   hare = f(hare)
   inc μ

 result = (λ, μ)

when isMainModule:

 func f(x: int): int = (x * x + 1) mod 255

 let x0 = 3
 let (λ, μ) = brent(f, x0)
 echo "Cycle length: ", λ
 echo "Cycle start index: ", μ
 var cycle: seq[int]
 var x = x0
 for i in 0..<λ+μ:
   if i >= μ: cycle.add x
   x = f(x)
 echo "Cycle: ", cycle.join(" ")</lang>

Output:

Cycle length: 6
Cycle start index: 2
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

Translation of: Raku

<lang perl>use utf8;

sub cyclical_function { ($_[0] * $_[0] + 1) % 255 }

sub brent {

   my($f, $x0) = @_;
   my $power = 1;
   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 0..$λ-1;

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

}

my ( $l, $s ) = brent( \&cyclical_function, 3 );

sub show_range {

   my($start,$stop) = @_;
   my $result;
   my $x = 3;
   for my $n (0..$stop) {
       $result .= "$x " if $n >= $start;
       $x = cyclical_function($x);
   }
   $result;

}

print show_range(0,19) . "\n"; print "Cycle length $l\n"; print "Cycle start index $s\n"; print show_range($s,$s+$l-1) . "\n";</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

Phix

Translation of the Wikipedia code, but using the more descriptive len and pos, instead of lambda and mu, and adding a limit.

function f(integer x)
    return mod(x*x+1,255)
end function
 
function brent(integer x0)
integer pow2 = 1, len = 1, pos = 1
integer tortoise = x0,
        hare = f(x0)
sequence s = {tortoise,hare} -- (kept for output only)
 
    -- main phase: search successive powers of two
    while tortoise!=hare do
        if pow2 = len then
            tortoise = hare
            pow2 *= 2
            if pow2>16 then
                return {s,0,0}
            end if
            len = 0
        end if
        hare = f(hare)
        s &= hare
        len += 1
    end while
 
    -- Find the position of the first repetition of length len
    tortoise = x0
    hare = x0
    for i=1 to len do
        hare = f(hare)
    end for
    -- The distance between the hare and tortoise is now len.
 
    -- Next, the hare and tortoise move at same speed until they agree
    while tortoise<>hare do
        tortoise = f(tortoise)
        hare = f(hare)
        pos += 1
    end while
 
    return {s,len,pos}
end function
 
sequence s
integer len, pos, x0 = 3
{s,len,pos} = brent(x0)
printf(1," Brent's algorithm\n")
?s
printf(1," Cycle starts at position: %d (1-based)\n",{pos})
printf(1," The length of the Cycle = %d\n",{len})
?s[pos..pos+len-1]

Output:

 Brent's algorithm
{3,10,101,2,5,26,167,95,101,2,5,26,167,95}
 Cycle starts at position: 3 (1-based)
 The length of the Cycle = 6
{101,2,5,26,167,95}

PL/M

<lang plm>100H: /* CP/M CALLS */ BDOS: PROCEDURE (FN, ARG); DECLARE FN BYTE, ARG ADDRESS; GO TO 5; END BDOS; EXIT: PROCEDURE; CALL BDOS(0,0); END EXIT; PRINT: PROCEDURE (S); DECLARE S ADDRESS; CALL BDOS(9,S); END PRINT;

/* THIS IS A HACK TO CALL A FUNCTION GIVEN ITS ADDRESS */ APPLY: PROCEDURE (FN, ARG) ADDRESS;

   DECLARE (FN, ARG) ADDRESS;
   INNER: PROCEDURE (X) ADDRESS; /* THIS SETS UP THE ARGUMENTS IN MEMORY */
       DECLARE X ADDRESS;
       GO TO FN;                 /* THEN WE JUMP TO THE ADDRESS GIVEN */
   END INNER;
   RETURN INNER(ARG);

END APPLY;

/* THIS IS ANOTHER HACK - THE SINGLE 0 (WHICH IS AN 8080 NOP) WILL BE

  STORED AS A CONSTANT RIGHT BEFORE THE FUNCTION.
  PL/M-80 DOES NOT ALLOW THE PROGRAMMER TO DIRECTLY GET THE ADDRESS
  OF A FUNCTION, BUT THE ADDRESS OF THIS CONSTANT WILL BE RIGHT IN
  FRONT OF IT AND CAN BE USED INSTEAD. */

DECLARE F$ADDR DATA (0); /* F(X) FROM THE TASK */ F: PROCEDURE (X) ADDRESS;

   DECLARE X ADDRESS;
   RETURN (X*X + 1) MOD 255;

END F;

/* PRINT A NUMBER */ PRINT$NUMBER: PROCEDURE (N);

   DECLARE S (7) BYTE INITIAL ('..... $');
   DECLARE (N, P) ADDRESS, C BASED P BYTE;
   P = .S(5);

DIGIT:

   P = P - 1;
   C = N MOD 10 + '0';
   N = N / 10;
   IF N > 0 THEN GO TO DIGIT;
   CALL PRINT(P);

END PRINT$NUMBER;

/* PRINT F^M(X0) TO F^N(X0) */ PRINT$RANGE: PROCEDURE (FN, X0, M, N);

   DECLARE (FN, X0, M, N, I) ADDRESS;
   DO I=0 TO N-1;
       IF I>=M THEN CALL PRINT$NUMBER(X0);
       X0 = APPLY(FN,X0);
   END;
   CALL PRINT(.(13,10,'$'));

END PRINT$RANGE;

/* BRENT'S ALGORITHM */ BRENT: PROCEDURE (FN, X0, MU$A, LAM$A);

   DECLARE (FN, X0, MU$A, LAM$A) ADDRESS;
   DECLARE (MU BASED MU$A, LAM BASED LAM$A) ADDRESS;
   DECLARE (TORT, HARE, POW, I) ADDRESS;

   POW, LAM = 1;
   TORT = X0;
   HARE = APPLY(FN,X0);
   DO WHILE TORT <> HARE;
       IF POW = LAM THEN DO;
           TORT = HARE;
           POW = SHL(POW, 1);
           LAM = 0;
       END;
       HARE = APPLY(FN,HARE);
       LAM = LAM + 1;
   END;
   
   TORT, HARE = X0;
   DO I=1 TO LAM;
       HARE = APPLY(FN,HARE);
   END;

   MU = 0;
   DO WHILE TORT <> HARE;
       TORT = APPLY(FN,TORT);
       HARE = APPLY(FN,HARE);
       MU = MU + 1;
   END;

END BRENT;

DECLARE (MU, LAM) ADDRESS;

/* PRINT THE FIRST 20 VALUES IN THE SEQUENCE */ CALL PRINT$RANGE(.F$ADDR, 3, 0, 20); /* FIND THE CYCLE */ CALL BRENT(.F$ADDR, 3, .MU, .LAM); CALL PRINT(.'LENGTH: $'); CALL PRINT$NUMBER(LAM); CALL PRINT(.'START: $'); CALL PRINT$NUMBER(MU); CALL PRINT(.(13,10,'$')); /* PRINT THE CYCLE */ CALL PRINT$RANGE(.F$ADDR, 3, MU, MU+LAM); CALL EXIT; EOF</lang>

Output:

3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95 
LENGTH: 6 START: 2 
101 2 5 26 167 95 

Python

Procedural

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]

Functional

In functional terms, the problem lends itself at first sight (see the Haskell version), to a reasonably compact recursive definition of a cycleLength function:

Translation of: Haskell

Works with: Python version 3.7

<lang python>Cycle detection by recursion.

from itertools import (chain, cycle, islice) from operator import (eq)

1. cycleFound :: Eq a => [a] -> Maybe ([a], Int, Int)

def cycleFound(xs):

   Just the first cycle found, with its length
      and start index, or Nothing if no cycle seen.
   
   return bind(cycleLength(xs))(
       lambda n: bind(
           findIndex(uncurry(eq))(zip(xs, xs[n:]))
       )(lambda iStart: Just(
           (xs[iStart:iStart + n], n, iStart)
       ))
   )

1. cycleLength :: Eq a => [a] -> Maybe Int

def cycleLength(xs):

   Just the length of the first cycle found,
      or Nothing if no cycle seen.
   
   def go(pwr, lng, x, ys):
       if ys:
           y, *yt = ys
           return Just(lng) if x == y else (
               go(2 * pwr, 1, y, yt) if (
                   lng == pwr
               ) else go(pwr, 1 + lng, x, yt)
           )
       else:
           return Nothing()

   return go(1, 1, xs[0], xs[1:]) if xs else Nothing()

1. TEST ----------------------------------------------------
2. main :: IO ()

def main():

   Reports of any cycle detection.

   print(
       fTable(
           'First cycle detected, if any:\n'
       )(fst)(maybe('No cycle found')(
           showCycle
       ))(
           compose(cycleFound)(snd)
       )([
           (
               'cycle([1, 2, 3])',
               take(1000)(cycle([1, 2, 3]))
           ), (
               '[0..100] + cycle([1..8])',
               take(1000)(
                   chain(
                       enumFromTo(0)(100),
                       cycle(enumFromTo(1)(8))
                   )
               )
           ), (
               '[1..500]',
               enumFromTo(1)(500)
           ), (
               'f(x) = (x*x + 1) modulo 255',
               take(1000)(iterate(
                   lambda x: (1 + (x * x)) % 255
               )(3))
           )
       ])
   )

1. DISPLAY -------------------------------------------------

1. showList :: [a] -> String

def showList(xs):

   'Compact stringification of a list,
       (no spaces after commas).
   
   return .join(repr(xs).split())

1. showCycle :: ([a], Int, Int) -> String

def showCycle(cli):

   Stringification of cycleFound tuple.
   c, lng, iStart = cli
   return showList(c) + ' (from:' + str(iStart) + (
       ', length:' + str(lng) + ')'
   )

1. GENERIC -------------------------------------------------

1. Just :: a -> Maybe a

def Just(x):

   Constructor for an inhabited Maybe (option type) value.
   return {'type': 'Maybe', 'Nothing': False, 'Just': x}

1. Nothing :: Maybe a

def Nothing():

   Constructor for an empty Maybe (option type) value.
   return {'type': 'Maybe', 'Nothing': True}

1. bind (>>=) :: Maybe a -> (a -> Maybe b) -> Maybe b

def bind(m):

   bindMay provides the mechanism for composing a
      sequence of (a -> Maybe b) functions.
      If m is Nothing, it is passed straight through.
      If m is Just(x), the result is an application
      of the (a -> Maybe b) function (mf) to x.
   return lambda mf: (
       m if m.get('Nothing') else mf(m.get('Just'))
   )

1. compose (<<<) :: (b -> c) -> (a -> b) -> a -> c

def compose(g):

   Right to left function composition.
   return lambda f: lambda x: g(f(x))

1. enumFromTo :: (Int, Int) -> [Int]

def enumFromTo(m):

   Integer enumeration from m to n.
   return lambda n: list(range(m, 1 + n))

1. findIndex :: (a -> Bool) -> [a] -> Maybe Int

def findIndex(p):

   Just the first index at which an
      element in xs matches p,
      or Nothing if no elements match.
   def go(xs):
       try:
           return Just(next(
               i for i, x in enumerate(xs) if p(x)
           ))
       except StopIteration:
           return Nothing()
   return lambda xs: go(xs)

1. fst :: (a, b) -> a

def fst(tpl):

   First member of a pair.
   return tpl[0]

1. fTable :: String -> (a -> String) ->
2. (b -> String) ->
3. (a -> b) -> [a] -> String

def fTable(s):

   Heading -> x display function -> fx display function ->
         f -> value list -> tabular string.
   def go(xShow, fxShow, f, xs):
       w = max(map(compose(len)(xShow), xs))
       return s + '\n' + '\n'.join([
           xShow(x).rjust(w, ' ') + ' -> ' + fxShow(f(x)) for x in xs
       ])
   return lambda xShow: lambda fxShow: (
       lambda f: lambda xs: go(
           xShow, fxShow, f, xs
       )
   )

1. iterate :: (a -> a) -> a -> Gen [a]

def iterate(f):

   An infinite list of repeated
      applications of f to x.
   def go(x):
       v = x
       while True:
           yield v
           v = f(v)
   return lambda x: go(x)

1. maybe :: b -> (a -> b) -> Maybe a -> b

def maybe(v):

   Either the default value v, if m is Nothing,
      or the application of f to x,
      where m is Just(x).
   return lambda f: lambda m: v if m.get('Nothing') else (
       f(m.get('Just'))
   )

1. snd :: (a, b) -> b

def snd(tpl):

   Second member of a pair.
   return tpl[1]

1. take :: Int -> [a] -> [a]
2. take :: Int -> String -> String

def take(n):

   The prefix of xs of length n,
      or xs itself if n > length xs.
   return lambda xs: (
       xs[0:n]
       if isinstance(xs, list)
       else list(islice(xs, n))
   )

1. uncurry :: (a -> b -> c) -> ((a, b) -> c)

def uncurry(f):

   A function over a tuple
      derived from a default or
      curried function.
   return lambda xy: f(xy[0], xy[1])

1. concat :: a -> [a]
2. concat :: [String] -> String

def concat(xxs):

   The concatenation of all the elements in a list.
   xs = list(chain.from_iterable(xxs))
   unit =  if isinstance(xs, str) else []
   return unit if not xs else (
       .join(xs) if isinstance(xs[0], str) else xs
   )

1. MAIN ---

if __name__ == '__main__':

   main()</lang>

Output:

First cycle detected, if any:

           cycle([1, 2, 3]) -> [1,2,3] (from:0, length:3)
   [0..100] + cycle([1..8]) -> [1,2,3,4,5,6,7,8] (from:101, length:8)
                   [1..500] -> No cycle found
f(x) = (x*x + 1) modulo 255 -> [101,2,5,26,167,95] (from:2, length:6)

But recursion scales poorly in Python, and the version above, while good for lists of a few hundred elements, will need reworking for longer lists and better use of space.

If we start by refactoring the recursion into the form of a higher order (but still recursive) until p f x function, we can then reimplement the internals of until itself to avoid recursion, without losing the benefits of compositional structure:

Recursive until: <lang python># until :: (a -> Bool) -> (a -> a) -> a -> a def until(p):

   The result of repeatedly applying f until p holds.
      The initial seed value is x.
   def go(f, x):
       return x if p(x) else go(f, f(x))
   return lambda f: lambda x: go(f, x)</lang>

cycleLength refactored in terms of until: <lang python># cycleLength :: Eq a => [a] -> Maybe Int def cycleLength(xs):

   Just the length of the first cycle found,
      or Nothing if no cycle seen.

   # f :: (Int, Int, Int, [Int]) -> (Int, Int, Int, [Int])
   def f(tpl):
       pwr, lng, x, ys = tpl
       y, *yt = ys
       return (2 * pwr, 1, y, yt) if (
           lng == pwr
       ) else (pwr, 1 + lng, x, yt)

   # p :: (Int, Int, Int, [Int]) -> Bool
   def p(tpl):
       _, _, x, ys = tpl
       return (not ys) or x == ys[0]

   if xs:
       _, lng, x, ys = until(p)(f)(
           (1, 1, xs[0], xs[1:])
       )
       return (
           Just(lng) if (x == ys[0]) else Nothing()
       ) if ys else Nothing()
   else:
       return Nothing()</lang>

Iterative reimplementation of until: <lang python># until_ :: (a -> Bool) -> (a -> a) -> a -> a def until_(p):

   The result of repeatedly applying f until p holds.
      The initial seed value is x.
   def go(f, x):
       v = x
       while not p(v):
           v = f(v)
       return v
   return lambda f: lambda x: go(f, x)</lang>

and now it all works again, with the structure conserved but recursion removed. The Python no longer falls out of the tree at the sight of an ouroboros, and we can happily search for cycles in lists of several thousand items:

Works with: Python version 3.7

<lang python>Cycle detection without recursion.

from itertools import (chain, cycle, islice) from operator import (eq)

1. cycleFound :: Eq a => [a] -> Maybe ([a], Int, Int)

def cycleFound(xs):

   Just the first cycle found, with its length
      and start index, or Nothing if no cycle seen.
   
   return bind(cycleLength(xs))(
       lambda n: bind(
           findIndex(uncurry(eq))(zip(xs, xs[n:]))
       )(lambda iStart: Just(
           (xs[iStart:iStart + n], n, iStart)
       ))
   )

1. cycleLength :: Eq a => [a] -> Maybe Int

def cycleLength(xs):

   Just the length of the first cycle found,
      or Nothing if no cycle seen.

   # f :: (Int, Int, Int, [Int]) -> (Int, Int, Int, [Int])
   def f(tpl):
       pwr, lng, x, ys = tpl
       y, *yt = ys
       return (2 * pwr, 1, y, yt) if (
           lng == pwr
       ) else (pwr, 1 + lng, x, yt)

   # p :: (Int, Int, Int, [Int]) -> Bool
   def p(tpl):
       _, _, x, ys = tpl
       return (not ys) or x == ys[0]

   if xs:
       _, lng, x, ys = until(p)(f)(
           (1, 1, xs[0], xs[1:])
       )
       return (
           Just(lng) if (x == ys[0]) else Nothing()
       ) if ys else Nothing()
   else:
       return Nothing()

1. TEST ----------------------------------------------------
2. main :: IO ()

def main():

   Reports of any cycle detection.

   print(
       fTable(
           'First cycle detected, if any:\n'
       )(fst)(maybe('No cycle found')(
           showCycle
       ))(
           compose(cycleFound)(snd)
       )([
           (
               'cycle([1, 2, 3])',
               take(10000)(cycle([1, 2, 3]))
           ), (
               '[0..10000] + cycle([1..8])',
               take(20000)(
                   chain(
                       enumFromTo(0)(10000),
                       cycle(enumFromTo(1)(8))
                   )
               )
           ), (
               '[1..10000]',
               enumFromTo(1)(10000)
           ), (
               'f(x) = (x*x + 1) modulo 255',
               take(10000)(iterate(
                   lambda x: (1 + (x * x)) % 255
               )(3))
           )
       ])
   )

1. DISPLAY -------------------------------------------------

1. showList :: [a] -> String

def showList(xs):

   'Compact stringification of a list,
       (no spaces after commas).
   
   return .join(repr(xs).split())

1. showCycle :: ([a], Int, Int) -> String

def showCycle(cli):

   Stringification of cycleFound tuple.
   c, lng, iStart = cli
   return showList(c) + ' (from:' + str(iStart) + (
       ', length:' + str(lng) + ')'
   )

1. GENERIC -------------------------------------------------

1. Just :: a -> Maybe a

def Just(x):

   Constructor for an inhabited Maybe (option type) value.
   return {'type': 'Maybe', 'Nothing': False, 'Just': x}

1. Nothing :: Maybe a

def Nothing():

   Constructor for an empty Maybe (option type) value.
   return {'type': 'Maybe', 'Nothing': True}

1. bind (>>=) :: Maybe a -> (a -> Maybe b) -> Maybe b

def bind(m):

   bindMay provides the mechanism for composing a
      sequence of (a -> Maybe b) functions.
      If m is Nothing, it is passed straight through.
      If m is Just(x), the result is an application
      of the (a -> Maybe b) function (mf) to x.
   return lambda mf: (
       m if m.get('Nothing') else mf(m.get('Just'))
   )

1. compose (<<<) :: (b -> c) -> (a -> b) -> a -> c

def compose(g):

   Right to left function composition.
   return lambda f: lambda x: g(f(x))

1. concat :: a -> [a]
2. concat :: [String] -> String

def concat(xxs):

   The concatenation of all the elements in a list.
   xs = list(chain.from_iterable(xxs))
   unit =  if isinstance(xs, str) else []
   return unit if not xs else (
       .join(xs) if isinstance(xs[0], str) else xs
   )

1. enumFromTo :: (Int, Int) -> [Int]

def enumFromTo(m):

   Integer enumeration from m to n.
   return lambda n: list(range(m, 1 + n))

1. findIndex :: (a -> Bool) -> [a] -> Maybe Int

def findIndex(p):

   Just the first index at which an
      element in xs matches p,
      or Nothing if no elements match.
   def go(xs):
       try:
           return Just(next(
               i for i, x in enumerate(xs) if p(x)
           ))
       except StopIteration:
           return Nothing()
   return lambda xs: go(xs)

1. fst :: (a, b) -> a

def fst(tpl):

   First member of a pair.
   return tpl[0]

1. fTable :: String -> (a -> String) ->
2. (b -> String) ->
3. (a -> b) -> [a] -> String

def fTable(s):

   Heading -> x display function -> fx display function ->
         f -> value list -> tabular string.
   def go(xShow, fxShow, f, xs):
       w = max(map(compose(len)(xShow), xs))
       return s + '\n' + '\n'.join([
           xShow(x).rjust(w, ' ') + ' -> ' + fxShow(f(x)) for x in xs
       ])
   return lambda xShow: lambda fxShow: (
       lambda f: lambda xs: go(
           xShow, fxShow, f, xs
       )
   )

1. iterate :: (a -> a) -> a -> Gen [a]

def iterate(f):

   An infinite list of repeated
      applications of f to x.
   def go(x):
       v = x
       while True:
           yield v
           v = f(v)
   return lambda x: go(x)

1. maybe :: b -> (a -> b) -> Maybe a -> b

def maybe(v):

   Either the default value v, if m is Nothing,
      or the application of f to x,
      where m is Just(x).
   return lambda f: lambda m: v if m.get('Nothing') else (
       f(m.get('Just'))
   )

1. snd :: (a, b) -> b

def snd(tpl):

   Second member of a pair.
   return tpl[1]

1. take :: Int -> [a] -> [a]
2. take :: Int -> String -> String

def take(n):

   The prefix of xs of length n,
      or xs itself if n > length xs.
   return lambda xs: (
       xs[0:n]
       if isinstance(xs, list)
       else list(islice(xs, n))
   )

1. uncurry :: (a -> b -> c) -> ((a, b) -> c)

def uncurry(f):

   A function over a tuple
      derived from a default or
      curried function.
   return lambda xy: f(xy[0], xy[1])

1. until :: (a -> Bool) -> (a -> a) -> a -> a

def until(p):

   The result of repeatedly applying f until p holds.
      The initial seed value is x.
   def go(f, x):
       v = x
       while not p(v):
           v = f(v)
       return v
   return lambda f: lambda x: go(f, x)

1. MAIN ---

if __name__ == '__main__':

   main()</lang>

Output:

First cycle detected, if any:

           cycle([1, 2, 3]) -> [1,2,3] (from:0, length:3)
 [0..10000] + cycle([1..8]) -> [1,2,3,4,5,6,7,8] (from:10001, length:8)
                 [1..10000] -> No cycle found
f(x) = (x*x + 1) modulo 255 -> [101,2,5,26,167,95] (from:2, length:6)

Racket

I feel a bit bad about overloading λ, but it’s in the spirit of the algorithm. <lang racket>

1. lang racket/base

returns (values lambda mu)

(define (brent f x0)

 ;; main phase: search successive powers of two
 (define λ
   (let main-phase ((power 1)
                    (λ 1)
                    (tortoise x0)
                    (hare (f x0))) ;; f(x0) is the element/node next to x0.
     (cond [(= hare tortoise) λ]
           ;; time to start a new power of two?
           [(= power λ) (main-phase (* power 2) 1 hare (f hare))]
           [else (main-phase power (add1 λ) tortoise (f hare))])))
 
 (values
  λ
  ;; Find the position of the first repetition of length λ
  (let race ((µ 0)
             (tortoise x0)
             ;; The distance between the hare and tortoise is now λ.
             (hare (for/fold ((hare x0)) ((_ (in-range λ))) (f hare))))
    ;; Next, the hare and tortoise move at same speed until they agree
    (if (= tortoise hare) µ (race (add1 µ) (f tortoise) (f hare))))))

(module+ test

 (require rackunit racket/generator)
 (define (f x) (modulo (+ (* x x) 1) 255))
 (define (make-generator f x0)
   (generator () (let loop ((x x0)) (yield x) (loop (f x)))))
 
 (define g (make-generator f 3))
 
 (define l (for/list ((_ 20)) (g)))
 (check-equal? l '(3 10 101 2 5 26 167 95 101 2 5 26 167 95 101 2 5 26 167 95))  
 (displayln l)
 (let-values (([µ λ] (brent f 3)))
   (printf "Cycle length = ~a~%Start Index = ~a~%" µ λ)))

</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
Start Index = 2

Raku

(formerly 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)

REXX

Brent's algorithm

<lang rexx>/*REXX program detects a cycle in an iterated function [F] using Brent's algorithm. */ init= 3; $= init

                                 do  until length($)>79;    $= $  f( word($, words($) ) )
                                 end   /*until*/

say ' original list=' $ ... /*display original number list*/ call Brent init /*invoke Brent algorithm for F*/ parse var result cycle idx /*get 2 values returned from F*/ say 'numbers in list=' words($) /*display number of numbers. */ say ' cycle length =' cycle /*display the cycle to term*/ say ' start index =' idx " ◄─── zero index" /* " " index " " */ say 'cycle 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 value of 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   when using the defaults:

  original list= 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 ...
numbers in list= 27
  cycle length = 6
  start index  = 2   ◄─── zero index
cycle sequence = 101 2 5 26 167 95

sequential search algorithm

<lang rexx>/*REXX program detects a cycle in an iterated function [F] using a sequential search. */ x= 3; $= x /*initial couple of variables*/

                           do  until cycle\==0;   x= f(x) /*calculate another number.  */
                           cycle= wordpos(x, $)           /*This could be a repeat.    */
                           $= $  x                        /*append number to   $  list.*/
                           end   /*until*/

say ' original list=' $ ... /*display the sequence. */ say 'numbers in list=' words($) /*display number of numbers. */ say ' cycle length =' words($) - cycle /*display the cycle to term*/ say ' start index =' cycle - 1 " ◄─── zero based" /* " " index " " */ say 'cycle sequence =' subword($, cycle, words($)-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:

  original list= 3 10 101 2 5 26 167 95 101 ...
numbers in list= 9
  cycle length = 6
  start index  = 2   ◄─── zero based
cycle sequence = 101 2 5 26 167 95

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 program detects a cycle in an iterated function [F] using a hash table. */ !.= .;  !.x= 1; x= 3; $= x /*assign initial value. */

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

say ' original list=' $ ... /*maybe display the list.*/ say 'numbers in list=' # /*display number of nums.*/ say ' cycle length =' # - !.x /* " " cycle. */ say ' start index =' !.x - 1 " ◄─── zero based" /* " " index. */ say 'cycle sequence =' subword($, !.x, # - !.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 identical to the 2^nd REXX version.

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 49^th power   was chosen;   it
generates a cyclic sequence that contains over 1.5 million numbers. <lang rexx>/*REXX program detects a cycle in an iterated function [F] using a robust hashing. */ parse arg power . /*obtain optional args from C.L. */ if power== | power="," then power=8 /*Not specified? Use the default*/ numeric digits 500 /*be able to handle big numbers. */ divisor= 2**power - 1 /*compute the divisor, power of 2*/ numeric digits max(9, length(divisor) * 2 + 1) /*allow for the square plus one*/ say ' power =' power /*display the power to the term*/ say ' divisor =' "{2**"power'}-1 = ' divisor /* " " divisor. " " " */ say x=3; $=x; $$=; m=100;  !.=.;  !.x=1 /*M: maximum numbers to display.*/

              do n=1+words($);    x= f(x);  $$=$$ x
              if n//2000==0  then do;   $=$ $$;   $$=;  end   /*Is a 2000th N?  Rejoin.*/
              if !.x\==.     then leave               /*is this number a repeat?  Leave*/
              !.x= n
              end   /*n*/                             /*N:   is the size of   $   list.*/

$= space($ $$) /*append residual numbers to $ */ if n<m then say ' original list=' $ ... /*maybe display the list to term.*/

            say 'numbers in list=' n                  /*display number of numbers.     */
            say '  cycle length =' n - !.x            /*display the cycle to the term. */
            say '  start index  =' !.x - 1   "  ◄─── zero based"      /*show the index.*/

if n<m then say 'cycle sequence =' subword($, !.x, n- !.x) /*maybe display the 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   when the input (power of two) used is:     49

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

  original list= 3 10 101 2 5 26 167 95 101 ...
numbers in list= 9
  cycle length = 6
  start index  = 2   ◄─── zero based
cycle sequence = 101 2 5 26 167 95 

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

Scala

Procedural

Output:

Best seen in running your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).

<lang Scala>object CycleDetection extends App {

 def brent(f: Int => Int, x0: Int): (Int, Int) = {
   // main phase: search successive powers of two
   // f(x0) is the element/node next to x0.
   var (power, λ, μ, tortoise, hare) = (1, 1, 0, x0, f(x0))

   while (tortoise != hare) {
     if (power == λ) { // time to start a new power of two?
       tortoise = hare
       power *= 2
       λ = 0
     }
     hare = f(hare)
     λ += 1
   }

   // Find the position of the first repetition of length 'λ'
   tortoise = x0
   hare = x0
   for (i <- 0 until λ) hare = f(hare)

   // The distance between the hare and tortoise is now 'λ'.
   // Next, the hare and tortoise move at same speed until they agree
   while (tortoise != hare) {
     tortoise = f(tortoise)
     hare = f(hare)
     μ += 1
   }
   (λ, μ)
 }

 def cycle = loop.slice(μ, μ + λ)

 def f = (x: Int) => (x * x + 1) % 255

 // Generator for the first terms of the sequence starting from 3
 def loop: LazyList[Int] = 3 #:: loop.map(f(_))

 val (λ, μ) = brent(f, 3)
 println(s"Cycle length = $λ")
 println(s"Start index  = $μ")
 println(s"Cycle        = ${cycle.force}")

}</lang>

Functional

<lang Scala> import scala.annotation.tailrec

object CycleDetection {

 def brent(f: Int => Int, x0: Int): (Int, Int) = {
   val lambda = findLambda(f, x0)
   val mu = findMu(f, x0, lambda)
   (lambda, mu)
 }

 def cycle(f: Int => Int, x0: Int): Seq[Int] = {
   val (lambda, mu) = brent(f, x0)
   (1 until mu + lambda)
     .foldLeft(Seq(x0))((list, _) => f(list.head) +: list)
     .reverse
     .drop(mu)
 }

 def findLambda(f: Int => Int, x0: Int): Int = {
   findLambdaRec(f, tortoise = x0, hare = f(x0), power = 1, lambda = 1)
 }

 def findMu(f: Int => Int, x0: Int, lambda: Int): Int = {
   val hare = (0 until lambda).foldLeft(x0)((x, _) => f(x))
   findMuRec(f, tortoise = x0, hare, mu = 0)
 }

 @tailrec
 private def findLambdaRec(f: Int => Int, tortoise: Int, hare: Int, power: Int, lambda: Int): Int = {
   if (tortoise == hare) {
     lambda
   } else {
     val (newTortoise, newPower, newLambda) = if (power == lambda) {
       (hare, power * 2, 0)
     } else {
       (tortoise, power, lambda)
     }
     findLambdaRec(f, newTortoise, f(hare), newPower, newLambda + 1)
   }
 }

 @tailrec
 private def findMuRec(f: Int => Int, tortoise: Int, hare: Int, mu: Int): Int = {
   if (tortoise == hare) {
     mu
   } else {
     findMuRec(f, f(tortoise), f(hare), mu + 1)
   }
 }

 def main(args: Array[String]): Unit = {
   val f = (x: Int) => (x * x + 1) % 255
   val x0 = 3
   val (lambda, mu) = brent(f, x0)
   val list = cycle(f, x0)

   println("Cycle length = " + lambda)
   println("Start index  = " + mu)
   println("Cycle        = " + list.mkString(","))
 }

} </lang>

Output:

Cycle length = 6
Start index  = 2
Cycle        = 101,2,5,26,167,95

Sidef

Translation of: Raku

<lang ruby>func brent (f, x0) {

   var power = 1
   var λ = 1
   var tortoise = x0
   var hare = f(x0)

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

   var μ = 0
   tortoise = x0
   hare = x0
   { hare = f(hare) } * λ

   while (tortoise != hare) {
       tortoise = f(tortoise)
       hare = f(hare)
       μ += 1
   }

   return (λ, μ)

}

func cyclical_function(x) { (x*x + 1) % 255 }

var (l, s) = brent(cyclical_function, 3)

var seq = gather {

   var x = 3
   { take(x); x = cyclical_function(x) } * 20

}

say seq.join(', ')+', ...'

say "Cycle length #{l}."; say "Cycle start index #{s}." say [seq[s .. (s + l - 1)]]</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]

Visual Basic .NET

Translation of: C#

<lang vbnet>Module Module1

   Function FindCycle(Of T As IEquatable(Of T))(x0 As T, yielder As Func(Of T, T)) As Tuple(Of Integer, Integer)
       Dim power = 1
       Dim lambda = 1
       Dim tortoise As T
       Dim hare As T

       tortoise = x0
       hare = yielder(x0)

       ' Find lambda, the cycle length
       While Not tortoise.Equals(hare)
           If power = lambda Then
               tortoise = hare
               power *= 2
               lambda = 0
           End If
           hare = yielder(hare)
           lambda += 1
       End While

       ' Find mu, the zero-based index of the start of the cycle
       Dim mu = 0
       tortoise = x0
       hare = x0
       For times = 1 To lambda
           hare = yielder(hare)
       Next

       While Not tortoise.Equals(hare)
           tortoise = yielder(tortoise)
           hare = yielder(hare)
           mu += 1
       End While

       Return Tuple.Create(lambda, mu)
   End Function

   Sub Main()
       ' A recurrence relation to use in testing
       Dim sequence = Function(_x As Integer) (_x * _x + 1) Mod 255

       ' Display the first 41 numbers in the test series
       Dim x = 3
       Console.Write(x)
       For times = 0 To 39
           x = sequence(x)
           Console.Write(",{0}", x)
       Next
       Console.WriteLine()

       ' Test the FindCycle method
       Dim cycle = FindCycle(3, sequence)
       Console.WriteLine("Cycle length = {0}", cycle.Item1)
       Console.WriteLine("Start index = {0}", cycle.Item2)
   End Sub

End Module</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

Wren

Working from the code in the Wikipedia article: <lang ecmascript>var brent = Fn.new { |f, x0|

   var lam = 1
   var power = 1
   var tortoise = x0
   var hare = f.call(x0)
   while (tortoise != hare) {
       if (power == lam) {
           tortoise = hare
           power = power * 2
           lam = 0
       }
       hare = f.call(hare)
       lam = lam + 1
   }
   tortoise = hare = x0
   for (i in 0...lam) hare = f.call(hare)
   var mu = 0
   while (tortoise != hare) {
       tortoise = f.call(tortoise)
       hare = f.call(hare)
       mu = mu + 1
   }
   return [lam, mu]

}

var f = Fn.new { |x| (x*x + 1) % 255 } var x0 = 3 var x = x0 var seq = List.filled(21, 0) // limit to first 21 terms say for (i in 0..20) {

   seq[i] = x
   x = f.call(x)

} var res = brent.call(f, x0) var lam = res[0] var mu = res[1] System.print("Sequence = %(seq)") System.print("Cycle length = %(lam)") System.print("Start index = %(mu)") System.print("Cycle = %(seq[mu...mu+lam])")</lang>

Output:

Sequence     = [3, 10, 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
Cycle        = [101, 2, 5, 26, 167, 95]

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