# Copyright: Erich Kaltofen, 2006
cycle_indices := proc(x0, nxt, equal)
# Finds the preperiod length mu and the period length lambda of an eventually 
# cycling function nxt with respect to the predicate equal.

   local X, Y, Z, i, n, lambda, mu;

      # R. Floyd's trick:  First find n such that
      #
      #      equal(nxt^(n)(x0), nxt^(2*n)(x0) ).
      #
      #   Then n is the first exponent in the period that is divisible
      #   by lambda.
      if printlevel >= 2 then print(x0); # , nxt, equal);
      fi;

      X := x0;  Y := x0;
      for i from 0 by 1 do
          # At this point X = nxt^(i)(x0), Y = nxt^(2*i)(x0), and
          #   not( equal( nxt^(j)(x0), nxt^(2j)(x0) ) ) for all j < i. *)

          X := nxt(X); Y := nxt( nxt(Y) );
          if printlevel >= 2 then print("X=", X, "Y=", Y); fi;
          if equal(X, Y) then n := i; break; fi;
      od;
      if printlevel >= 2 then print("n = ", n); fi;

      # Search for shortest repetion in period
      Z := X; lambda := 0;
      for i from 0 by 1 while( lambda = 0) do
          # At this point Z  = nxt^(n+i)(x0)

          Z := nxt(Z);
          if equal(X, Z) then lambda := i+1; fi;
      od;   

      # Search for the preperiod
      Y := x0; Z := x0; 
      for i from 0 by 1 while(i < lambda) do
          Z := nxt(Z);
      od;
      
      for mu from 0 by 1 while(not( equal(Y, Z) )) do
          # At this point we have Y = nxt^(mu)(x0) and
          # Z = nxt^(mu+lambda)(x0).

          Y := nxt(Y); Z := nxt(Z);
      od;

      RETURN([mu, lambda]);
end:

readlib(procmake): readlib(procbody):
pollard_rho := proc(n)
# Demonstration of the Pollard rho method for integer factoring
   local x,y,g, randomf, mygcd, pl, mu_lam, x_0;
      randomf := proc(x)
                 return [x[1]^2 + 1 mod x[2], x[2]];
                 end;
      mygcd := proc(x, y)
               # assumes x[2] = y[2] = N
               local g;
               g := igcd(x[1]-y[1], x[2]);
               if g > 1
                  then printf(`GCD =%d\n`,g); RETURN(true);
                  else RETURN(false);
               fi;
              end;
      pl := printlevel; printlevel := 3;
      x_0:=3; if nargs > 1 then x_0 := args[2]; fi;
      mu_lam := cycle_indices([x_0,n], randomf, mygcd);
      printlevel := pl;
      RETURN(mu_lam);
end:

# this code creates the online documentation for cycle_indices
`help/text/cycle_indices` := TEXT(
`FUNCTION: cycle_indices - compute the lengths of preperiod and period`,
`                          of a user supplied function`,
`   `,
`CALLING SEQUENCES: cycle_indices(x0, fn, eq)`,
`   `,
`PARAMETERS: x0 - a starting value`,
`            fn - a single argument function`,
`            eq - an equality test`,
`   `,
`SYNOPSIS:   `,
`- The call cycle_indices(x0, fn, eq) returns a pair [mu, lambda]`,
`  where mu is the length of the preperiod and lambda is the`,
`  length of the period of the sequence x0, fn(x0), fn(fn(x0)),...`,
`  where equality is tested by the function eq`,
`   `,
`EXAMPLE:   `,
`> nxt := proc(x) x^2 + 1 mod 101; end:`,
`> eq := proc(x, y) x = y; end:`,
`> cycle_indices(3, nxt, eq);`,
`   `,
`                                    [12, 9]`
):

`help/text/pollard_rho` := TEXT(
`FUNCTION: pollard_rho - demonstrate the Pollard rho`,
`                        factorization algorithm`,
`   `,
`CALLING SEQUENCES: pollard_rho(N)`,
`   `,
`PARAMETERS: N - an integer to be factored`,
`   `,
`SYNOPSIS:   `,
`- The call pollard_rho(N) returns a pair [mu, lambda]`,
`  where mu is the length of the preperiod and lambda is the`,
`  length of the period of the Pollard rho factoring algorithm`,
`  Intermediate results including the factor found are traced`,
`   `,
`EXAMPLE:   `,
`> pollard_rho(ithprime(100)*ithprime(101));`,
`         3, proc(x) (x^2+1) mod 295927 end, proc(x,y)`,
`                                            local g;`,
`                                                g := igcd(x-y,295927);`,
`                                                if 1 < g then`,
`                                                    printf(``GCD =%d``,g);`,
`                                                    RETURN(true)`,
`                                                else RETURN(false)`,
`                                                fi`,
`                                            end`,
` `,
`X := 10 Y := 101`,
`X := 101 Y := 210428`,
`X := 10202 Y := 198611`,
`X := 210428 Y := 75543`,
`X := 90248 Y := 266320`,
`X := 198611 Y := 217919`,
`X := 148003 Y := 111456`,
`X := 75543 Y := 71335`,
`X := 88582 Y := 151922`,
`X := 266320 Y := 8966`,
`X := 38676 Y := 2563`,
`X := 217919 Y := 170139`,
`X := 101164 Y := 47605`,
`GCD =541`,
`n = 12`,
`GCD =541`,
`GCD =541`,
`                                    [1, 13]`
):