On Wed, 4 Sept 2024 at 19:21, Tom Lane <t...@sss.pgh.pa.us> wrote:
>
> >> I'm not sure why you are doing the testing for special values (NaN etc.)
> >> yourself when the C library function already does it. For example, if I
> >> remove the isnan(arg1) check in your dlgamma(), then it still behaves
> >> the same in all tests.
>
> > It's useful to do that so that we don't need to assume that every
> > platform conforms to the POSIX standard, and because it simplifies the
> > later overflow checks. This is consistent with the approach taken in
> > other functions, such as dexp(), dsin(), dcos(), etc.
>
> dexp() and those other functions date from the late stone age, before
> it was safe to assume that libm's behavior matched the POSIX specs.
> Today I think we can assume that, at least till proven differently.
> There's not necessarily anything wrong with what you've coded, but
> I don't buy this argument for it.
>
OK, thinking about this some more, I think we should reserve overflow
errors for genuine overflows, which I take to mean cases where the
exact mathematical result should be finite, but is too large to be
represented in a double.
In particular, this means that zero and negative integer inputs are
not genuine overflows, but should return NaN or +/-Inf, as described
in the POSIX spec.
Doing that, and assuming that tgamma() and lgamma() behave according
to spec, leads to the attached, somewhat simpler patch.
Regards,
Dean
diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
new file mode 100644
index 461fc3f..b2a3368
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -1399,6 +1399,27 @@ SELECT NOT(ROW(table.*) IS NOT NULL) FRO
<row>
<entry role="func_table_entry"><para role="func_signature">
<indexterm>
+ <primary>gamma</primary>
+ </indexterm>
+ <function>gamma</function> ( <type>double precision</type> )
+ <returnvalue>double precision</returnvalue>
+ </para>
+ <para>
+ Gamma function
+ </para>
+ <para>
+ <literal>gamma(0.5)</literal>
+ <returnvalue>1.772453850905516</returnvalue>
+ </para>
+ <para>
+ <literal>gamma(6)</literal>
+ <returnvalue>120</returnvalue>
+ </para></entry>
+ </row>
+
+ <row>
+ <entry role="func_table_entry"><para role="func_signature">
+ <indexterm>
<primary>gcd</primary>
</indexterm>
<function>gcd</function> ( <replaceable>numeric_type</replaceable>, <replaceable>numeric_type</replaceable> )
@@ -1436,6 +1457,23 @@ SELECT NOT(ROW(table.*) IS NOT NULL) FRO
</para></entry>
</row>
+ <row>
+ <entry role="func_table_entry"><para role="func_signature">
+ <indexterm>
+ <primary>lgamma</primary>
+ </indexterm>
+ <function>lgamma</function> ( <type>double precision</type> )
+ <returnvalue>double precision</returnvalue>
+ </para>
+ <para>
+ Natural logarithm of the absolute value of the gamma function
+ </para>
+ <para>
+ <literal>lgamma(1000)</literal>
+ <returnvalue>5905.220423209181</returnvalue>
+ </para></entry>
+ </row>
+
<row>
<entry role="func_table_entry"><para role="func_signature">
<indexterm>
diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
new file mode 100644
index f709c21..00da60b
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2787,6 +2787,63 @@ derfc(PG_FUNCTION_ARGS)
}
+/* ========== GAMMA FUNCTIONS ========== */
+
+
+/*
+ * dgamma - returns the gamma function of arg1
+ */
+Datum
+dgamma(PG_FUNCTION_ARGS)
+{
+ float8 arg1 = PG_GETARG_FLOAT8(0);
+ float8 result;
+
+ /* note: the POSIX/C99 gamma function is called "tgamma", not "gamma" */
+ errno = 0;
+ result = tgamma(arg1);
+
+ /*
+ * If an ERANGE error occurs, it means there is an overflow. This may also
+ * happen if the input is +/-0, which is not a genuine overflow, and the
+ * result should be +/-infinity.
+ */
+ if (errno == ERANGE && arg1 != 0)
+ float_overflow_error();
+
+ PG_RETURN_FLOAT8(result);
+}
+
+
+/*
+ * dlgamma - natural logarithm of absolute value of gamma of arg1
+ */
+Datum
+dlgamma(PG_FUNCTION_ARGS)
+{
+ float8 arg1 = PG_GETARG_FLOAT8(0);
+ float8 result;
+
+ /*
+ * Note: lgamma may not be thread-safe because it may write to a global
+ * variable signgam, which may not be thread-local. However, this doesn't
+ * matter to us, since we don't use signgam.
+ */
+ errno = 0;
+ result = lgamma(arg1);
+
+ /*
+ * If an ERANGE error occurs, it means there is an overflow. This also
+ * happens if the input is zero or a negative integer, which are not
+ * genuine overflows, and the result should be infinity.
+ */
+ if (errno == ERANGE && !(arg1 <= 0 && floor(arg1) == arg1))
+ float_overflow_error();
+
+ PG_RETURN_FLOAT8(result);
+}
+
+
/*
* =========================
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
new file mode 100644
index ff5436a..cee8398
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -3490,6 +3490,13 @@
proname => 'erfc', prorettype => 'float8', proargtypes => 'float8',
prosrc => 'derfc' },
+{ oid => '8702', descr => 'gamma function',
+ proname => 'gamma', prorettype => 'float8', proargtypes => 'float8',
+ prosrc => 'dgamma' },
+{ oid => '8703', descr => 'natural logarithm of absolute value of gamma function',
+ proname => 'lgamma', prorettype => 'float8', proargtypes => 'float8',
+ prosrc => 'dlgamma' },
+
{ oid => '1618',
proname => 'interval_mul', prorettype => 'interval',
proargtypes => 'interval float8', prosrc => 'interval_mul' },
diff --git a/src/test/regress/expected/float8.out b/src/test/regress/expected/float8.out
new file mode 100644
index 344d6b7..d7052e3
--- a/src/test/regress/expected/float8.out
+++ b/src/test/regress/expected/float8.out
@@ -830,6 +830,128 @@ FROM (VALUES (float8 '-infinity'),
(22 rows)
RESET extra_float_digits;
+-- gamma functions
+-- we run these with extra_float_digits = -1, to get consistently rounded
+-- results on all platforms.
+SET extra_float_digits = -1;
+SELECT x,
+ gamma(x),
+ lgamma(x)
+FROM (VALUES (0.5), (1), (2), (3), (4), (5)) AS t(x);
+ x | gamma | lgamma
+-----+-----------------+------------------
+ 0.5 | 1.7724538509055 | 0.5723649429247
+ 1 | 1 | 0
+ 2 | 1 | 0
+ 3 | 2 | 0.69314718055995
+ 4 | 6 | 1.7917594692281
+ 5 | 24 | 3.1780538303479
+(6 rows)
+
+-- test special cases for gamma functions
+SELECT gamma(float8 '-infinity');
+ gamma
+-------
+ NaN
+(1 row)
+
+SELECT lgamma(float8 '-infinity');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '-1000000.5');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '-1000000.5');
+ lgamma
+------------------
+ -12815524.147684
+(1 row)
+
+SELECT gamma(float8 '-1000000');
+ gamma
+-------
+ NaN
+(1 row)
+
+SELECT lgamma(float8 '-1000000');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '-1');
+ gamma
+-------
+ NaN
+(1 row)
+
+SELECT lgamma(float8 '-1');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '-0');
+ gamma
+-----------
+ -Infinity
+(1 row)
+
+SELECT lgamma(float8 '-0');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '0');
+ gamma
+----------
+ Infinity
+(1 row)
+
+SELECT lgamma(float8 '0');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '1000000');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '1000000');
+ lgamma
+-----------------
+ 12815504.569148
+(1 row)
+
+SELECT lgamma(float8 '1e308');
+ERROR: value out of range: overflow
+SELECT gamma(float8 'infinity');
+ gamma
+----------
+ Infinity
+(1 row)
+
+SELECT lgamma(float8 'infinity');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 'nan');
+ gamma
+-------
+ NaN
+(1 row)
+
+SELECT lgamma(float8 'nan');
+ lgamma
+--------
+ NaN
+(1 row)
+
+RESET extra_float_digits;
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
ERROR: "10e400" is out of range for type double precision
diff --git a/src/test/regress/sql/float8.sql b/src/test/regress/sql/float8.sql
new file mode 100644
index 98e9926..36fd7b8
--- a/src/test/regress/sql/float8.sql
+++ b/src/test/regress/sql/float8.sql
@@ -245,6 +245,36 @@ FROM (VALUES (float8 '-infinity'),
RESET extra_float_digits;
+-- gamma functions
+-- we run these with extra_float_digits = -1, to get consistently rounded
+-- results on all platforms.
+SET extra_float_digits = -1;
+SELECT x,
+ gamma(x),
+ lgamma(x)
+FROM (VALUES (0.5), (1), (2), (3), (4), (5)) AS t(x);
+-- test special cases for gamma functions
+SELECT gamma(float8 '-infinity');
+SELECT lgamma(float8 '-infinity');
+SELECT gamma(float8 '-1000000.5');
+SELECT lgamma(float8 '-1000000.5');
+SELECT gamma(float8 '-1000000');
+SELECT lgamma(float8 '-1000000');
+SELECT gamma(float8 '-1');
+SELECT lgamma(float8 '-1');
+SELECT gamma(float8 '-0');
+SELECT lgamma(float8 '-0');
+SELECT gamma(float8 '0');
+SELECT lgamma(float8 '0');
+SELECT gamma(float8 '1000000');
+SELECT lgamma(float8 '1000000');
+SELECT lgamma(float8 '1e308');
+SELECT gamma(float8 'infinity');
+SELECT lgamma(float8 'infinity');
+SELECT gamma(float8 'nan');
+SELECT lgamma(float8 'nan');
+RESET extra_float_digits;
+
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
