/*
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).
% This is from
% http://picat-lang.org/euler/p14.pi
% (slightly adjusted)
% testing all numbers: 1.7s
% only the odd numbers: 1.5s
euler14 =>
max_chain(N,_Chain,Len),
printf("N=%w,Len=%w%n",N,Len).
table (-,-,max)
max_chain(N,Chain,Len) =>
% between(2,999999,N), % Original
between(3,2,999999,N), % just checking the odd numbers
gen(N,Chain,Len).
table (+,-,-)
gen(1,Chain,Len) => Chain=[1], Len=1.
gen(N,Chain,Len), N mod 2 ==0 =>
gen(N div 2,Chain1,Len1),
Chain=[N|Chain1],
Len=Len1+1.
gen(N,Chain,Len) =>
gen(3*N+1,Chain1,Len1),
Chain=[N|Chain1],
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.
% 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.