#!/usr/bin/perl # sdice # Col. George Sicherman. 1999-08-06. # # Print all pairs of n-sided dice that give the same sums # as standard n-sided dice. # Algorithm: # # 1. Factor x^n - 1 into cyclotomic polynomials. # 2. Recursively generate all partitions of two copies of # the primary factors x^([p-1][p^j]) + x^([p-2][p^j]) + ... + 1 # into two sets of equal size for each p. # 3. Recursivly generate all partitions of two copies of the # unitary factors (cyclotomics for composites) # into two sets, not necessarily of equal size. # 4. At the bottom of the recursion, the products of the two sets # of polynomials generate the spot assignments for the two dice. # Polynomials are represented with the high-order coefficient first. sub usage { die "usage: sdice \n" } use integer; &usage unless @ARGV==1; \$n = shift; &usage unless \$n && \$n !~ /\D/; # Storage: # # \$unitary[\$k] is the kth unitary factor. # \$primary[\$p][\$k] is the kth primary factor for prime \$p. for (\$i=2; \$i<=\$n; \$i++) { # Drop the x-1 factor. next if \$n % \$i; @poly = &cyclotomic(\$i); \$sum = 0; map {\$sum += \$_} @poly; if (1 == \$sum) { push @unitary, [@poly] } else { # Add the prime to the hash if it's not already there. \$phash{\$sum} = keys %phash unless defined \$phash{\$sum}; push @{\$primary[\$phash{\$sum}]}, [@poly]; } } &recurse(0,0,0,1,0,[1],[1]); # Main recursion. # # 0. Index of current prime. # 1. Number of factors assigned to first die for current prime. # 2. Number of factors assigned to second die for current prime. # 3. =1 if factors so far have been assigned evenly, 0 else. # 4. Number of pairs of unitary factors assigned. # 5. Cumulative product of polynomials for first die. # 6. Cumulative product of polynomials for second die. sub recurse { my (\$p, \$p1, \$p2, \$e, \$u, \$y1, \$y2) = @_; my (\$c, \$f, \$s, \$t, @r1, @r2); \$c = @{\$primary[\$p]}; # Number of distinct factors for p. if (\$p1 + \$p2 < 2*\$c) { # Assign next factor for this prime. \$f = (\$p1 + \$p2) / 2; # Next factor to apply. if (\$p1+2 <= \$c) { @r1 = &mul(\$y1, \$primary[\$p][\$f]); @r1 = &mul(\@r1, \$primary[\$p][\$f]); &recurse(\$p, \$p1+2, \$p2, 0, \$u, \@r1, \$y2); } if (\$p1+1 <= \$c && \$p2+1 <= \$c) { @r1 = &mul(\$y1, \$primary[\$p][\$f]); @r2 = &mul(\$y2, \$primary[\$p][\$f]); &recurse(\$p, \$p1+1, \$p2+1, \$e, \$u, \@r1, \@r2); } if (!\$e && \$p2+2 <= \$c) { @r2 = &mul(\$y2, \$primary[\$p][\$f]); @r2 = &mul(\@r2, \$primary[\$p][\$f]); &recurse(\$p, \$p1, \$p2+2, 0, \$u, \$y1, \@r2); } } elsif (\$p < @primary) { # Advance to next prime. &recurse(\$p+1, 0, 0, \$e, \$u, \$y1, \$y2); } elsif (\$u < @unitary) { # Primary factors have all have been assigned. # Assign next pair of unitary factors. for (\$s=0; \$s+\$e<=2; \$s++) { @r1 = @\$y1; for (\$t=\$s; \$t<2; \$t++) { @r1 = &mul(\@r1, \$unitary[\$u]); } @r2 = @\$y2; for (\$t=0; \$t<\$s; \$t++) { @r2 = &mul(\@r2, \$unitary[\$u]); } &recurse(\$p, \$p1, \$p2, \$e, \$u+1, \@r1, \@r2); } } else { # Unitary factors have all been assigned. # Print the result if legal. return if grep {\$_ < 0} (@\$y1, @\$y2); &dieprint(\$y1); &dieprint(\$y2); print "\n"; } } sub dieprint { my (\$i, \$j, \$sep, \$val, \$poly); \$poly = shift; \$sep = ""; for (\$i=@\$poly-1; \$i>=0; \$i--) { \$val = @\$poly-\$i; for (\$j=0; \$j<\$poly->[\$i]; \$j++) { print \$sep . \$val; \$sep = " "; } } print "\n"; } sub polyprint { # for debugging my (\$i, \$coeff, \$ans); \$ans = shift; for (\$i=0; \$i<@\$ans; \$i++) { next unless \$coeff = \$ans->[\$i]; print \$coeff < 0? "- ": \$i? "+ ": ""; print abs(\$coeff) if abs(\$coeff) != 1 || \$i+1==@\$ans; print "x" if \$i+1<@\$ans; print "^" . (@\$ans - \$i - 1) if \$i+2<@\$ans; print " " if \$i+1<@\$ans; } print "\n"; } sub cyclotomic { my (\$k, \$m, \$n, @answer, @numer, @denom); \$n = shift; # Closed form: for each divisor k of n, # multiply or divide by x^k - 1 according as mu(n/k) is + or -. for (\$k=\$n; \$k>=1; \$k--) { next if \$n % \$k; next unless \$m = &mu(\$n/\$k); if (\$m > 0) { push @numer, [1, (0) x (\$k-1), -1] } else { push @denom, [1, (0) x (\$k-1), -1] } } @answer = (1); for (@numer) { @answer = &mul(\@answer, \$_) } for (@denom) { @answer = &div(\@answer, \$_) } return @answer; } sub mul { my @ret = (); for (my \$p=0; \$p<@{\$_[0]}; \$p++) { for (my \$q=0; \$q<@{\$_[1]}; \$q++) { \$ret[\$p+\$q] += \$_[0][\$p] * \$_[1][\$q]; } } return @ret; } sub div { my @ret = (); for (my \$p = @{\$_[0]} - @{\$_[1]}; \$p >= 0; \$p--) { my \$q = \$_[0][0] / \$_[1][0]; push @ret, \$q; for (my \$s = 0; \$s < @{\$_[1]}; \$s++) { \$_[0][\$s] -= \$q * \$_[1][\$s]; } shift @{\$_[0]}; } return @ret; } sub mu { my \$n = shift; my \$last = 1; my \$ret = 1; for (my \$i=2; \$i<=\$n; \$i++) { last if \$n==1; next if \$n % \$i; return 0 if \$last == \$i; \$last = \$i; \$n /= \$i; \$ret = -\$ret; redo; } return \$ret; }