Re: Greatest Common Divisor

diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
index 57a1539506..e2b7d6d240 100644
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -870,6 +870,32 @@
        <entry><literal>-43</literal></entry>
       </row>
 
+      <row>
+       <entry>
+        <indexterm>
+         <primary>gcd</primary>
+        </indexterm>
+        <literal><function>gcd(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</function></literal>
+       </entry>
+       <entry>(same as argument types)</entry>
+       <entry>greatest common divisor</entry>
+       <entry><literal>gcd(1071, 462)</literal></entry>
+       <entry><literal>21</literal></entry>
+      </row>
+
+      <row>
+       <entry>
+        <indexterm>
+         <primary>lcm</primary>
+        </indexterm>
+        <literal><function>lcm(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</function></literal>
+       </entry>
+       <entry>(same as argument types)</entry>
+       <entry>least common multiple</entry>
+       <entry><literal>lcm(1071, 462)</literal></entry>
+       <entry><literal>23562</literal></entry>
+      </row>
+
       <row>
        <entry>
         <indexterm>
diff --git a/src/backend/utils/adt/int.c b/src/backend/utils/adt/int.c
index 583ce71e66..5244b891dc 100644
--- a/src/backend/utils/adt/int.c
+++ b/src/backend/utils/adt/int.c
@@ -1196,6 +1196,226 @@ int2abs(PG_FUNCTION_ARGS)
 	PG_RETURN_INT16(result);
 }
 
+/*
+ * Greatest Common Divisor
+ *
+ * Special cases:
+ *   - gcd(x, 0) = gcd(0, x) = abs(x)
+ *   		because 0 is divisible by anything
+ *   - gcd(0, 0) = 0
+ *   		complies with the previous definition and is a common convention
+ *
+ * The following cases involving INT_MIN have two possible results.  They could
+ * return INT_MIN because INT_MIN is a valid divisor of INT_MIN, or they could
+ * throw an exception because the result is negative.
+ * The consensus is to throw an exception.
+ *
+ *   - gcd(INT_MIN, 0)
+ *   - gcd(INT_MIN, INT_MIN)
+ *
+ * Any other value with INT_MIN will be a positive value representable within
+ * the data type.
+ */
+static int32
+int4gcd_internal(int32 arg1, int32 arg2)
+{
+	int32	swap;
+	int32	a1, a2;
+
+	/*
+	 * Put the greater absolute value in arg1.
+	 *
+	 * This would happen automatically in the loop below, but avoids an
+	 * expensive modulo simulation on some architectures.
+	 *
+	 * We do this in negative space in order to handle INT_MIN.
+	 */
+	a1 = (arg1 < 0) ? arg1 : -arg1;
+	a2 = (arg2 < 0) ? arg2 : -arg2;
+	if (a1 > a2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	/* Special care needs to be taken with INT_MIN.  See comments above. */
+	if (arg1 == PG_INT32_MIN)
+	{
+		if (arg2 == 0 || arg2 == PG_INT32_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("integer out of range")));
+		/*
+		 * Making arg2 positive avoids the INT_MIN % -1 issue described in
+		 * int4mod.
+		 */
+		arg2 = -arg2;
+	}
+
+	/* Find GCD using the basic Euclidean algorithm */
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	/*
+	 * Make sure the result is positive. (We know we don't have INT_MIN
+	 * anymore).
+	 */
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	return arg1;
+}
+
+static int16
+int2gcd_internal(int16 arg1, int16 arg2)
+{
+	/* See int4gcd_internal for commented version. */
+	int16	swap;
+	int16	a1, a2;
+
+	a1 = (arg1 < 0) ? arg1 : -arg1;
+	a2 = (arg2 < 0) ? arg2 : -arg2;
+	if (a1 > a2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	if (arg1 == PG_INT16_MIN)
+	{
+		if (arg2 == 0 || arg2 == PG_INT16_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("smallint out of range")));
+		arg2 = -arg2;
+	}
+
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	return arg1;
+}
+
+Datum
+int4gcd(PG_FUNCTION_ARGS)
+{
+	int32	arg1 = PG_GETARG_INT32(0);
+	int32	arg2 = PG_GETARG_INT32(1);
+	int32	result;
+
+	result = int4gcd_internal(arg1, arg2);
+
+	PG_RETURN_INT32(result);
+}
+
+Datum
+int2gcd(PG_FUNCTION_ARGS)
+{
+	int16	arg1 = PG_GETARG_INT16(0);
+	int16	arg2 = PG_GETARG_INT16(1);
+	int16	result;
+
+	result = int2gcd_internal(arg1, arg2);
+
+	PG_RETURN_INT16(result);
+}
+
+/*
+ * Least Common Multiple
+ */
+
+Datum
+int4lcm(PG_FUNCTION_ARGS)
+{
+	int32	arg1 = PG_GETARG_INT32(0);
+	int32	arg2 = PG_GETARG_INT32(1);
+	int32	gcd;
+	int32	result;
+
+	/* lcm(n, 0) = lcm(0, n) = 0 */
+	if (arg1 == 0 || arg2 == 0)
+		PG_RETURN_INT32(0);
+
+	/* lcm(n, n) = abs(n) */
+	if (arg1 == arg2)
+	{
+		if (arg1 == PG_INT32_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("integer out of range")));
+
+		PG_RETURN_INT32((arg1 < 0) ? -arg1 : arg1);
+	}
+
+	/* lcm(m, n) = abs(m / gcd(m, n) * n) */
+	gcd = int4gcd_internal(arg1, arg2);
+	arg1 = arg1 / gcd;
+
+	if (unlikely(pg_mul_s32_overflow(arg1, arg2, &result)))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("integer out of range")));
+
+	/* If the result is INT_MIN, it cannot be represented. */
+	if (result == PG_INT32_MIN)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("integer out of range")));
+
+	PG_RETURN_INT32((result < 0) ? -result : result);
+}
+
+Datum
+int2lcm(PG_FUNCTION_ARGS)
+{
+	/* See int4lcm for commented version. */
+	int16	arg1 = PG_GETARG_INT16(0);
+	int16	arg2 = PG_GETARG_INT16(1);
+	int16	gcd;
+	int16	result;
+
+	if (arg1 == 0 || arg2 == 0)
+		PG_RETURN_INT16(0);
+
+	if (arg1 == arg2)
+	{
+		if (arg1 == PG_INT16_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("smallint out of range")));
+
+		PG_RETURN_INT16((arg1 < 0) ? -arg1 : arg1);
+	}
+
+	gcd = int2gcd_internal(arg1, arg2);
+	arg1 = arg1 / gcd;
+
+	if (unlikely(pg_mul_s16_overflow(arg1, arg2, &result)))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("smallint out of range")));
+
+	if (result == PG_INT16_MIN)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("smallint out of range")));
+
+	PG_RETURN_INT16((result < 0) ? -result : result);
+}
+
 Datum
 int2larger(PG_FUNCTION_ARGS)
 {
diff --git a/src/backend/utils/adt/int8.c b/src/backend/utils/adt/int8.c
index fcdf77331e..96c65688a5 100644
--- a/src/backend/utils/adt/int8.c
+++ b/src/backend/utils/adt/int8.c
@@ -667,6 +667,104 @@ int8mod(PG_FUNCTION_ARGS)
 	PG_RETURN_INT64(arg1 % arg2);
 }
 
+/*
+ * Greatest Common Divisor
+ *
+ * See int4gcd_internal for commented version.
+ */
+
+static int64
+int8gcd_internal(int64 arg1, int64 arg2)
+{
+	int64	swap;
+	int64	a1, a2;
+
+	a1 = (arg1 < 0) ? arg1 : -arg1;
+	a2 = (arg2 < 0) ? arg2 : -arg2;
+	if (a1 > a2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	if (arg1 == PG_INT64_MIN)
+	{
+		if (arg2 == 0 || arg2 == PG_INT64_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("bigint out of range")));
+		arg2 = -arg2;
+	}
+
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	return arg1;
+}
+
+Datum
+int8gcd(PG_FUNCTION_ARGS)
+{
+	int64	arg1 = PG_GETARG_INT64(0);
+	int64	arg2 = PG_GETARG_INT64(1);
+	int64	result;
+
+	result = int8gcd_internal(arg1, arg2);
+
+	PG_RETURN_INT64(result);
+}
+
+/*
+ * Least Common Multiple
+ *
+ * See int4lcm for commented version.
+ */
+
+Datum
+int8lcm(PG_FUNCTION_ARGS)
+{
+	int64	arg1 = PG_GETARG_INT64(0);
+	int64	arg2 = PG_GETARG_INT64(1);
+	int64	gcd;
+	int64	result;
+
+	if (arg1 == 0 || arg2 == 0)
+		PG_RETURN_INT64(0);
+
+	if (arg1 == arg2)
+	{
+		if (arg1 == PG_INT64_MIN)
+			ereport(ERROR,
+					(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+					 errmsg("bigint out of range")));
+
+		PG_RETURN_INT64((arg1 < 0) ? -arg1 : arg1);
+	}
+
+	gcd = int8gcd_internal(arg1, arg2);
+	arg1 = arg1 / gcd;
+
+	if (unlikely(pg_mul_s64_overflow(arg1, arg2, &result)))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("bigint out of range")));
+
+	if (result == PG_INT64_MIN)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("bigint out of range")));
+
+	PG_RETURN_INT64((result < 0) ? -result : result);
+}
+
 
 Datum
 int8inc(PG_FUNCTION_ARGS)
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 14054272c8..d286004db0 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -220,6 +220,8 @@ struct NumericData
 	 | ((n)->choice.n_short.n_header & NUMERIC_SHORT_WEIGHT_MASK)) \
 	: ((n)->choice.n_long.n_weight))
 
+#define NUMERIC_IS_INTEGRAL(n)	(NUMERIC_NDIGITS(n) <= (NUMERIC_WEIGHT(n) + 1))
+
 /* ----------
  * NumericVar is the format we use for arithmetic.  The digit-array part
  * is the same as the NumericData storage format, but the header is more
@@ -2691,6 +2693,141 @@ numeric_div_trunc(PG_FUNCTION_ARGS)
 }
 
 
+/*
+ * Greatest Common Divisor
+ *
+ * We define here that gcd(n, 0) = n and therefore gdc(0, 0) = 0.
+ * See the comments on int[24]gcd for more details.
+ */
+static void
+gcd_var(NumericVar arg1, NumericVar arg2, NumericVar *result)
+{
+	NumericVar	swap;
+	int			cmp;
+
+	/*
+	 * Unlike the integer types, there are no negative numerics that cannot be
+	 * represented positively, so just abs them from the beginning.
+	 */
+	arg1.sign = NUMERIC_POS;
+	arg2.sign = NUMERIC_POS;
+
+	/* gcd(a, a) = a */
+	cmp = cmp_var(&arg1, &arg2);
+	if (cmp == 0)
+	{
+		set_var_from_var(&arg1, result);
+		return;
+	}
+
+	init_var(&swap);
+
+	/* Save ourselves a call to mod_var if arg1 < arg2 */
+	if (cmp == -1)
+	{
+		set_var_from_var(&arg1, &swap);
+		set_var_from_var(&arg2, &arg1);
+		set_var_from_var(&swap, &arg2);
+	}
+
+	/* Use basic Euclidean algorithm for GCD */
+	while (arg2.ndigits != 0)
+	{
+		/* this loop can take awhile, so allow it to be interrupted */
+		CHECK_FOR_INTERRUPTS();
+
+		mod_var(&arg1, &arg2, &swap);
+		set_var_from_var(&arg2, &arg1);
+		set_var_from_var(&swap, &arg2);
+	}
+
+	free_var(&swap);
+
+	set_var_from_var(&arg1, result);
+}
+
+Datum
+numeric_gcd(PG_FUNCTION_ARGS)
+{
+	Numeric		num1 = PG_GETARG_NUMERIC(0);
+	Numeric		num2 = PG_GETARG_NUMERIC(1);
+	Numeric		res;
+	NumericVar	arg1;
+	NumericVar	arg2;
+	NumericVar	result;
+
+	if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
+		PG_RETURN_NUMERIC(make_result(&const_nan));
+
+	init_var_from_num(num1, &arg1);
+	init_var_from_num(num2, &arg2);
+
+	init_var(&result);
+
+	gcd_var(arg1, arg2, &result);
+	res = make_result_opt_error(&result, NULL);
+
+	free_var(&result);
+
+	PG_RETURN_NUMERIC(res);
+}
+
+/*
+ * Least Common Multiple
+ */
+
+Datum
+numeric_lcm(PG_FUNCTION_ARGS)
+{
+	Numeric		num1 = PG_GETARG_NUMERIC(0);
+	Numeric		num2 = PG_GETARG_NUMERIC(1);
+	Numeric		res;
+	NumericVar	arg1;
+	NumericVar	arg2;
+	NumericVar	result;
+
+	/* If we get a NaN, just return NaN */
+	if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
+		PG_RETURN_NUMERIC(make_result(&const_nan));
+
+	/* Values must be integral */
+	if (!NUMERIC_IS_INTEGRAL(num1) || !NUMERIC_IS_INTEGRAL(num2))
+		ereport(ERROR,
+				(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+				 errmsg("parameters to lcm must be integral")));
+
+	/* If either argument is 0, the result is 0 */
+	if (NUMERIC_NDIGITS(num1) == 0 || NUMERIC_NDIGITS(num2) == 0)
+		PG_RETURN_NUMERIC(make_result(&const_zero));
+
+	init_var_from_num(num1, &arg1);
+	init_var_from_num(num2, &arg2);
+
+	/*
+	 * Make arguments positive.  We can do this because all negative numbers
+	 * are representable in the positive space.
+	 */
+	arg1.sign = NUMERIC_POS;
+	arg2.sign = NUMERIC_POS;
+
+	/* If both arguments are the same, just return it */
+	if (cmp_var(&arg1, &arg2) == 0)
+		PG_RETURN_NUMERIC(num1);
+
+	init_var(&result);
+
+	/* result = arg1 / gcd(arg1, arg2) * arg2 */
+	gcd_var(arg1, arg2, &result);
+	div_var(&arg1, &result, &result, 0, false);
+	mul_var(&arg2, &result, &result, 0);
+
+	res = make_result_opt_error(&result, NULL);
+
+	free_var(&result);
+
+	PG_RETURN_NUMERIC(res);
+}
+
 /*
  * numeric_mod() -
  *
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
index 0b6045acb1..334fe30d44 100644
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -10729,4 +10729,32 @@
   proname => 'pg_partition_root', prorettype => 'regclass',
   proargtypes => 'regclass', prosrc => 'pg_partition_root' },
 
+# greatest common divisor
+{ oid => '8463', descr => 'greatest common divisor',
+  proname => 'gcd', prorettype => 'int2', proargtypes => 'int2 int2',
+  prosrc => 'int2gcd' },
+{ oid => '8464', descr => 'greatest common divisor',
+  proname => 'gcd', prorettype => 'int4', proargtypes => 'int4 int4',
+  prosrc => 'int4gcd' },
+{ oid => '8465', descr => 'greatest common divisor',
+  proname => 'gcd', prorettype => 'int8', proargtypes => 'int8 int8',
+  prosrc => 'int8gcd' },
+{ oid => '8466', descr => 'greatest common divisor',
+  proname => 'gcd', prorettype => 'numeric', proargtypes => 'numeric numeric',
+  prosrc => 'numeric_gcd' },
+
+# least common multiple
+{ oid => '8467', descr => 'least common multiple',
+  proname => 'lcm', prorettype => 'int2', proargtypes => 'int2 int2',
+  prosrc => 'int2lcm' },
+{ oid => '8468', descr => 'least common multiple',
+  proname => 'lcm', prorettype => 'int4', proargtypes => 'int4 int4',
+  prosrc => 'int4lcm' },
+{ oid => '8469', descr => 'least common multiple',
+  proname => 'lcm', prorettype => 'int8', proargtypes => 'int8 int8',
+  prosrc => 'int8lcm' },
+{ oid => '8470', descr => 'least common multiple',
+  proname => 'lcm', prorettype => 'numeric', proargtypes => 'numeric numeric',
+  prosrc => 'numeric_lcm' },
+
 ]
diff --git a/src/test/regress/expected/int2.out b/src/test/regress/expected/int2.out
index 8c255b9e4d..a3dc104565 100644
--- a/src/test/regress/expected/int2.out
+++ b/src/test/regress/expected/int2.out
@@ -306,3 +306,73 @@ FROM (VALUES (-2.5::numeric),
   2.5 |          3
 (7 rows)
 
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (24948::int2, 4914::int2)) AS v(a, b);
+ gcd | gcd | gcd | gcd | gcd 
+-----+-----+-----+-----+-----
+ 378 | 378 | 378 | 378 | 378
+(1 row)
+
+SELECT gcd((-32768)::int2, 16384::int2);
+  gcd  
+-------
+ 16384
+(1 row)
+
+SELECT gcd((-32768)::int2, (-1)::int2);
+ gcd 
+-----
+   1
+(1 row)
+
+SELECT gcd((-32768)::int2, 0::int2); -- fail
+ERROR:  smallint out of range
+SELECT gcd((-32768)::int2, (-32768)::int2); -- fail
+ERROR:  smallint out of range
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int2, 462::int2)) AS v(a, b);
+ lcm  | lcm  | lcm  | lcm  | lcm  
+------+------+------+------+------
+ 2310 | 2310 | 2310 | 2310 | 2310
+(1 row)
+
+SELECT lcm(42::int2, 42::int2);
+ lcm 
+-----
+  42
+(1 row)
+
+SELECT lcm(42::int2, 0::int2);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int2, 42::int2);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int2, 0::int2);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm((-32768)::int2, 0::int2);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int2, (-32768)::int2);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(32767::int2, 32766::int2); -- fail
+ERROR:  smallint out of range
diff --git a/src/test/regress/expected/int4.out b/src/test/regress/expected/int4.out
index bda7a8daef..83f3ed2ab3 100644
--- a/src/test/regress/expected/int4.out
+++ b/src/test/regress/expected/int4.out
@@ -403,3 +403,73 @@ FROM (VALUES (-2.5::numeric),
   2.5 |          3
 (7 rows)
 
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (61866666::int4, 6410818::int4)) AS v(a, b);
+ gcd  | gcd  | gcd  | gcd  | gcd  
+------+------+------+------+------
+ 1466 | 1466 | 1466 | 1466 | 1466
+(1 row)
+
+SELECT gcd((-2147483648)::int4, 1073741824::int4);
+    gcd     
+------------
+ 1073741824
+(1 row)
+
+SELECT gcd((-2147483648)::int4, (-1)::int4);
+ gcd 
+-----
+   1
+(1 row)
+
+SELECT gcd((-2147483648)::int4, 0::int4); -- fail
+ERROR:  integer out of range
+SELECT gcd((-2147483648)::int4, (-2147483648)::int4); -- fail
+ERROR:  integer out of range
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int4, 462::int4)) AS v(a, b);
+ lcm  | lcm  | lcm  | lcm  | lcm  
+------+------+------+------+------
+ 2310 | 2310 | 2310 | 2310 | 2310
+(1 row)
+
+SELECT lcm(42::int4, 42::int4);
+ lcm 
+-----
+  42
+(1 row)
+
+SELECT lcm(42::int4, 0::int4);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int4, 42::int4);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int4, 0::int4);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm((-2147483648)::int4, 0::int4);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int4, (-2147483648)::int4);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(2147483647::int4, 2147483646::int4); -- fail
+ERROR:  integer out of range
diff --git a/src/test/regress/expected/int8.out b/src/test/regress/expected/int8.out
index 8447a28c3d..fe6680ef2c 100644
--- a/src/test/regress/expected/int8.out
+++ b/src/test/regress/expected/int8.out
@@ -886,3 +886,73 @@ FROM (VALUES (-2.5::numeric),
   2.5 |          3
 (7 rows)
 
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (288484263558::int8, 29893644334::int8)) AS v(a, b);
+   gcd   |   gcd   |   gcd   |   gcd   |   gcd   
+---------+---------+---------+---------+---------
+ 6835958 | 6835958 | 6835958 | 6835958 | 6835958
+(1 row)
+
+SELECT gcd((-9223372036854775808)::int8, 4611686018427387904::int8);
+         gcd         
+---------------------
+ 4611686018427387904
+(1 row)
+
+SELECT gcd((-9223372036854775808)::int8, (-1)::int8);
+ gcd 
+-----
+   1
+(1 row)
+
+SELECT gcd((-9223372036854775808)::int8, 0::int8); -- fail
+ERROR:  bigint out of range
+SELECT gcd((-9223372036854775808)::int8, (-9223372036854775808)::int8); -- fail
+ERROR:  bigint out of range
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int8, 462::int8)) AS v(a, b);
+ lcm  | lcm  | lcm  | lcm  | lcm  
+------+------+------+------+------
+ 2310 | 2310 | 2310 | 2310 | 2310
+(1 row)
+
+SELECT lcm(42::int8, 42::int8);
+ lcm 
+-----
+  42
+(1 row)
+
+SELECT lcm(42::int8, 0::int8);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int8, 42::int8);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int8, 0::int8);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm((-9223372036854775808)::int8, 0::int8);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::int8, (-9223372036854775808)::int8);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(9223372036854775807::int8, 9223372036854775806::int8); -- fail
+ERROR:  bigint out of range
diff --git a/src/test/regress/expected/numeric.out b/src/test/regress/expected/numeric.out
index 1cb3c3bfab..84330d68a4 100644
--- a/src/test/regress/expected/numeric.out
+++ b/src/test/regress/expected/numeric.out
@@ -2094,3 +2094,71 @@ SELECT SUM((-9999)::numeric) FROM generate_series(1, 100000);
  -999900000
 (1 row)
 
+--
+-- Tests for GCD()
+--
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (330::numeric, 462::numeric)) AS v(a, b);
+ gcd | gcd | gcd | gcd | gcd 
+-----+-----+-----+-----+-----
+  66 |  66 |  66 |  66 |  66
+(1 row)
+
+SELECT gcd(0::numeric, 0::numeric);
+ gcd 
+-----
+   0
+(1 row)
+
+SELECT gcd(330.3::numeric, 462::numeric);
+ gcd 
+-----
+ 0.3
+(1 row)
+
+SELECT gcd(330::numeric, 462.5::numeric);
+ gcd 
+-----
+ 2.5
+(1 row)
+
+--
+-- Tests for LCM()
+--
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::numeric, 462::numeric)) AS v(a, b);
+ lcm  | lcm  | lcm  | lcm  | lcm  
+------+------+------+------+------
+ 2310 | 2310 | 2310 | 2310 | 2310
+(1 row)
+
+SELECT lcm(42::numeric, 42::numeric);
+ lcm 
+-----
+  42
+(1 row)
+
+SELECT lcm(42::numeric, 0::numeric);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::numeric, 42::numeric);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(0::numeric, 0::numeric);
+ lcm 
+-----
+   0
+(1 row)
+
+SELECT lcm(330.3::numeric, 462::numeric); -- fails
+ERROR:  parameters to lcm must be integral
+SELECT lcm(330::numeric, 462.5::numeric); -- fails
+ERROR:  parameters to lcm must be integral
+SELECT lcm(9999 * (10::numeric)^131068 + (10::numeric^131068 - 1), 2); -- fails
+ERROR:  value overflows numeric format
diff --git a/src/test/regress/sql/int2.sql b/src/test/regress/sql/int2.sql
index 7dbafb6dac..18b068e23f 100644
--- a/src/test/regress/sql/int2.sql
+++ b/src/test/regress/sql/int2.sql
@@ -112,3 +112,22 @@ FROM (VALUES (-2.5::numeric),
              (0.5::numeric),
              (1.5::numeric),
              (2.5::numeric)) t(x);
+
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (24948::int2, 4914::int2)) AS v(a, b);
+SELECT gcd((-32768)::int2, 16384::int2);
+SELECT gcd((-32768)::int2, (-1)::int2);
+SELECT gcd((-32768)::int2, 0::int2); -- fail
+SELECT gcd((-32768)::int2, (-32768)::int2); -- fail
+
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int2, 462::int2)) AS v(a, b);
+SELECT lcm(42::int2, 42::int2);
+SELECT lcm(42::int2, 0::int2);
+SELECT lcm(0::int2, 42::int2);
+SELECT lcm(0::int2, 0::int2);
+SELECT lcm((-32768)::int2, 0::int2);
+SELECT lcm(0::int2, (-32768)::int2);
+SELECT lcm(32767::int2, 32766::int2); -- fail
diff --git a/src/test/regress/sql/int4.sql b/src/test/regress/sql/int4.sql
index f014cb2d32..0a0e2796b2 100644
--- a/src/test/regress/sql/int4.sql
+++ b/src/test/regress/sql/int4.sql
@@ -155,3 +155,22 @@ FROM (VALUES (-2.5::numeric),
              (0.5::numeric),
              (1.5::numeric),
              (2.5::numeric)) t(x);
+
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (61866666::int4, 6410818::int4)) AS v(a, b);
+SELECT gcd((-2147483648)::int4, 1073741824::int4);
+SELECT gcd((-2147483648)::int4, (-1)::int4);
+SELECT gcd((-2147483648)::int4, 0::int4); -- fail
+SELECT gcd((-2147483648)::int4, (-2147483648)::int4); -- fail
+
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int4, 462::int4)) AS v(a, b);
+SELECT lcm(42::int4, 42::int4);
+SELECT lcm(42::int4, 0::int4);
+SELECT lcm(0::int4, 42::int4);
+SELECT lcm(0::int4, 0::int4);
+SELECT lcm((-2147483648)::int4, 0::int4);
+SELECT lcm(0::int4, (-2147483648)::int4);
+SELECT lcm(2147483647::int4, 2147483646::int4); -- fail
diff --git a/src/test/regress/sql/int8.sql b/src/test/regress/sql/int8.sql
index e890452236..def0c17f4f 100644
--- a/src/test/regress/sql/int8.sql
+++ b/src/test/regress/sql/int8.sql
@@ -225,3 +225,22 @@ FROM (VALUES (-2.5::numeric),
              (0.5::numeric),
              (1.5::numeric),
              (2.5::numeric)) t(x);
+
+-- gcd
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (288484263558::int8, 29893644334::int8)) AS v(a, b);
+SELECT gcd((-9223372036854775808)::int8, 4611686018427387904::int8);
+SELECT gcd((-9223372036854775808)::int8, (-1)::int8);
+SELECT gcd((-9223372036854775808)::int8, 0::int8); -- fail
+SELECT gcd((-9223372036854775808)::int8, (-9223372036854775808)::int8); -- fail
+
+-- lcm
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::int8, 462::int8)) AS v(a, b);
+SELECT lcm(42::int8, 42::int8);
+SELECT lcm(42::int8, 0::int8);
+SELECT lcm(0::int8, 42::int8);
+SELECT lcm(0::int8, 0::int8);
+SELECT lcm((-9223372036854775808)::int8, 0::int8);
+SELECT lcm(0::int8, (-9223372036854775808)::int8);
+SELECT lcm(9223372036854775807::int8, 9223372036854775806::int8); -- fail
diff --git a/src/test/regress/sql/numeric.sql b/src/test/regress/sql/numeric.sql
index a939412359..73e35d0173 100644
--- a/src/test/regress/sql/numeric.sql
+++ b/src/test/regress/sql/numeric.sql
@@ -1043,3 +1043,25 @@ select scale(-13.000000000000000);
 -- cases that need carry propagation
 SELECT SUM(9999::numeric) FROM generate_series(1, 100000);
 SELECT SUM((-9999)::numeric) FROM generate_series(1, 100000);
+
+--
+-- Tests for GCD()
+--
+SELECT gcd(a, b), gcd(a, -b), gcd(-a, b), gcd(-a, -b), gcd(b, a)
+FROM (VALUES (330::numeric, 462::numeric)) AS v(a, b);
+SELECT gcd(0::numeric, 0::numeric);
+SELECT gcd(330.3::numeric, 462::numeric);
+SELECT gcd(330::numeric, 462.5::numeric);
+
+--
+-- Tests for LCM()
+--
+SELECT lcm(a, b), lcm(a, -b), lcm(-a, b), lcm(-a, -b), lcm(b, a)
+FROM (VALUES (330::numeric, 462::numeric)) AS v(a, b);
+SELECT lcm(42::numeric, 42::numeric);
+SELECT lcm(42::numeric, 0::numeric);
+SELECT lcm(0::numeric, 42::numeric);
+SELECT lcm(0::numeric, 0::numeric);
+SELECT lcm(330.3::numeric, 462::numeric); -- fails
+SELECT lcm(330::numeric, 462.5::numeric); -- fails
+SELECT lcm(9999 * (10::numeric)^131068 + (10::numeric^131068 - 1), 2); -- fails