use Math::Libm qw(erf erfc M_SQRT2);
sub compare_bell_curves {
  my ($self,$m1,$sd1,$m2,$sd2) = @_;
  if ($sd1 > $sd2) {
    ($m1,$sd1,$m2,$sd2)=($m2,$sd2,$m1,$sd1);
  } elsif ($sd1 == $sd2) {
    # stupid corner case
    my $dist = abs($m1-$m2)/$sd1;
    return erfc($dist/2/M_SQRT2);
  }
  $m2 -= $m1;
  $m1 = 0;
  # Some terms omitted since $m1 = 0
  my $sd2s= $sd2*$sd2;
  my $sd1s= $sd1*$sd1;
  my $A = ($sd2s - $sd1s);
  my $B = 2*($m2*$sd1s);
  my $C = 2*(log($sd1)-log($sd2))*$sd1s*$sd2s - $m2*$m2*$sd1s;
  my $disc = $B*$B - 4*$A*$C;
  my $rdisc = sqrt($disc);
  my $lower = (-$B - $rdisc)/(2*$A);
  my $upper = (-$B + $rdisc)/(2*$A);
  my $p1 = 0.5 + erf(($lower-$m2)/$sd2/M_SQRT2)/2;
  my $p2 = (erf($upper/$sd1/M_SQRT2)-erf($lower/$sd1/M_SQRT2))/2;
  my $p3 = erfc(($upper-$m2)/$sd2/M_SQRT2)/2;
  $p1+$p2+$p3;
}


##</code><code>##

--
@/=map{[/./g]}qw/.h_nJ Xapou cets krht ele_ r_ra/;
map{y/X_/\n /;print}map{pop@$_}@/for@/