/* 

  Limerick primes in Picat.

  From John D. Cook: Limerick Primes
  http://www.johndcook.com/blog/2011/03/08/limerick-primes
  """
  The other day, Futility Closet posted this observation:

      10102323454577 is the smallest 14-digit prime number that 
      follows the rhyme scheme of a Shakespearean sonnet 
      (ababcdcdefefgg).

  I posted this on AlgebraFact and got a lot of responses. One was from 
  Matt Parker who replied that 11551 was the smallest prime with a limerick 
  rhyme scheme.
  
  So how many limerick primes are there? Well, there aren’t many candidates. 
  A limerick prime has to have the form AABBA where A is an odd digit and B 
  is any digit other than A. So for each of five choices of A, there are 
  nine possible B’s. 
  """

  Below there are some other bases. Note: 0 is removed as a possible digit.

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

*/

import cp.


main => go.

%
% There are 8 different Limerick primes:
%
% [11551,[1,1,5,5,1]]
% [33113,[3,3,1,1,3]]
% [33223,[3,3,2,2,3]]
% [33773,[3,3,7,7,3]]
% [77447,[7,7,4,4,7]]
% [77557,[7,7,5,5,7]]
% [99119,[9,9,1,1,9]]
% [99559,[9,9,5,5,9]]

%
go ?=>
  Base = 10,
  All = find_all([P,X],limerick_prime(Base, P,X)),
  foreach(A in All)
    println(A)
  end,
  println(len=All.len),
  nl.

go => true.


%
% Sample runs for Bases 2..36
%
% [2 = no_limerick]
% [3 = no_limerick]
% [base = 4,983,[3,3,1,1,3],33113]
% [base = 5,811,[1,1,2,2,1],11221]
% [base = 6,1597,[1,1,2,2,1],11221]
% [base = 7,8291,[3,3,1,1,3],33113]
% [base = 8,23333,[5,5,4,4,5],55445]
% [base = 9,51307,[7,7,3,3,7],77337]
% [base = 10,33223,[3,3,2,2,3],33223]
% [base = 11,80657,[5,5,6,6,5],55665]
% [base = 12,157411,[7,7,1,1,7],77117]
% [base = 13,154523,[5,5,4,4,5],55445]
% [base = 14,454031,[11,11,6,6,11],BB66B]
% [base = 15,379927,[7,7,8,8,7],77887]
% [base = 16,629417,[9,9,10,10,9],99AA9]
% [base = 17,621799,[7,7,9,9,7],77997]
% [base = 18,777373,[7,7,5,5,7],77557]
% [base = 19,1238429,[9,9,10,10,9],99AA9]
% [base = 20,1852211,[11,11,10,10,11],BBAAB]
% [base = 21,2653741,[13,13,11,11,13],DDBBD]
% [base = 22,2210723,[9,9,13,13,9],99DD9]
% [base = 23,3217619,[11,11,10,10,11],BBAAB]
% [base = 24,4501213,[13,13,14,14,13],DDEED]
% [base = 25,4470061,[11,11,2,2,11],BB22B]
% [base = 26,5229911,[11,11,14,14,11],BBEEB]
% [base = 27,7176721,[13,13,16,16,13],DDGGD]
% [base = 28,9562127,[15,15,16,16,15],FFGGF]
% [base = 29,9523903,[13,13,14,14,13],DDEED]
% [base = 30,10894963,[13,13,15,15,13],DDFFD]
% [base = 31,14318543,[15,15,19,19,15],FFJJF]
% [base = 32,18397649,[17,17,14,14,17],HHEEH]
% [base = 33,23231029,[19,19,14,14,19],JJEEJ]
% [base = 34,26163359,[19,19,22,22,19],JJMMJ]
% [base = 35,26258417,[17,17,15,15,17],HHFFH]
% [base = 36,32823163,[19,19,18,18,19],JJIIJ]
%
go2 ?=>
  S = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ",

  foreach(Base in 2..36)
    if limerick_prime2(Base, P,X) then
      println([base=Base,P,X,[S[I+1] : I in X]])
    else
      println([Base=no_limerick])
    end
  end,
  nl.

%
% Count solutions for different bases
% Base #solutions
% ---------------
% 4 = 1
% 5 = 1
% 6 = 5
% 7 = 4
% 8 = 7
% 9 = 7
% 10 = 8
% 11 = 11
% 12 = 10
% 13 = 14
% 14 = 12
% 15 = 12
% 16 = 29
% 17 = 21
% 18 = 25
% 19 = 27
% 20 = 27
% 21 = 26
% 22 = 32
% 23 = 25
% 24 = 35
% 25 = 37
% 26 = 37
% 27 = 36
% 28 = 47
% 29 = 47
% 30 = 57
% 31 = 52
% 32 = 70
% 33 = 70
% 34 = 67
% 35 = 67
% 36 = 77
%
go3 => 
  foreach(Base in 4..36)
    Count = count_all(limerick_prime(Base, _P,_X)),
    println(Base=Count)
  end,
  nl.
 

%
% First version. (See below for an alternate version.)
%
limerick_prime(Base, P,X) =>
  N = 5,
  
  P :: Base**(N-1)..Base**N-1,

  X = new_list(N),
  X :: 1..Base-1, % skip 0s

  is_prime(P),
  to_num(X,Base,P),


  % """
  % A limerick prime has to have the form AABBA where A is an 
  % odd digit and B is any digit other than A. 
  % """
  X[1] mod 2 #= 1,
  % AA..A
  X[1] #= X[2], X[1] #= X[5],
  % ..BB.
  X[3] #= X[4],
  X[1] #!= X[3],

  Vars = X ++ [P],
  solve([degree,split],Vars).
  % solve([rand],Vars).  

%
% Alternative approach
%
limerick_prime2(Base, P,X) =>
  N = 5,

  P :: Base**(N-1)..Base**N-1,

  X = new_list(N),
  X :: 1..Base-1, % skip 0s

  A :: 1..Base-1,
  B :: 1..Base-1,

  X = [A,A,B,B,A],
  A mod 2 #= 1,
  A #!= B,

  is_prime(P),
  to_num(X,Base,P),
  Vars = X ++ [P,A,B],
  solve([ff,updown],Vars).
  

is_prime(X) =>
   X #> 1 ,
   Max = fd_max(X),
   foreach(I in 2..1+ceiling(sqrt(Max)))
       I #< X #=> X mod I #> 0
   end.


%
% converts a number Num to/from a list of integer List given a base Base
%
to_num(List, Base, Num) =>
   Len = length(List),
   Num #= sum([List[I]*Base**(Len-I) : I in 1..Len]).