I ran your script and got 0.00693150741379456. I made a few minor changes and got 0.0069140819065811. Here's what I tried:

#!/usr/bin/perl 

use strict;
use warnings;

sub logfact {
   return gammln(shift(@_) + 1);
}

sub hypergeom {
   my ($n, $m, $N, $i) = @_;
   my $loghyp1 = logfact($m)+logfact($n)+logfact($N)+logfact($m+$n-$N)
+;
   my $loghyp2 = logfact($i)+logfact($n-$i)+logfact($m+$i-$N)+logfact(
+$N-$i)+logfact($m+$n);
   return exp($loghyp1 - $loghyp2);
}

sub gammln {
  my $xx = shift;
  my @cof = (76.180091729471457, -86.505320329416776,
             24.014098240830911, -1.231739572450155,
             0.12086509738661, -5.395239384953e-06);
  my $y = my $x = $xx;
  my $tmp = $x + 5.5;
  $tmp -= ($x + 0.5) * log($tmp);
  my $ser = 1.0000000001900149;
  foreach my $j (0 .. 5) {
     $ser += $cof[$j] / ++$y;
  }
  -$tmp + log(2.5066282746310007*$ser/$x);
}

print hypergeom(300,700,100,40),"\n";
[download]

I think that you'll need Math::GSL for Math::GSL::Randist. One other module that I like for getting gamma and log_gamma is Math::GammaFunction.