#!/usr/bin/perl

$, = " "; 
$\ = "\n"; 

##                 PA  MN  NM  NH  IA  WI  OH  FL  NV  CO  NC  VA  WV  AR  MO  AZ  ME  OR  WA  MI
my @state_evs   = (21, 10,  5,  4,  7, 10, 20, 27,  5,  9, 15, 13,  5,  6, 11, 10,  4,  7, 11, 17);
my @state_polls = ( 6,  0,  2,  2,  0,  0, -2, -2, -7, -5, -3, -4, -5, -5, -6, -7,  6,  7,  7,  5);
my @state_moes  = ( 4,  4,  4,4.5,  4,  4,3.5,2.8,4.1,4.5,  4,3.9,  4,  4,3.9,  4,3.9,  4,  4,  4);
my $ev_dem_safe = 168; 

my $number_states = $#state_evs + 1;
my $possibilities = (1 << $number_states);

my @stdevs = (); 
my @state_prob  = (); 
foreach my $state_index (0..($number_states-1)) {
  my $diff = $state_polls[$state_index]; 
  my $moe = $state_moes[$state_index]; 
  my $standard_deviations = $diff/$moe; 
  push @stdevs, $standard_deviations; 
  push @state_prob, N($standard_deviations); 
}

my @state_space = ();
my @prob = (1) x $possibilities;
foreach my $state_index (0..($number_states-1)) {
   my $state_bit = 1 << $state_index;
   my $state_ev = $state_evs[$state_index];
   my $p = $state_prob[$state_index]; 
   foreach (0..($possibilities-1)) {
      $state_space[$_] += ($_ & $state_bit)? $state_ev: 0; 
      $prob[$_] *= ($_ & $state_bit)? $p: 1-$p; 
   }
}

my ($prob_lose, $prob_tie, $prob_win) = (0,0,0); 
foreach (0..$#state_space) {
   $prob_lose += $prob[$_] if $ev_dem_safe + $state_space[$_] < 269; 
   $prob_tie  += $prob[$_] if $ev_dem_safe + $state_space[$_] == 269; 
   $prob_win  += $prob[$_] if $ev_dem_safe + $state_space[$_] > 269; 
}
printf "lose %4.2f%%, tie %4.2f%%, win %4.2f%%\n", 
	100*$prob_lose, 100*$prob_tie, 100*$prob_win; 


###########################################################
###########################################################

sub N {
  my $x1 = $_[0]; 

  my $x = abs($x1); 
  my $index = ($x1 < 0)? -1: 0; 
  my $sqrt_two_pi = 2.50662827463100024161235523934;
  my $ph = exp(-$x*$x/2.0) / $sqrt_two_pi;

  if( $x <= 3.0 ) {
     $w = $x;
     $sd = 10.0 * (1.0 + $x);
     $s = int($sd);
     for( ; $s > 0; $s--) {
         $w = $w*$x*$x/(2*$s+1) + $x;
     }

     $r = $ph * $w + 0.5;
  } else {
    if( $x <= 10.0 ) {
       $sd = -10.0 + 85.0/$x + 255.0/$x/$x;
       $s = int($sd);
       $w = $x * $x;
       for( ; $s > 0; $s--) {
           $w = $x*$x + 4*$s - 1 - 2*$s * (2*$s+1)/$w;
       }
       $r = 1 + $ph * (1/$w-1)/$x;
    } else {
     $r = 1;
    }
  }
  return ($index)? 1.0-$r: $r;
}

