#!/usr/bin/perl
#require 5.008;
our $VERSION = "0.00"; #Time-stamp: <2011-11-18T10:59:41Z>

use strict;
use warnings;
use utf8; # Japanese

use Encode;

our $CHECK_LOG = 0;

sub Perm {
  my ($n, $k) = @_;
  my $r = 1;
  for (my $i = 0; $i < $k; $i++) {
    $r *= ($n - $i);
  }
  return $r;
}

sub Comb {
  my ($n, $k) = @_;
  if ($n - $k < $k) {
    $k = $n - $k;
  }
  return Perm($n, $k) / Perm($k, $k);
}

sub bits {
  my ($n) = @_;
  my $r = 0;
  while ($n) {
    if ($n % 2) {
      $r++;
    }
    $n = $n >> 1;
  }
  return $r;
}

MAIN:
{
  our @CARDS = (0 .. 21);
  splice(@CARDS, 13, 1);
  splice(@CARDS, 0, 1);
  my $N = scalar @CARDS;

  my @p;

  for (my $odd = 0; $odd < 7; $odd++) {
    my $even = 6 - $odd;
    $p[$odd] = {};
    if ($odd == 0 || $even == 0) {
      $p[$odd]->{ex8ex11} = 0;
    } else {
      $p[$odd]->{ex8ex11} = Perm($N/2 - 1, $odd - 1) * Perm($N/2 - 1, $even - 1)
	/ Perm($N - 2, 4) * ($even / $N) * ($odd / ($N - 1));
    }

    $p[$odd]->{ne8ne11} = Perm($N/2 - 1, $odd) * Perm($N/2 - 1, $even)
      / Perm($N, 6);

    if ($even == 0) {
      $p[$odd]->{ex8ne11} = 0;
    } else {
      $p[$odd]->{ex8ne11} = Perm($N/2 - 1, $odd) * Perm($N/2 - 1, $even - 1)
	/ Perm($N - 1, 5) * ($even / $N);
    }

    if ($odd == 0) {
      $p[$odd]->{ne8ex11} = 0;
    } else {
      $p[$odd]->{ne8ex11} = Perm($N/2 - 1, $odd - 1) * Perm($N/2 - 1, $even)
	/ Perm($N - 1, 5) * ($odd / $N);
    }
  }

  {
    my $sum = 0;
    for (my $i = 0; $i < 7; $i++) {
      my $lsum = 0;
      foreach my $k (sort keys %{$p[$i]}) {
	print "p[$i]->{$k} == $p[$i]->{$k}\n";
	$sum += $p[$i]->{$k};
	$lsum += $p[$i]->{$k};
      }
      print "p[$i] == $lsum\n";
    }
    print "$sum\n";
    print "" . (Perm(10, 3) * Perm(10, 3) / Perm(20, 6)) . "\n";
    print "" . (Perm(10, 2) * Perm(10, 4) / Perm(20, 6)) . "\n";
    print "" . (Perm(10, 1) * Perm(10, 5) / Perm(20, 6)) . "\n";
    print "" . (Perm(10, 6) / Perm(20, 6)) . "\n";
    print "" . Comb(6, 3) * (Perm(10, 3) * Perm(10, 3) / Perm(20, 6)) . "\n";
    print "" . Comb(6, 2) * (Perm(10, 2) * Perm(10, 4) / Perm(20, 6)) . "\n";
    print "" . Comb(6, 1) * (Perm(10, 1) * Perm(10, 5) / Perm(20, 6)) . "\n";
    print "" . Comb(6, 0) * (Perm(10, 6) / Perm(20, 6)) . "\n";
  }

  my @r;

  $r[3] = ($p[3]->{ex8ex11} * 2 * Comb(6, 3)
	   + $p[3]->{ne8ne11} * 2 * Comb(6, 3)
	   + $p[3]->{ex8ne11} * 1 * Comb(6, 3)
	   + $p[3]->{ne8ex11} * 1 * Comb(6, 3)) / 2 / Comb(6, 3);
  $r[2] = ($p[2]->{ex8ex11} * 2 * Comb(6, 2)
	   + $p[2]->{ne8ne11} * 2 * Comb(6, 2)
	   + $p[3]->{ne8ex11} * 1 * Comb(6, 3)
	   + $p[2]->{ex8ne11} * 2 * Comb(6, 2)) / 2 / Comb(6, 2);
  $r[1] = ($p[1]->{ex8ex11} * 2 * Comb(6, 1)
	   + $p[1]->{ne8ne11} * 2 * Comb(6, 1)
	   + $p[2]->{ne8ex11} * 2 * Comb(6, 2)
	   + $p[1]->{ex8ne11} * 2 * Comb(6, 1)) / 2 / Comb(6, 1);
  $r[0] = ($p[0]->{ex8ex11} * 2 * Comb(6, 0)
	   + $p[0]->{ne8ne11} * 2 * Comb(6, 0)
	   + $p[1]->{ne8ex11} * 2 * Comb(6, 1)
	   + $p[0]->{ne8ex11} * 2 * Comb(6, 0)
	   + $p[0]->{ex8ne11} * 2 * Comb(6, 0)) / 2 / Comb(6, 0);
  $r[4] = $r[2];
  $r[5] = $r[1];
  $r[6] = $r[0];

  my @r2;
  my @b2;

  if ($CHECK_LOG) {
    open(my $fh, "<", "calc_prob_2.log") or die $!;
    while (<$fh>) {
      if (/^([01-9a-f][01-9a-f]): ([01-9\.]+)/) {
	$r2[unpack("c", pack("H2", $1))] = 0 + $2;
      } elsif (/^([01-9]): ([01-9\.]+)/) {
	$b2[int($1)] = $2;
      }
    }
    close($fh);
  }

  for (my $i = 0; $i < 64; $i++) {
    if ($CHECK_LOG && $r2[$i] - $r[bits($i)] > 0.000001) {
      print "ERROR $i $r2[$i]\n";
    }
    printf("%02x: %g\n", $i, $r[bits($i)]);
  }

  for (my $i = 0; $i < 7; $i++) {
    printf("%01d: %g\n", $i, $r[$i] * Comb(6, $i));
  }
}