#!/usr/bin/perl

#	sdice <n>
#	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>\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;
}