# Derived from the program provided by _Richard J. Mathar_ in the Mathar/Hasler link. 

fop := fopen("D:/SEQFAN/Maple/ppartable4.txt", WRITE, TEXT); 

ppartabl := proc(n, c) 
  local i, j, k, l, m, pless, p1, p2, n1, n2, pset1, pset2, 
        alln1n2done, foundp1p2, largestprime, nodd, nprimes, 
        primeparts; 

  # construct set of primes < n in pless. 
  pless := {}; 
  for i from 2 to n-1 do 
    if isprime(i) then 
      pless := `union`(pless, {i}); 
    end if; 
  end do: 

  if `mod`(n, 2) = 1 then 
    nodd := true; 
    largestprime := pless[-1]; 
  else 
    nodd := false; 
    # determine all prime 2-partitions of n. 
    nprimes := numelems(pless); 
    primeparts := {}; 
    for k from 2 to nprimes do 
      for l from nprimes by -1 to 2 do 
        if pless[k] <> pless[l] and n = pless[k]+pless[l] then 
          primeparts := `union`(primeparts, 
                                {pless[k], pless[l]}); 
          break; 
        end if; 
      end do; 
    end do; 
  end if; 

  j := 0; 
  pset2 := {}; 
  # loop over all non-empty subsets of the primes, P1 
  for pset1 in combinat[choose](pless, c) do 
    if 2 <= nops(pset1) then 
      # prevent duplication. 
      if pset1 = pset2 then 
        break; 
      end if; 
      if nodd = true then 
        if n = 2+largestprime and 
           (member(2, pset1) = true and 
            member(largestprime, pset1) = false) or 
           (member(2, pset1) = false and 
            member(largestprime, pset1) = true)) then 
          next; 
        end if; 
      else 
        for m to numelems(pset1) do 
          if (member(pset1[m], primeparts) = true and 
              member(pset1[-m], primeparts) = false) or 
             (member(pset1[m], primeparts) = false and 
              member(pset1[-m], primeparts) = true) then 
            next; 
          end if; 
        end do; 
      end if; 
      pset2 := `minus`(pless, pset1); 
      alln1n2done := true; 
      for n1 to n-1 do 
        n2 := n-n1; 
        foundp1p2 := false; 
        for p1 in pset1 do 
          if igcd(n1, p1) <> 1 then 
            foundp1p2 := true; 
            break; 
          end if; 
          for p2 in pset2 do 
            if igcd(n2, p2) <> 1 then 
              foundp1p2 := true; 
              break; 
            end if; 
          end do; 
          if foundp1p2 = true then 
            break; 
          end if; 
        end do; 
        if foundp1p2 = false then 
          alln1n2done := false; 
          break; 
        end if; 
      end do; 
      if alln1n2done = true then 
        j := j+1; 
        if j = 1 then 
          fprintf(fop, "%d ", n); 
        end if; 
        fprintf(fop, " %d %a %a ", j, pset1, pset2); 
      end if; 
    end if; 
  end do: 
end proc: 

with(numtheory): 

c := 2; 
lo := 3; 
hi := 1200; 

for n from lo to hi do 
  # If n-1 is a prime power then n cannot be prime-partionable. 
  if `not`(numelems(factorset(n-1)) = 1) then 
    ppartabl(n, c); 
  end if; 
end do: 

fclose(fop);