v5-0001-Add-support-for-gamma-and-lgamma-functions.patch

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