Re: binomial() challenge (+thanks)

["Ruud H.G. van Tol" <rvtol@isolution.nl>]

On 2022-11-28 14:42, Bill Allombert wrote:
> On Mon, Nov 28, 2022 at 12:58:57PM +0100, Ruud H.G. van Tol wrote:
> > On 2022-11-26 19:09, Bill Allombert wrote:
> > > On Sat, Nov 26, 2022 at 01:17:28PM +0100, Ruud H.G. van Tol wrote:

> > > > PARI binomial() challenge:
> > > >
> > > > binom_1(n,k) = gamma(n+1) / gamma(k+1) / gamma(n-k+1)
> > > > binom_2(n,k) = exp(lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1))
> > > Ah! You had a one-time-in-a-life opportunity to use the factorial() function
> > > and you blew it!
> > Heheh, I use factorials all the time,
> > for stuff like this:
> > https://en.wikipedia.org/wiki/Factorial_number_system
>
> I am not speaking about the factorials, but about the little-known factorial()
> GP function!

ACK, thanks for insisting.

(gp) ? factorial(2^55)
%7 = 5.319626231150072416775431621676913199 E580869065836896632


(Perl)
sub log_factorial {
    my $n= shift;

    if ( !defined $n or $n < 0 ) {
        return undef;
    } elsif ( $n <= 1 ) {  # 0, 1
        return 0;
    } elsif ( $n <= 127 ) {
        state @log_factorial;
        if ( !@log_factorial ) {  # build at first call
            @log_factorial= ( 0, 0 );
            my $f= 1;
            for my $i ( 2 .. 127 ) {
                $f *= $i;
                $log_factorial[ $i ]= log( $f );
            }
        }
        return $log_factorial[ $n ];
    }

    my $x= $n + 1;
    return ( $x - 0.5 ) * log( $x ) - $x + log( 2 * PI ) / 2
             + 1 / (   12 * ( $x ** 1 ) )
             - 1 / (  360 * ( $x ** 3 ) )
             + 1 / ( 1260 * ( $x ** 5 ) )
    ;
}
say exp(log_factorial(170));
7.25741561530806e+306


(gp) ? factorial(170)
%12 = 7.2574156153079989673967282111292631147 E306

- - - - - - - - -

That simplistic perl-version needs Math::BigFloat for n>170,
and is of course much slower.

- - - - - - - - -

$ perl -Mv5.22 -wE'
  use Math::Pari qw( factorial  setprecision );
  setprecision(40);  # weird issue: fails for 38
  say factorial(2**55);
'
5.31962623115007241677543162167691319...E580869065836896632

:)

-- Ruud