/* 

  Euler #14 in Picat.

  Problem 14
  """
  The following iterative sequence is defined for the set of positive integers:

  n n/2 (n is even)
  n 3n + 1 (n is odd)

  Using the rule above and starting with 13, we generate the following 
  sequence:
  13 40 20 10 5 16 8 4 2 1

  It can be seen that this sequence (starting at 13 and finishing at 1) 
  contains 
  10 terms. Although it has not been proved yet (Collatz Problem), it is 
  thought that all starting numbers finish at 1.

  Which starting number, under one million, produces the longest chain?

  NOTE: Once the chain starts the terms are allowed to go above one million.)
  """

  This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/


main => go.

go => time(euler14).


%
% Same as euler14d/0 (NFZ's version) but the Chain argument is skipped
% since it's not really needed for this problem.
%
% It's significantly faster than euler14d: 0.693s vs 1.232s.
% (Thanks to Jan Wielemaker for the idea!)
%
% testing all numbers: 0.9s
% only the odd numbers: 0.7s
%
euler14 =>
    garbage_collect(300_000_000),
    max_chain(N,Len),
    printf("N=%w,Len=%w%n",N,Len).

table (-,max)
max_chain(N,Len) =>
    % between(2,999999,N), % Original
    between(3,2,999999,N), % just checking the odd numbers
    % member(N,3..2..999999), % slightly slower
    gen(N,Len).

table (+,-)
gen(1,Len) => Len=1.
gen(N,Len), N mod 2 ==0 => 
    gen(N div 2,Len1),
    Len=Len1+1.
gen(N,Len) =>
    gen(3*N+1,Len1),
    Len=Len1+1.


    
between(From,Step,To,N) =>
  between(From,To div Step,Tmp),
  N = (Tmp-From)*Step+From.



% This is 2.6s (1.96 with only odd numbers
euler14a => 
   [MaxLen, MaxN] = longest_seq(999999),
   writeln([maxLen=MaxLen, maxN=MaxN]),
   nl.

% 8.9s
euler14b => 
   [MaxLen, MaxN] = longest_seq2(999999),
   writeln([maxLen=MaxLen, maxN=MaxN]),
   nl.

% 2.78s
euler14c =>
   L=[collLength(N)=N : N in 1..10**6],
   Max=max(L),
   writeln(max=Max),
   nl.




% This is from
% http://picat-lang.org/euler/p14.pi
% (slightly adjusted)
% testing all numbers: 1.4s
% only the odd numbers: 1.2s
euler14d =>
    garbage_collect(300_000_000),
    max_chain_d(N,_Chain,Len),
    printf("N=%w,Len=%w%n",N,Len).

table (-,-,max)
max_chain_d(N,Chain,Len) =>
    % between(2,999999,N), % Original
    between(3,2,999999,N), % just checking the odd numbers
    % member(N,3..2..999999), % slightly slower
    gen_d(N,Chain,Len).

table (+,-,-)
gen_d(1,Chain,Len) => Chain=[1], Len=1.
gen_d(N,Chain,Len), N mod 2 ==0 => 
    gen_d(N div 2,Chain1,Len1),
    Chain=[N|Chain1],
    Len=Len1+1.
gen_d(N,Chain,Len) =>
    gen_d(3*N+1,Chain1,Len1),
    Chain=[N|Chain1],
    Len=Len1+1.


% We just care about the lengths so 
% map will do fine.
% Testing all numbers 2..Limit: 2.5s
% Testing just the odd numbers: 1.96s
longest_seq(Limit) = [MaxLen, MaxN] =>

   Lens = new_map(),
   MaxLen = 0,
   MaxN = 1,
   foreach(N in 3..2..Limit) 
      M = N,
      CLen = 1,
      while (M > 1) 
        (
         Lens.has_key(M) -> 
            CLen := CLen + Lens.get(M) - 1,
            M := 1
         ; 
            M := hailstone1(M), % 2.5s
            % M := hailstone2(M),  % 3.1s
            % Without the call to hailstone1/1: 2.7s
            % (
            %    M mod 2 == 0 ->
            %    % M /\ 1 == 0 ->
            %       M := M // 2
            %       % M := M >> 1
            %    ;
            %       M := 3*M+1
            % ),
            CLen := CLen + 1
         )
      end,
      Lens.put(N, CLen),
      (
      CLen > MaxLen ->
         MaxLen := CLen,
         MaxN := N
        ;
          true
      )
   end.
   


% Without caching the lengths
longest_seq2(Limit) = [MaxLen, MaxN] =>
   MaxLen = 0,
   MaxN = 1,
   N = 1,
   while( N < Limit) 
      M = N,
      CLen = 1,
      while (M > 1) 
         M := hailstone1(M),
         CLen := CLen + 1
      end,
      if CLen > MaxLen then
         MaxLen := CLen,
         MaxN := N
      end,
      N := N + 1 
   end.


% With table: 6.2s
% Without table: 2.5s
% table
hailstone1(N) = N//2, N mod 2 == 0 => true.
hailstone1(N) = 3*N+1 => true.


% Alternative, a bit slower
% table
hailstone2(N) = H => 
   (N mod 2 == 0 ->
      H = N // 2
   ;
      H := 3*N+1
   ).

table
hailstoneseq1(N) = Seq =>
   Seq := [N],
   while (N > 1)
      N := hailstone1(N),
      Seq := Seq ++ [N]
   end.

table
hailstoneseq2(N) = Seq =>
   Seq := [N],
   while (N > 1)
      N := hailstone2(N),
      Seq := Seq ++ [N]
   end.


table
collLength(1) = 1.
collLength(2) = 2.
collLength(N) = 1+ cond(N mod 2 == 0,
                        collLength(N div 2),
                        collLength(3*N +1)).