Re: PARI/GP wrong on 7!!+7+1 being not prime, or oeis.org/Mathematica?

[Karim Belabas <Karim.Belabas@u-bordeaux.fr>]

* Denis Simon [2025-09-09 14:03]:
> Another possibility is to add Bill's suggestion in the help ?factorial
> or ??factorial.

Rather to ?? _!

factorial(x) is a legacy transcendental function in PARI/GP, functionally
equivalent to gamma(x + 1) with restricted inputs (factorial is limited
to non-negative integers) :

? gamma(5)
%1 = 24.000000000000000000000000000000000000
? factorial(4)
%2 = 24.000000000000000000000000000000000000
? gamma(Pi)
%3 = 2.2880377953400324179595889090602339229
? factorial(Pi-1)
  ***   at top-level: factorial(Pi-1)
  ***                           ^-----
  ***   incorrect type in gtos [integer expected] (t_REAL).

So the natural extension of factorial() would rather be

  double_factorial(z) = sqrt(2/Pi)*2^(z/2)*gamma(z/2+1)

? double_factorial(7)
%4 = 105.00000000000000000000000000000000000

In my opinion, the current factorial() function is quite useless and
could be replaced by a generalized factorial.

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/

> ----- Mail original -----
> > De: "Bill Allombert" <Bill.Allombert@math.u-bordeaux.fr>
> > À: "pari-users" <pari-users@pari.math.u-bordeaux.fr>
> > Envoyé: Mardi 9 Septembre 2025 10:09:45
> > Objet: Re: PARI/GP wrong on 7!!+7+1 being not prime, or oeis.org/Mathematica?
> 
> > On Thu, Aug 28, 2025 at 02:37:46PM -0400, James Cloos wrote:
> >> BA> Adding this one would break backward compatibility, so this is not
> >> BA> an option.
> >> 
> >> does pari have a generalized factorial function?
> >> 
> >> maxima's doc notes this:
> >> 
> >>      For an argument ‘n’ which is not an integer, float, or rational,
> >>      ‘n!!’ yields a noun form ‘genfact (n, n/2, 2)’.
> >> 
> >> perhaps pari's docs could suggest using such a generalized factorial
> >> when a double factorial is desired?
> > 
> > You can use
> > doublefact(n)=vecprod([2*i-1|i<-[1..(n+1)\2]])
> > 
> > Cheers,
> %6 = 105.00000000000000000000000000000000000
> Bill