/* 

  Partition function in Picat.

  https://rosettacode.org/wiki/Partition_function_P
  """
  The Partition Function P, often notated P(n) is the number of solutions where n∈ℤ can be 
  expressed as the sum of a set of positive integers.

  Example

  P(4) = 5   because   4 = Σ(4) = Σ(3,1) = Σ(2,2) = Σ(2,1,1) = Σ(1,1,1,1)


  P(n) can be expressed as the recurrence relation:

  P(n) = P(n-1) +P(n-2) -P(n-5) -P(n-7) +P(n-12) +P(n-15) -P(n-22) -P(n-26) +P(n-35) +P(n-40) ...

  The successive numbers in the above equation have the differences:   
     1, 3, 2, 5, 3, 7, 4, 9, 5, 11, 6, 13, 7, 15, 8 ...

  This task may be of popular interest because Mathologer made the video, The hardest 
  "What comes next?" (Euler's pentagonal formula), where he asks the programmers among his 
  viewers to calculate P(666). The video has been viewed more than 100,000 times in 
  the first couple of weeks since its release.

  In Wolfram Language, this function has been implemented as PartitionsP.


  Task

  - Write a function which returns the value of PartitionsP(n). Solutions can be iterative or recursive.

  - Bonus task: show how long it takes to compute PartitionsP(6666). 
  """


  Output:

  [1,2,3,5,7,11,15,22,30,42,56,77,101,135,176,231,297,385,490,627]

  p(666) = 11956824258286445517629485

  CPU time 0.005 seconds.

  p(1000) = 24061467864032622473692149727991

  CPU time 0.009 seconds.

  p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

  CPU time 0.22 seconds.

  p(10000) = 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144

  CPU time 0.447 seconds.


  p(66666) = 931283431869095717238416063534148471363928685651267074563360050977549251436058076515262033789845457356081278451246751375974038318319810923102769109446980055567090089060954748065131666952830623286286824837240058805177319865230913417587625830803675380262765598632742903192085254824621

  CPU time 17.594 seconds.


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

*/

% import v3_utils.
% import util.
% import cp.


main => go.

go ?=>
  println([partition1(I) : I in 1..20]),
  nl,
  
  time(println('p(666)'=partition1(666))),
  
  initialize_table, % clear the cache
  time(println('p(1000)'=partition1(1000))),
  
  initialize_table,
  time(println('p(6666)'=partition1(6666))),
  
  initialize_table,
  time(println('p(10000)'=partition1(10000))),    

  initialize_table,
  time(println('p(66666)'=partition1(66666))),    

  nl.
go => true.

% Inspired by the Maple approach
table
partition1(0) = 1.
partition1(N) = P =>
  S = 0,
  K = 1,
  M = (K*(3*K-1)) // 2,  
  while (M <= N)
     M := (K*(3*K-1)) // 2,  
     S := S - ((-1)**K)*partition1(N-M),
     K := K + 1
  end,
  K := 1,
  M := (K*(3*K+1)) // 2,
  while (M <= N)
     M := (K*(3*K+1)) // 2,  
     S := S - ((-1)**K)*partition1(N-M),
     K := K + 1
  end,
  P = S.