Add functions gcd() and lcm() for integer and numeric types.

These compute the greatest common divisor and least common multiple of
a pair of numbers using the Euclidean algorithm.

Vik Fearing, reviewed by Fabien Coelho.

Discussion: https://postgr.es/m/adbd3e0b-e3f1-5bbc-21db-03caf1cef0f7@2ndquadrant.com

src/backend/utils/adt/int8.c

@@ -667,6 +667,132 @@ int8mod(PG_FUNCTION_ARGS)
 	PG_RETURN_INT64(arg1 % arg2);
 }
 
+/*
+ * Greatest Common Divisor
+ *
+ * Returns the largest positive integer that exactly divides both inputs.
+ * 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
+ *
+ * Special care must be taken if either input is INT64_MIN ---
+ * gcd(0, INT64_MIN), gcd(INT64_MIN, 0) and gcd(INT64_MIN, INT64_MIN) are
+ * all equal to abs(INT64_MIN), which cannot be represented as a 64-bit signed
+ * integer.
+ */
+static int64
+int8gcd_internal(int64 arg1, int64 arg2)
+{
+	int64	swap;
+	int64	a1, a2;
+
+	/*
+	 * Put the greater absolute value in arg1.
+	 *
+	 * This would happen automatically in the loop below, but avoids an
+	 * expensive modulo operation, and simplifies the special-case handling
+	 * for INT64_MIN below.
+	 *
+	 * We do this in negative space in order to handle INT64_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 INT64_MIN.  See comments above. */
+	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")));
+
+		/*
+		 * Some machines throw a floating-point exception for INT64_MIN % -1,
+		 * which is a bit silly since the correct answer is perfectly
+		 * well-defined, namely zero.  Guard against this and just return the
+		 * result, gcd(INT64_MIN, -1) = 1.
+		 */
+		if (arg2 == -1)
+			return 1;
+	}
+
+	/* Use the Euclidean algorithm to find the GCD */
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	/*
+	 * Make sure the result is positive. (We know we don't have INT64_MIN
+	 * anymore).
+	 */
+	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
+ */
+Datum
+int8lcm(PG_FUNCTION_ARGS)
+{
+	int64	arg1 = PG_GETARG_INT64(0);
+	int64	arg2 = PG_GETARG_INT64(1);
+	int64	gcd;
+	int64	result;
+
+	/*
+	 * Handle lcm(x, 0) = lcm(0, x) = 0 as a special case.  This prevents a
+	 * division-by-zero error below when x is zero, and an overflow error from
+	 * the GCD computation when x = INT64_MIN.
+	 */
+	if (arg1 == 0 || arg2 == 0)
+		PG_RETURN_INT64(0);
+
+	/* lcm(x, y) = abs(x / gcd(x, y) * y) */
+	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 the result is INT64_MIN, it cannot be represented. */
+	if (unlikely(result == PG_INT64_MIN))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("bigint out of range")));
+
+	if (result < 0)
+		result = -result;
+
+	PG_RETURN_INT64(result);
+}
 
 Datum
 int8inc(PG_FUNCTION_ARGS)