gamma() and lgamma() functions
Attached is a patch adding support for the gamma and log-gamma
functions. See, for example:
https://en.wikipedia.org/wiki/Gamma_function
I think these are very useful general-purpose mathematical functions.
They're part of C99, and they're commonly included in other
mathematical libraries, such as the python math module, so I think
it's useful to make them available from SQL.
The error-handling for these functions seems to be a little trickier
than most, so that might need further discussion.
Regards,
Dean
Attachments:
gamma-and-lgamma.patchtext/x-patch; charset=US-ASCII; name=gamma-and-lgamma.patchDownload
diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
new file mode 100644
index 5a16910..6464a8b
--- 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 cbbb8ae..e43aa42
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2779,6 +2779,91 @@ 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;
+
+ /*
+ * Per the POSIX spec, return NaN if the input is NaN, +/-infinity if the
+ * input is +/-0, NaN if the input is -infinity, and infinity if the input
+ * is infinity.
+ */
+ if (isnan(arg1))
+ PG_RETURN_FLOAT8(get_float8_nan());
+ if (arg1 == 0)
+ {
+ if (signbit(arg1))
+ PG_RETURN_FLOAT8(-get_float8_infinity());
+ else
+ PG_RETURN_FLOAT8(get_float8_infinity());
+ }
+ if (isinf(arg1))
+ {
+ if (arg1 < 0)
+ PG_RETURN_FLOAT8(get_float8_nan());
+ else
+ PG_RETURN_FLOAT8(get_float8_infinity());
+ }
+
+ /*
+ * Note that the POSIX/C99 gamma function is called "tgamma", not "gamma".
+ *
+ * For negative integer inputs, it may raise a domain error or a pole
+ * error, or return NaN or +/-infinity (the POSIX spec indicates that this
+ * may change in future implementations). Since we can't reliably detect
+ * integer inputs above a certain size, and we can't distinguish this from
+ * any other overflow error, just treat them all the same.
+ */
+ errno = 0;
+ result = tgamma(arg1);
+ if (errno != 0 || isinf(result) || isnan(result))
+ 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;
+
+ /* Per the POSIX spec, return NaN if the input is NaN */
+ if (isnan(arg1))
+ PG_RETURN_FLOAT8(get_float8_nan());
+
+ errno = 0;
+ result = lgamma(arg1);
+
+ /*
+ * If an ERANGE error occurs, it means there is an overflow. This can
+ * happen if the input is 0, a negative integer, -infinity, or infinity.
+ * In each of those cases, the correct result is infinity. However, it
+ * can also happen in other cases where the correct result is finite, but
+ * too large to be represented (e.g., if the input is larger than around
+ * 2.5e305). There seems to be no way to distinguish those cases, so we
+ * just return infinity for them too. That's not ideal, but in practice,
+ * lgamma() is much less likely to overflow than tgamma(), so it seems OK.
+ */
+ if (unlikely(errno == ERANGE))
+ result = get_float8_infinity();
+
+ 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 6a5476d..acaf4fb
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -3484,6 +3484,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..f16d670
--- a/src/test/regress/expected/float8.out
+++ b/src/test/regress/expected/float8.out
@@ -830,6 +830,50 @@ 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 (float8 '-infinity'),
+ ('-0'), (0), (0.5), (1), (2), (3), (4), (5),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+ x | gamma | lgamma
+-----------+-----------------+------------------
+ -Infinity | NaN | Infinity
+ -0 | -Infinity | Infinity
+ 0 | Infinity | Infinity
+ 0.5 | 1.7724538509055 | 0.5723649429247
+ 1 | 1 | 0
+ 2 | 1 | 0
+ 3 | 2 | 0.69314718055995
+ 4 | 6 | 1.7917594692281
+ 5 | 24 | 3.1780538303479
+ Infinity | Infinity | Infinity
+ NaN | NaN | NaN
+(11 rows)
+
+-- invalid inputs for gamma()
+SELECT gamma(-1);
+ERROR: value out of range: overflow
+SELECT gamma(1000000);
+ERROR: value out of range: overflow
+-- valid for lgamma()
+SELECT lgamma(-1);
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT lgamma(1000000);
+ lgamma
+-----------------
+ 12815504.569148
+(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..68c68f0
--- a/src/test/regress/sql/float8.sql
+++ b/src/test/regress/sql/float8.sql
@@ -245,6 +245,25 @@ 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 (float8 '-infinity'),
+ ('-0'), (0), (0.5), (1), (2), (3), (4), (5),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+-- invalid inputs for gamma()
+SELECT gamma(-1);
+SELECT gamma(1000000);
+-- valid for lgamma()
+SELECT lgamma(-1);
+SELECT lgamma(1000000);
+
+RESET extra_float_digits;
+
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
On Mon, Jul 1, 2024 at 5:33 PM Dean Rasheed <dean.a.rasheed@gmail.com>
wrote:
Attached is a patch adding support for the gamma and log-gamma
functions. See, for example:https://en.wikipedia.org/wiki/Gamma_function
I think these are very useful general-purpose mathematical functions.
They're part of C99, and they're commonly included in other
mathematical libraries, such as the python math module, so I think
it's useful to make them available from SQL.The error-handling for these functions seems to be a little trickier
than most, so that might need further discussion.Regards,
Dean
Hi! The patch file seems broken.
git apply gamma-and-lgamma.patch
error: git apply: bad git-diff — exptec /dev/null in line 2
Best regards, Stepan Neretin.
On 1 Jul 2024, at 16:20, Stepan Neretin <sncfmgg@gmail.com> wrote:
The patch file seems broken.
git apply gamma-and-lgamma.patch error: git apply: bad git-diff — exptec /dev/null in line 2
It's a plain patch file, if you apply it with patch and not git it will work fine:
$ patch -p1 < gamma-and-lgamma.patch
patching file 'doc/src/sgml/func.sgml'
patching file 'src/backend/utils/adt/float.c'
patching file 'src/include/catalog/pg_proc.dat'
patching file 'src/test/regress/expected/float8.out'
patching file 'src/test/regress/sql/float8.sql'
--
Daniel Gustafsson
On Mon, Jul 1, 2024 at 5:33 PM Dean Rasheed <dean.a.rasheed@gmail.com>
wrote:
Attached is a patch adding support for the gamma and log-gamma
functions. See, for example:https://en.wikipedia.org/wiki/Gamma_function
I think these are very useful general-purpose mathematical functions.
They're part of C99, and they're commonly included in other
mathematical libraries, such as the python math module, so I think
it's useful to make them available from SQL.The error-handling for these functions seems to be a little trickier
than most, so that might need further discussion.Regards,
Dean
I tried to review the patch without applying it. It looks good to me, but I
have one notice:
ERROR: value out of range: overflow. I think we need to add information
about the available ranges in the error message
On 2024-Jul-01, Stepan Neretin wrote:
I have one notice:
ERROR: value out of range: overflow. I think we need to add information
about the available ranges in the error message
I think this is a project of its own. The error comes from calling
float_overflow_error(), which is a generic routine used in several
functions which have different overflow conditions. It's not clear to
me that throwing errors listing the specific range for each case is
worth the effort.
--
Álvaro Herrera Breisgau, Deutschland — https://www.EnterpriseDB.com/
"Escucha y olvidarás; ve y recordarás; haz y entenderás" (Confucio)
On Mon, 1 Jul 2024 at 16:04, Alvaro Herrera <alvherre@alvh.no-ip.org> wrote:
On 2024-Jul-01, Stepan Neretin wrote:
I have one notice:
ERROR: value out of range: overflow. I think we need to add information
about the available ranges in the error messageI think this is a project of its own. The error comes from calling
float_overflow_error(), which is a generic routine used in several
functions which have different overflow conditions. It's not clear to
me that throwing errors listing the specific range for each case is
worth the effort.
Right. It's also worth noting that gamma() has several distinct ranges
of validity for which it doesn't overflow, so it'd be hard to codify
that succinctly in an error message.
Something that bothers me about float.c is that there is precious
little consistency as to whether functions raise overflow errors or
return infinity. For example, exp() gives an overflow error for
sufficiently large (finite) inputs, whereas sinh() and cosh() (very
closely related) return infinity. I think raising an error is the more
technically correct thing to do, but returning infinity is sometimes
perhaps more practically useful.
However, given how much more quickly the result from gamma() explodes,
I think it's more useful for it to raise an error so that you know you
have a problem, and that you probably want to use lgamma() instead.
(As an aside, I've seen people (and ChatGPT) write algorithms to solve
the Birthday Problem using the gamma() function. That doesn't work
because gamma(366) overflows, so you have to use lgamma().)
Regards,
Dean
On 01.07.24 12:33, Dean Rasheed wrote:
Attached is a patch adding support for the gamma and log-gamma
functions. See, for example:https://en.wikipedia.org/wiki/Gamma_function
I think these are very useful general-purpose mathematical functions.
They're part of C99, and they're commonly included in other
mathematical libraries, such as the python math module, so I think
it's useful to make them available from SQL.
What are examples of where this would be useful in a database context?
The error-handling for these functions seems to be a little trickier
than most, so that might need further discussion.
Yeah, this is quite something.
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. However, the same does not happen in your
dgamma(). The overflow checks after the C library call are written
differently for the two functions. dgamma() does not check errno for
ERANGE for example. It might also be good if dgamma() checked errno for
EDOM, because other the result of gamma(-1) being "overflow" seems a bit
wrong.
I'm also not clear why you are turning a genuine result overflow into
infinity in lgamma().
I think this could use some kind of chart or something about all the
possible behaviors and how the C library reports them and what we intend
to do with them.
Btw., I'm reading that lgamma() in the C library is not necessarily
thread-safe. Is that a possible problem?
On Fri, 23 Aug 2024 at 13:40, Peter Eisentraut <peter@eisentraut.org> wrote:
What are examples of where this would be useful in a database context?
gamma() and lgamma() are the kinds of functions that are generally
useful for a variety of tasks like statistical analysis and
combinatorial computations, and having them in the database allows
those sort of computations to be performed without the need to export
the data to an external tool. We discussed adding them in a thread
last year [1]/messages/by-id/CA+hUKGKJAcB8Q5qziKTTSnkA4Mnv_6f+7-_XUgbh9jFjSdEFQg@mail.gmail.com, and there has been at least a little prior interest
[2]: https://stackoverflow.com/questions/58884066/how-can-i-run-the-equivalent-of-excels-gammaln-function-in-postgres
[1]: /messages/by-id/CA+hUKGKJAcB8Q5qziKTTSnkA4Mnv_6f+7-_XUgbh9jFjSdEFQg@mail.gmail.com
[2]: https://stackoverflow.com/questions/58884066/how-can-i-run-the-equivalent-of-excels-gammaln-function-in-postgres
Of course, there's a somewhat fuzzy line between what is generally
useful enough, and what is too specialised for core Postgres, but I
would argue that these qualify, since they are part of C99, and
commonly appear in other general-purpose math libraries like the
Python math module.
The error-handling for these functions seems to be a little trickier
than most, so that might need further discussion.Yeah, this is quite something.
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.
The overflow checks after the C library call are written
differently for the two functions. dgamma() does not check errno for
ERANGE for example. It might also be good if dgamma() checked errno for
EDOM, because other the result of gamma(-1) being "overflow" seems a bit
wrong.
They're intentionally different because the functions themselves are
different. In this case:
select gamma(-1);
ERROR: value out of range: overflow
it is correct to throw an error, because gamma(-1) is undefined (it
goes to -Inf as x goes to -1 from above, and +Inf as x goes to -1 from
below, so there is no well-defined limit).
I've updated the patch to give a more specific error message for
negative integer inputs, as opposed to other overflow cases.
Relying on errno being ERANGE or EDOM doesn't seem possible though,
because the docs indicate that, while its behaviour is one thing
today, per POSIX, that will change in the future.
By contrast, lgamma() does not raise an error for such inputs:
select lgamma(-1);
lgamma
----------
Infinity
This is correct because lgamma() is the log of the absolute value of
the gamma function, so the limit is +Inf as x goes to -1 from both
sides.
I'm also not clear why you are turning a genuine result overflow into
infinity in lgamma().
OK, I've changed that to only return infinity if the input is
infinite, zero, or a negative integer. Otherwise, it now throws an
overflow error.
Btw., I'm reading that lgamma() in the C library is not necessarily
thread-safe. Is that a possible problem?
It's not quite clear what to do about that. Some platforms may define
the lgamma_r() function, but not all. Do we have a standard pattern
for dealing with functions for which there is no thread-safe
alternative?
Regards,
Dean
Attachments:
v2-gamma-and-lgamma.patchtext/x-patch; charset=US-ASCII; name=v2-gamma-and-lgamma.patchDownload
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..4694ea2
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2787,6 +2787,97 @@ 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;
+
+ /*
+ * Per the POSIX spec, return NaN if the input is NaN, +/-infinity if the
+ * input is +/-0, NaN if the input is -infinity, and infinity if the input
+ * is infinity.
+ */
+ if (isnan(arg1))
+ PG_RETURN_FLOAT8(get_float8_nan());
+ if (arg1 == 0)
+ {
+ if (signbit(arg1))
+ PG_RETURN_FLOAT8(-get_float8_infinity());
+ else
+ PG_RETURN_FLOAT8(get_float8_infinity());
+ }
+ if (isinf(arg1))
+ {
+ if (arg1 < 0)
+ PG_RETURN_FLOAT8(get_float8_nan());
+ else
+ PG_RETURN_FLOAT8(get_float8_infinity());
+ }
+
+ /*
+ * Note that the POSIX/C99 gamma function is called "tgamma", not "gamma".
+ *
+ * For negative integer inputs, it may raise a domain error or a pole
+ * error, or return NaN or +/-infinity (the POSIX spec indicates that this
+ * may change in future implementations). Try to distinguish these cases
+ * from other overflow cases, and give a more specific error message.
+ */
+ errno = 0;
+ result = tgamma(arg1);
+ if (errno != 0 || isinf(result) || isnan(result))
+ {
+ if (arg1 < 0 && floor(arg1) == arg1)
+ ereport(ERROR,
+ errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+ errmsg("the gamma function is undefined for negative integer inputs"));
+ 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;
+
+ /* Per the POSIX spec, return NaN if the input is NaN */
+ if (isnan(arg1))
+ PG_RETURN_FLOAT8(get_float8_nan());
+
+ errno = 0;
+ result = lgamma(arg1);
+
+ /*
+ * If an ERANGE error occurs, it means there is an overflow. This can
+ * happen if the input is 0, a negative integer, -infinity, or infinity.
+ * In each of those cases, the correct result is infinity. Otherwise, it
+ * is a genuine overflow error.
+ */
+ if (unlikely(errno == ERANGE))
+ {
+ if (isinf(arg1) || (arg1 <= 0 && floor(arg1) == arg1))
+ result = get_float8_infinity();
+ else
+ 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..a31624b
--- a/src/test/regress/expected/float8.out
+++ b/src/test/regress/expected/float8.out
@@ -830,6 +830,53 @@ 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 (float8 '-infinity'),
+ ('-0'), (0), (0.5), (1), (2), (3), (4), (5),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+ x | gamma | lgamma
+-----------+-----------------+------------------
+ -Infinity | NaN | Infinity
+ -0 | -Infinity | Infinity
+ 0 | Infinity | Infinity
+ 0.5 | 1.7724538509055 | 0.5723649429247
+ 1 | 1 | 0
+ 2 | 1 | 0
+ 3 | 2 | 0.69314718055995
+ 4 | 6 | 1.7917594692281
+ 5 | 24 | 3.1780538303479
+ Infinity | Infinity | Infinity
+ NaN | NaN | NaN
+(11 rows)
+
+-- invalid inputs for gamma()
+SELECT gamma(-1);
+ERROR: the gamma function is undefined for negative integer inputs
+SELECT gamma(1000000);
+ERROR: value out of range: overflow
+-- valid for lgamma()
+SELECT lgamma(-1);
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT lgamma(1000000);
+ lgamma
+-----------------
+ 12815504.569148
+(1 row)
+
+-- inavlid input for lgamma()
+SELECT lgamma(1e308);
+ERROR: value out of range: overflow
+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..de2656a
--- a/src/test/regress/sql/float8.sql
+++ b/src/test/regress/sql/float8.sql
@@ -245,6 +245,27 @@ 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 (float8 '-infinity'),
+ ('-0'), (0), (0.5), (1), (2), (3), (4), (5),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+-- invalid inputs for gamma()
+SELECT gamma(-1);
+SELECT gamma(1000000);
+-- valid for lgamma()
+SELECT lgamma(-1);
+SELECT lgamma(1000000);
+-- inavlid input for lgamma()
+SELECT lgamma(1e308);
+
+RESET extra_float_digits;
+
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
Dean Rasheed <dean.a.rasheed@gmail.com> writes:
On Fri, 23 Aug 2024 at 13:40, Peter Eisentraut <peter@eisentraut.org> wrote:
What are examples of where this would be useful in a database context?
Of course, there's a somewhat fuzzy line between what is generally
useful enough, and what is too specialised for core Postgres, but I
would argue that these qualify, since they are part of C99, and
commonly appear in other general-purpose math libraries like the
Python math module.
Yeah, I think any math function that's part of C99 or POSIX is
arguably of general interest.
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.
Btw., I'm reading that lgamma() in the C library is not necessarily
thread-safe. Is that a possible problem?
It's not quite clear what to do about that.
Per the Linux man page, the reason lgamma() isn't thread-safe is
The lgamma(), lgammaf(), and lgammal() functions return the natural
logarithm of the absolute value of the Gamma function. The sign of the
Gamma function is returned in the external integer signgam declared in
<math.h>. It is 1 when the Gamma function is positive or zero, -1 when
it is negative.
AFAICS this patch doesn't inspect signgam, so whether it gets
overwritten by a concurrent thread wouldn't matter. However,
it'd be a good idea to add a comment noting the hazard.
(Presumably, the reason POSIX says "need not be thread-safe"
is that an implementation could make it thread-safe by
making signgam thread-local, but the standard doesn't wish
to mandate that.)
regards, tom lane
I wrote:
AFAICS this patch doesn't inspect signgam, so whether it gets
overwritten by a concurrent thread wouldn't matter. However,
it'd be a good idea to add a comment noting the hazard.
Further to that ... I looked at POSIX issue 8 (I had been reading 7)
and found this illuminating discussion:
Earlier versions of this standard did not require lgamma(),
lgammaf(), and lgammal() to be thread-safe because signgam was a
global variable. They are now required to be thread-safe to align
with the ISO C standard (which, since the introduction of threads
in 2011, requires that they avoid data races), with the exception
that they need not avoid data races when storing a value in the
signgam variable. Since signgam is not specified by the ISO C
standard, this exception is not a conflict with that standard.
So the other reason to avoid using signgam is that it might
not exist at all in some libraries.
regards, tom lane
On Wed, 4 Sept 2024 at 19:21, Tom Lane <tgl@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
Attachments:
v3-gamma-and-lgamma.patchtext/x-patch; charset=US-ASCII; name=v3-gamma-and-lgamma.patchDownload
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');
On Fri, 6 Sept 2024 at 10:42, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
... assuming that tgamma() and lgamma() behave according to spec ...
Nope, that was too much to hope for. Let's see if this fares any better.
Regards,
Dean
Attachments:
v4-gamma-and-lgamma.patchtext/x-patch; charset=US-ASCII; name=v4-gamma-and-lgamma.patchDownload
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..6e3252e
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2787,6 +2787,80 @@ 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.
+ *
+ * On some platforms, tgamma() will not set errno but just return infinity
+ * or zero to report overflow/underflow; therefore, test both cases (note
+ * that, like the exponential function, the gamma function has no zeros).
+ */
+ if (errno == ERANGE && arg1 != 0)
+ {
+ if (result != 0.0)
+ float_overflow_error();
+ else
+ float_underflow_error();
+ }
+ else if (isinf(result) && arg1 != 0 && !isinf(arg1))
+ float_overflow_error();
+ else if (result == 0.0)
+ float_underflow_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.
+ *
+ * On some platforms, lgamma() will not set errno but just return infinity
+ * to report overflow, but it should never underflow.
+ */
+ if ((errno == ERANGE || isinf(result)) && !isinf(arg1) &&
+ !(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..5c91ef1
--- 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: underflow
+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');
Hi!
I think having gamma() and lgamma() functions in PostgreSQL would be
useful for some users, this is asked [1]https://stackoverflow.com/questions/58884066/how-can-i-run-the-equivalent-of-excels-gammaln-function-in-postgres.
I have a question regarding the current implementation of gamma()
function. Code block:
+ if (errno == ERANGE && arg1 != 0)
+ {
+ if (result != 0.0)
+ float_overflow_error();
+ else
+ float_underflow_error();
+ }
+ else if (isinf(result) && arg1 != 0 && !isinf(arg1))
+ float_overflow_error();
+ else if (result == 0.0)
+ float_underflow_error();Why in some cases (if arg1 is close to 0, but not 0) an error
(float_overflow_error) will be returned, but in the case of "arg1 = 0"
the value 'Infinity' will be returned?
For example:
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')?
Links:
[1]: https://stackoverflow.com/questions/58884066/how-can-i-run-the-equivalent-of-excels-gammaln-function-in-postgres
https://stackoverflow.com/questions/58884066/how-can-i-run-the-equivalent-of-excels-gammaln-function-in-postgres
--
With best regards,
Dmitry Koval
Postgres Professional: http://postgrespro.com
On Thu, 14 Nov 2024 at 16:28, Dmitry Koval <d.koval@postgrespro.ru> wrote:
I think having gamma() and lgamma() functions in PostgreSQL would be
useful for some users, this is asked [1].
Thanks for looking.
SELECT gamma(float8 '1e-320');
ERROR: value out of range: overflow
SELECT gamma(float8 '0');
gamma
----------
Infinity
(1 row)
That's correct since gamma(1e-320) is roughly 1e320, which overflows
the double precision type, but it's not actually infinite, whereas
gamma(0) is infinite.
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. In my experience, it's
uncommon to pass negative or even close-to-zero inputs to gamma(), but
it is common to fail to appreciate just how quickly the result grows
for positive inputs (it overflows for inputs just over 171.6), so it
seems better to be made aware of such problems (which might be solved
by using lgamma() instead).
So I think throwing overflow errors is the most useful thing to do
from a practical point of view, and if we do that for too-large
inputs, it makes sense to do the same for too-close-to-zero inputs.
OTOH, gamma(+/-0) = +/-Infinity is right from a mathematical point of
view, and consistent with the POSIX spec, and it's also consistent
with the functions cot() and cotd().
Regards,
Dean
On Thu, 14 Nov 2024 at 22:35, Dean Rasheed <dean.a.rasheed@gmail.com> wrote:
On Thu, 14 Nov 2024 at 16:28, Dmitry Koval <d.koval@postgrespro.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.
Attached is an update doing that.
Regards,
Dean
Attachments:
v5-0001-Add-support-for-gamma-and-lgamma-functions.patchtext/x-patch; charset=US-ASCII; name=v5-0001-Add-support-for-gamma-and-lgamma-functions.patchDownload
From 4613ddc7267343e74703991267a1adb45eff3403 Mon Sep 17 00:00:00 2001
From: Dean Rasheed <dean.a.rasheed@gmail.com>
Date: Sun, 2 Mar 2025 09:23:20 +0000
Subject: [PATCH v5] Add support for gamma() and lgamma() functions.
These are useful general-purpose math functions which are included in
POSIX and C99, and are commonly included in other math libraries, so
expose them as SQL-callable functions.
Author: Dean Rasheed <dean.a.rasheed@gmail.com>
Reviewed-by: Stepan Neretin <sncfmgg@gmail.com>
Reviewed-by: Tom Lane <tgl@sss.pgh.pa.us>
Reviewed-by: Peter Eisentraut <peter@eisentraut.org>
Reviewed-by: Dmitry Koval <d.koval@postgrespro.ru>
Discussion: https://postgr.es/m/CAEZATCXpGyfjXCirFk9au+FvM0y2Ah+2-0WSJx7MO368ysNUPA@mail.gmail.com
---
doc/src/sgml/func.sgml | 38 ++++++++++++
src/backend/utils/adt/float.c | 87 ++++++++++++++++++++++++++++
src/include/catalog/pg_proc.dat | 7 +++
src/test/regress/expected/float8.out | 57 ++++++++++++++++++
src/test/regress/sql/float8.sql | 23 ++++++++
5 files changed, 212 insertions(+)
diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
index 0e6c534965..143031e03e 100644
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -1396,6 +1396,27 @@ SELECT NOT(ROW(table.*) IS NOT NULL) FROM TABLE; -- detect at least one null in
</para></entry>
</row>
+ <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>
@@ -1436,6 +1457,23 @@ SELECT NOT(ROW(table.*) IS NOT NULL) FROM TABLE; -- detect at least one null in
</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
index 2bc31eabf2..a46ebffae4 100644
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -2786,6 +2786,93 @@ 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;
+
+ /*
+ * Handle NaN and Inf cases explicitly. This simplifies the overflow
+ * checks on platforms that do not set errno.
+ */
+ if (isnan(arg1))
+ result = arg1;
+ else if (isinf(arg1))
+ {
+ /* Per POSIX, an input of -Inf causes a domain error */
+ if (arg1 < 0)
+ {
+ float_overflow_error();
+ result = get_float8_nan(); /* keep compiler quiet */
+ }
+ else
+ result = arg1;
+ }
+ else
+ {
+ /*
+ * Note: the POSIX/C99 gamma function is called "tgamma", not "gamma".
+ *
+ * On some platforms, tgamma() will not set errno but just return Inf,
+ * NaN, or zero to report overflow/underflow; therefore, test those
+ * cases explicitly (note that, like the exponential function, the
+ * gamma function has no zeros).
+ */
+ errno = 0;
+ result = tgamma(arg1);
+
+ if (errno != 0 || isinf(result) || isnan(result))
+ {
+ if (result != 0.0)
+ float_overflow_error();
+ else
+ float_underflow_error();
+ }
+ else if (result == 0.0)
+ float_underflow_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.
+ *
+ * On some platforms, lgamma() will not set errno but just return infinity
+ * to report overflow, but it should never underflow.
+ */
+ if (errno == ERANGE || (isinf(result) && !isinf(arg1)))
+ float_overflow_error();
+
+ PG_RETURN_FLOAT8(result);
+}
+
+
/*
* =========================
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
index cd9422d0ba..2149fa4281 100644
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -3542,6 +3542,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
index 9ef9793fe9..10a5a6e1b6 100644
--- a/src/test/regress/expected/float8.out
+++ b/src/test/regress/expected/float8.out
@@ -829,6 +829,63 @@ FROM (VALUES (float8 '-infinity'),
NaN | NaN | NaN
(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),
+ (float8 'infinity'), (float8 'nan')) 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
+ Infinity | Infinity | Infinity
+ NaN | NaN | NaN
+(8 rows)
+
+-- test overflow/underflow handling
+SELECT gamma(float8 '-infinity');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '-infinity');
+ lgamma
+----------
+ Infinity
+(1 row)
+
+SELECT gamma(float8 '-1000.5');
+ERROR: value out of range: underflow
+SELECT lgamma(float8 '-1000.5');
+ lgamma
+------------------
+ -5914.4377011169
+(1 row)
+
+SELECT gamma(float8 '-1');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '-1');
+ERROR: value out of range: overflow
+SELECT gamma(float8 '0');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '0');
+ERROR: value out of range: overflow
+SELECT gamma(float8 '1000');
+ERROR: value out of range: overflow
+SELECT lgamma(float8 '1000');
+ lgamma
+-----------------
+ 5905.2204232092
+(1 row)
+
+SELECT lgamma(float8 '1e308');
+ERROR: value out of range: overflow
RESET extra_float_digits;
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
diff --git a/src/test/regress/sql/float8.sql b/src/test/regress/sql/float8.sql
index 81a35e0bf1..db8d5724c2 100644
--- a/src/test/regress/sql/float8.sql
+++ b/src/test/regress/sql/float8.sql
@@ -245,6 +245,29 @@ 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),
+ (float8 'infinity'), (float8 'nan')) AS t(x);
+-- test overflow/underflow handling
+SELECT gamma(float8 '-infinity');
+SELECT lgamma(float8 '-infinity');
+SELECT gamma(float8 '-1000.5');
+SELECT lgamma(float8 '-1000.5');
+SELECT gamma(float8 '-1');
+SELECT lgamma(float8 '-1');
+SELECT gamma(float8 '0');
+SELECT lgamma(float8 '0');
+SELECT gamma(float8 '1000');
+SELECT lgamma(float8 '1000');
+SELECT lgamma(float8 '1e308');
+RESET extra_float_digits;
+
-- test for over- and underflow
INSERT INTO FLOAT8_TBL(f1) VALUES ('10e400');
--
2.43.0
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.a.rasheed@gmail.com>
wrote:
On Thu, 14 Nov 2024 at 22:35, Dean Rasheed <dean.a.rasheed@gmail.com>
wrote:On Thu, 14 Nov 2024 at 16:28, Dmitry Koval <d.koval@postgrespro.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
Attachments:
On Tue, 25 Mar 2025 at 05:01, Alexandra Wang
<alexandra.wang.oss@gmail.com> wrote:
I reviewed the patch and it looks good to me. I’ve run some tests, and
everything works as expected.
Thanks for the careful review and testing.
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)
Thanks for that. It's very useful to see how tgamma() and lgamma()
behave on those different platforms.
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.
Hmm, I don't think "division by zero" would be the right error to use,
if we did that. That's probably exposing some internal implementation
detail, but I think it would look odd to users. Perhaps "input is out
of range" would work, but I don't think it's really necessary to
distinguish between these different types of errors.
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.
OK, thanks. I've updated that comment and committed it.
Regards,
Dean
Em qua., 26 de mar. de 2025 às 06:50, Dean Rasheed <dean.a.rasheed@gmail.com>
escreveu:
OK, thanks. I've updated that comment and committed it.
func.sgml doesn't mention lgamma function, so users will think just gamma
exists, is that correct ?
regards
Marcos