euler12.pi

/* 

  Euler #12 in Picat.

  Problem 12
  """
  The sequence of triangle numbers is generated by adding the natural numbers. 
  So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. 
  The first ten terms would be:

  1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

  Let us list the factors of the first seven triangle numbers:

       1: 1
       3: 1,3
       6: 1,2,3,6
      10: 1,2,5,10
      15: 1,3,5,15
      21: 1,3,7,21
      28: 1,2,4,7,14,28

  We can see that the 7th triangle number, 28, is the first triangle number 
  to have over five divisors.

  Which is the first triangle number to have over five-hundred divisors?")
  """

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

*/

import util.
main => time(go).

go =>  euler12.

% Using a map: 0.26s
euler12 => 
  garbage_collect(300_000_000),
  println("Problem 12: "),
  Len = 0,
  I = 0,
  TNum = 0,
  while (Len <= 500)
     I := I + 1,
     TNum := TNum + I,
     Map = new_map(),
     factorsMap(TNum,Map),
     Len := prod([V+1 : _=V in Map])
  end,
  println([i=I, tnum=TNum, len=Len]),
  nl.


% 0.313s
euler12a => 
  println("Problem 12: "),
  Len = 0,
  I = 0,
  TNum = 0,
  while (Len <= 500)
     I := I + 1,
     TNum := TNum + I,
     Len := prod([E+1 : E in collect4(factors(TNum))])
  end,
  
  println([i=I, tnum=TNum, len=Len]),
  nl.



% recursive version, slightly slower: 0.318s
euler12b => 
  e12b(Num,Len),
  println([num=Num,len=Len]).

e12b(Num,Len) =>
  e12b(0,0,Num,0,Len).

e12b(N,Num0,Num,Len0,Len) ?=>
  Len0 < 500,
  Num1 = Num0+N+1,
  Len1 = prod([E+1 : E in collect4(factors(Num1))]),
  e12b(N+1,Num1,Num,Len1,Len).

e12b(N,Num0,Num,Len0,Len) =>
   N > 500,
   Num=Num0,
   Len=Len0.

% table
collect2(A) = C => 
  C = [],
  foreach(I in A.remove_dups()) 
    % C := C ++ [[I, [J: J in A, J == I].length]  ]
    C := C ++ [[J: J in A, J == I].length]
  end.

% Variant (slightly slower)
collect3(A) = C => 
  M = new_map(),
  foreach(I in A) 
    M.put(I, M.get(I,0)+1)
  end,
  C = [[K,V] : K=V in M].

% a variant of collect2/1
% using sum/1 is slightly faster than using length/1
% Also, skipping the first variable in the result list.
% collect4(A) = [[I,[J: J in A, J == I].length] : I in A.remove_dups() ].
% collect4(A) = [[I,length([1: J in A, J == I])] : I in A.remove_dups() ].
collect4(A) = [sum([1: J in A, J == I]) : I in A.remove_dups() ].
% collect4(A) = [[1: J in A, J == I].len : I in A.remove_dups() ].

collect5(A) = [ [1 : V in A, V=UU ].len : UU in A.remove_dups()].



alldivisorsM(N,Div) = [Divisors,NewN] =>
   M = N,
   Divisors1 = [],
   while (M mod Div == 0) 
      Divisors1 := [Div|Divisors1],
      M := M div Div
   end,
   NewN = M,
   Divisors = Divisors1.

% table
factors(N) = Factors =>
     M = N,
     Factors1 = [],
     while (M mod 2 == 0) 
         Factors1 := Factors1 ++ [2],  
         M := M div 2 
     end,
     T = 3,
     while (M > 1, T < ceiling(sqrt(M)))
        if M mod T == 0 then
           [Divisors, NewM] = alldivisorsM(M, T),
           Factors1 := Factors1 ++ Divisors,
           M := NewM
        end,
        T := T + 2 
     end,
     if M > 1 then Factors1 := [M|Factors1] end,
     Factors = Factors1.


% factors as a map
factorsMap(N,Map) =>
     % Map = new_map(),
     M = N,
     while (M mod 2 == 0) 
         Map.put(2,Map.get(2,0)+1),
         M := M div 2 
     end,
     T = 3,
     while (M > 1, T < ceiling(sqrt(M)))
        while (M mod T == 0) 
          Map.put(T,Map.get(T,0)+1),
          M := M div T
        end,
        T := T + 2 
     end,
     if M > 1 then 
        Map.put(M,Map.get(M,0)+1) 
     end.



next_prime(N) = P =>
  M = N+2,
  while(not prime_cached(M))
    M := M+2
  end,
  P = M.

table
prime_cached(N) => prime(N).