sub PTn { my @row; for( 1 .. shift ) { push @row, 1; $row [$_] += $row [$_ - 1] for reverse 1 .. @row - 2; } return @row; } sub expN { my( $P, $N, $X ) = @_; my $xDIVn = $X / $N; my @coefs = PTn( $P+1 ); my $rv = 0; my $toggle = -1; for my $p ( reverse 1 .. $P ) { $rv += $toggle * $coefs[ $p ] * $N / $p * exp( -$p * $xDIVn ); $toggle *= -1; } return $rv + $X; }