(* :Title: WeirdDice and WeirdDiceNumber *) 
(* :Author: Christopher Moretti *) 
(* :Summary: WeirdDice[n] finds all pairs of weird n-sided dice.  WeirdDiceNumber[n] 
finds the number of pairs of weird n-sided dice, which is useful for exploring 
the distribution of how many there are. *) 
(* :History: Version 1.0 by Christopher Moretti, 1997 *) 
  

Clear[normaldie, dietolist, tack0, tack1, tack2, massmatrix, nonnegcheck, WeirdDice, 
WeirdDiceNumber] 

normaldie[n_] := normaldie[n] = Cancel[((x^n - 1)*x)/(x - 1)] 

dietolist[f_] := Block[{a, b}, a = CoefficientList[f, x]; 
    b = Table[Table[j, {i, 1, a[[j + 1]]}], {j, 0, Length[a] - 1}]; 
    If[nonnegcheck[a], Return[Flatten[b]], Return[{}]]] 
  
tack1[n_] := Flatten[Append[n, {1}]] 
tack0[n_] := Flatten[Append[n, {0}]] 
tack2[n_] := Flatten[Append[n, {2}]] 
  
massmatrix[1] := {{1}} 
massmatrix[n_] := massmatrix[n] = 
   Union[tack1 /@ massmatrix[n - 1], tack0 /@ massmatrix[n - 1], tack2 /@ massmatrix[n - 
1]] 

nonnegcheck[f_] :=  !MemberQ[Flatten[NonNegative[f]], False] 

WeirdDice[n_] := Block[{faclist, a, tot, 
b,die2faclist,die1faclist,coeff1,coeff2,die1,die2,i}, 
    faclist = Transpose[FactorList[normaldie[n]]][[1]]; a = Length[faclist]; tot = 
massmatrix[a]; 
     b = Length[tot]; Do[die1faclist = tot[[i]]; die2faclist = Table[2, {i, 1, a}] - die1faclist; 
       die1 = Expand[Product[faclist[[i]]^ die1faclist[[i]], {i, 1, a}]]; 
       die2 = Expand[Product[faclist[[i]]^ die2faclist[[i]], {i, 1, a}]]; dielist1 = 
dietolist[die1]; 
       dielist2 = dietolist[die2]; coeff1 = CoefficientList[die1, x]; 
       coeff2 = CoefficientList[die2, x]; 
       If[(die1 /. x -> 0) == 0 && (die1 /. x -> 1) == n && nonnegcheck[coeff1] && 
         (die2 /. x -> 0) == 0 && (die2 /. x -> 1) == n && nonnegcheck[coeff2] && 
          !NumberQ[die1 - die2], Print[dietolist[die1]]; Print[dietolist[die2]]; Print[" "]; ], 
      {i, 1, (b - 1)/2}]]; 

WeirdDiceNumber[n_] := 
   WeirdDiceNumber[n] = 
    Block[{faclist, a, tot, b, num,die2faclist,die1faclist,coeff1,coeff2,die1,die2,i}, faclist = 
Transpose[FactorList[normaldie[n]]][[1]]; 
      a = Length[faclist]; tot = massmatrix[a]; b = Length[tot]; num = 0; 
      Do[die1faclist = tot[[i]]; die2faclist = Table[2, {i, 1, a}] - die1faclist; 
        die1 = Expand[Product[faclist[[i]]^die1faclist[[i]], {i, 1, a}]]; 
        die2 = Expand[Product[faclist[[i]]^die2faclist[[i]], {i, 1, a}]]; 
        dielist1 = dietolist[die1]; dielist2 = dietolist[die2]; coeff1 = CoefficientList[die1, x]; 
        coeff2 = CoefficientList[die2, x]; 
        If[(die1 /. x -> 0) == 0 && (die1 /. x -> 1) == n && nonnegcheck[coeff2] && 
          (die2 /. x -> 0) == 0 && (die2 /. x -> 1) == n && nonnegcheck[coeff1] && 
           !NumberQ[die1 - die2], num = num + 1], {i, 1, (b - 1)/2}]; Return[num]];


WeirdDice::usage="WeirdDice[n] finds all pairs of weird dice with n sides."
WeirdDiceNumber::usage=
  "WeirdDiceNumber[n] finds the number of pairs of weird n-sided dice."
(* Examples
	WeirdDice[6]
WeirdDiceNumber[10]
*)