Re: gamma() and lgamma() functions

Hi Dean,

Thanks for the patch!

I reviewed the patch and it looks good to me. I’ve run some tests, and
everything works as expected.

I have one minor thought:

On Sun, Mar 2, 2025 at 2:30 AM Dean Rasheed <dean(dot)a(dot)rasheed(at)gmail(dot)com>
wrote:

> On Thu, 14 Nov 2024 at 22:35, Dean Rasheed <dean(dot)a(dot)rasheed(at)gmail(dot)com>
> wrote:
> >
> > On Thu, 14 Nov 2024 at 16:28, Dmitry Koval <d(dot)koval(at)postgrespro(dot)ru>
> wrote:
> > >
> > > >SELECT gamma(float8 '1e-320');
> > > ERROR: value out of range: overflow
> > >
> > > >SELECT gamma(float8 '0');
> > > gamma
> > > ----------
> > > Infinity
> > > (1 row)
> > >
> > > Perhaps it would be logical if the behavior in these cases was the same
> > > (either ERROR or 'Infinity')?
> >
> > In general, I think that having gamma() throw overflow errors is
> > useful for spotting real problems in user code.
> >
>
> Thinking about this afresh, I remain of the opinion that having the
> gamma function throw overflow errors is useful for inputs that are too
> large, like gamma(200). But then, if we're going to do that, maybe we
> should just do the same for other invalid inputs (zero, negative
> integers, and -Inf). That resolves the consistency issue with inputs
> very close to zero, and seems like a practical solution.
>

Raising an error instead of silently returning Inf for invalid inputs
like zero, negative integers, and -Inf makes sense to me.

I also wondered whether we should distinguish domain errors
(e.g., zero or negative integers) from range errors (e.g., overflow).
I tested tgamma() and lgamma() on Ubuntu 24.04, macOS 15.3.2, and
Windows 11 (code attached).

Here’s a summary of return values, errnos, and exceptions on these
platforms: (errno=33 is EDOM and errno=34 is ERANGE)

tgamma(0):
Ubuntu : inf, errno=34, FE_DIVBYZERO
macOS : inf, errno=0, FE_DIVBYZERO
Windows: inf, errno=34, FE_DIVBYZERO

lgamma(0):
Ubuntu : inf, errno=34, FE_DIVBYZERO
macOS : inf, errno=0
Windows: inf, errno=34, FE_DIVBYZERO

tgamma(-1):
Ubuntu : nan, errno=33, FE_INVALID
macOS : nan, errno=0, FE_INVALID
Windows: nan, errno=33, FE_INVALID

lgamma(-1):
Ubuntu : inf, errno=34, FE_DIVBYZERO
macOS : inf, errno=0, FE_DIVBYZERO
Windows: inf, errno=34, FE_DIVBYZERO

tgamma(-1000.5):
Ubuntu : -0.0, errno=34, FE_UNDERFLOW
macOS : -0.0, errno=0
Windows: 0.0, errno=34

lgamma(-1000.5):
Ubuntu : -5914.44, errno=0
macOS : -5914.44, errno=0
Windows: -5914.44, errno=0

tgamma(1000):
Ubuntu : inf, errno=34, FE_OVERFLOW
macOS : inf, errno=0, FE_OVERFLOW
Windows: inf, errno=34, FE_OVERFLOW

lgamma(1000):
Ubuntu : 5905.22, errno=0
macOS : 5905.22, errno=0
Windows: 5905.22, errno=0

tgamma(+inf):
Ubuntu : inf, errno=0
macOS : inf, errno=0
Windows: inf, errno=0

lgamma(+inf):
Ubuntu : inf, errno=0
macOS : inf, errno=0
Windows: inf, errno=0

tgamma(-inf):
Ubuntu : nan, errno=33, FE_INVALID
macOS : nan, errno=0, FE_INVALID
Windows: nan, errno=33, FE_INVALID

lgamma(-inf):
Ubuntu : inf, errno=0
macOS : inf, errno=0
Windows: inf, errno=0

In the current implementation, float_overflow_error() is used for both
domain and range (overflow) errors. If we did want to distinguish
them, maybe we could use float_zero_divide_error() for domain errors?
Not sure it adds much value, though - just a thought, and I’m fine
either way.

OTOH I wondered if we could leverage fetestexcept(), but then I saw
this comment in dcos():

* check for errors. POSIX allows an error to be reported either via
> * errno or via fetestexcept(), but currently we only support checking
> * errno. (fetestexcept() is rumored to report underflow unreasonably
> * early on some platforms, so it's not clear that believing it would be a
> * net improvement anyway.)

So again, I’m totally fine with keeping the current semantics as-is.
That said, it might be helpful to update the comment in dlgamma():
+ /*
+ * If an ERANGE error occurs, it means there is an overflow.
+ *
to clarify that an ERANGE error can also indicate a pole error (e.g.,
input 0 or a negative integer), not just overflow.

Thanks,
Alex