Skip to content

Commit 13661dd

Browse files
deanrasheeddeanrasheed
committed
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
1 parent 530609a commit 13661dd

File tree

12 files changed

+689
-1
lines changed

12 files changed

+689
-1
lines changed

‎doc/src/sgml/func.sgml

+34
Original file line numberDiff line numberDiff line change
@@ -870,6 +870,40 @@
870870
<entry><literal>-43</literal></entry>
871871
</row>
872872

873+
<row>
874+
<entry>
875+
<indexterm>
876+
<primary>gcd</primary>
877+
</indexterm>
878+
<literal><function>gcd(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</function></literal>
879+
</entry>
880+
<entry>(same as argument types)</entry>
881+
<entry>
882+
greatest common divisor (the largest positive number that divides both
883+
inputs with no remainder); returns <literal>0</literal> if both inputs
884+
are zero
885+
</entry>
886+
<entry><literal>gcd(1071, 462)</literal></entry>
887+
<entry><literal>21</literal></entry>
888+
</row>
889+
890+
<row>
891+
<entry>
892+
<indexterm>
893+
<primary>lcm</primary>
894+
</indexterm>
895+
<literal><function>lcm(<replaceable>a</replaceable>, <replaceable>b</replaceable>)</function></literal>
896+
</entry>
897+
<entry>(same as argument types)</entry>
898+
<entry>
899+
least common multiple (the smallest strictly positive number that is
900+
an integral multiple of both inputs); returns <literal>0</literal> if
901+
either input is zero
902+
</entry>
903+
<entry><literal>lcm(1071, 462)</literal></entry>
904+
<entry><literal>23562</literal></entry>
905+
</row>
906+
873907
<row>
874908
<entry>
875909
<indexterm>

‎src/backend/utils/adt/int.c

+126
Original file line numberDiff line numberDiff line change
@@ -1196,6 +1196,132 @@ int2abs(PG_FUNCTION_ARGS)
11961196
PG_RETURN_INT16(result);
11971197
}
11981198

1199+
/*
1200+
* Greatest Common Divisor
1201+
*
1202+
* Returns the largest positive integer that exactly divides both inputs.
1203+
* Special cases:
1204+
* - gcd(x, 0) = gcd(0, x) = abs(x)
1205+
* because 0 is divisible by anything
1206+
* - gcd(0, 0) = 0
1207+
* complies with the previous definition and is a common convention
1208+
*
1209+
* Special care must be taken if either input is INT_MIN --- gcd(0, INT_MIN),
1210+
* gcd(INT_MIN, 0) and gcd(INT_MIN, INT_MIN) are all equal to abs(INT_MIN),
1211+
* which cannot be represented as a 32-bit signed integer.
1212+
*/
1213+
static int32
1214+
int4gcd_internal(int32 arg1, int32 arg2)
1215+
{
1216+
int32 swap;
1217+
int32 a1, a2;
1218+
1219+
/*
1220+
* Put the greater absolute value in arg1.
1221+
*
1222+
* This would happen automatically in the loop below, but avoids an
1223+
* expensive modulo operation, and simplifies the special-case handling
1224+
* for INT_MIN below.
1225+
*
1226+
* We do this in negative space in order to handle INT_MIN.
1227+
*/
1228+
a1 = (arg1 < 0) ? arg1 : -arg1;
1229+
a2 = (arg2 < 0) ? arg2 : -arg2;
1230+
if (a1 > a2)
1231+
{
1232+
swap = arg1;
1233+
arg1 = arg2;
1234+
arg2 = swap;
1235+
}
1236+
1237+
/* Special care needs to be taken with INT_MIN. See comments above. */
1238+
if (arg1 == PG_INT32_MIN)
1239+
{
1240+
if (arg2 == 0 || arg2 == PG_INT32_MIN)
1241+
ereport(ERROR,
1242+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1243+
errmsg("integer out of range")));
1244+
1245+
/*
1246+
* Some machines throw a floating-point exception for INT_MIN % -1,
1247+
* which is a bit silly since the correct answer is perfectly
1248+
* well-defined, namely zero. Guard against this and just return the
1249+
* result, gcd(INT_MIN, -1) = 1.
1250+
*/
1251+
if (arg2 == -1)
1252+
return 1;
1253+
}
1254+
1255+
/* Use the Euclidean algorithm to find the GCD */
1256+
while (arg2 != 0)
1257+
{
1258+
swap = arg2;
1259+
arg2 = arg1 % arg2;
1260+
arg1 = swap;
1261+
}
1262+
1263+
/*
1264+
* Make sure the result is positive. (We know we don't have INT_MIN
1265+
* anymore).
1266+
*/
1267+
if (arg1 < 0)
1268+
arg1 = -arg1;
1269+
1270+
return arg1;
1271+
}
1272+
1273+
Datum
1274+
int4gcd(PG_FUNCTION_ARGS)
1275+
{
1276+
int32 arg1 = PG_GETARG_INT32(0);
1277+
int32 arg2 = PG_GETARG_INT32(1);
1278+
int32 result;
1279+
1280+
result = int4gcd_internal(arg1, arg2);
1281+
1282+
PG_RETURN_INT32(result);
1283+
}
1284+
1285+
/*
1286+
* Least Common Multiple
1287+
*/
1288+
Datum
1289+
int4lcm(PG_FUNCTION_ARGS)
1290+
{
1291+
int32 arg1 = PG_GETARG_INT32(0);
1292+
int32 arg2 = PG_GETARG_INT32(1);
1293+
int32 gcd;
1294+
int32 result;
1295+
1296+
/*
1297+
* Handle lcm(x, 0) = lcm(0, x) = 0 as a special case. This prevents a
1298+
* division-by-zero error below when x is zero, and an overflow error from
1299+
* the GCD computation when x = INT_MIN.
1300+
*/
1301+
if (arg1 == 0 || arg2 == 0)
1302+
PG_RETURN_INT32(0);
1303+
1304+
/* lcm(x, y) = abs(x / gcd(x, y) * y) */
1305+
gcd = int4gcd_internal(arg1, arg2);
1306+
arg1 = arg1 / gcd;
1307+
1308+
if (unlikely(pg_mul_s32_overflow(arg1, arg2, &result)))
1309+
ereport(ERROR,
1310+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1311+
errmsg("integer out of range")));
1312+
1313+
/* If the result is INT_MIN, it cannot be represented. */
1314+
if (unlikely(result == PG_INT32_MIN))
1315+
ereport(ERROR,
1316+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
1317+
errmsg("integer out of range")));
1318+
1319+
if (result < 0)
1320+
result = -result;
1321+
1322+
PG_RETURN_INT32(result);
1323+
}
1324+
11991325
Datum
12001326
int2larger(PG_FUNCTION_ARGS)
12011327
{

‎src/backend/utils/adt/int8.c

+126
Original file line numberDiff line numberDiff line change
@@ -667,6 +667,132 @@ int8mod(PG_FUNCTION_ARGS)
667667
PG_RETURN_INT64(arg1 % arg2);
668668
}
669669

670+
/*
671+
* Greatest Common Divisor
672+
*
673+
* Returns the largest positive integer that exactly divides both inputs.
674+
* Special cases:
675+
* - gcd(x, 0) = gcd(0, x) = abs(x)
676+
* because 0 is divisible by anything
677+
* - gcd(0, 0) = 0
678+
* complies with the previous definition and is a common convention
679+
*
680+
* Special care must be taken if either input is INT64_MIN ---
681+
* gcd(0, INT64_MIN), gcd(INT64_MIN, 0) and gcd(INT64_MIN, INT64_MIN) are
682+
* all equal to abs(INT64_MIN), which cannot be represented as a 64-bit signed
683+
* integer.
684+
*/
685+
static int64
686+
int8gcd_internal(int64 arg1, int64 arg2)
687+
{
688+
int64 swap;
689+
int64 a1, a2;
690+
691+
/*
692+
* Put the greater absolute value in arg1.
693+
*
694+
* This would happen automatically in the loop below, but avoids an
695+
* expensive modulo operation, and simplifies the special-case handling
696+
* for INT64_MIN below.
697+
*
698+
* We do this in negative space in order to handle INT64_MIN.
699+
*/
700+
a1 = (arg1 < 0) ? arg1 : -arg1;
701+
a2 = (arg2 < 0) ? arg2 : -arg2;
702+
if (a1 > a2)
703+
{
704+
swap = arg1;
705+
arg1 = arg2;
706+
arg2 = swap;
707+
}
708+
709+
/* Special care needs to be taken with INT64_MIN. See comments above. */
710+
if (arg1 == PG_INT64_MIN)
711+
{
712+
if (arg2 == 0 || arg2 == PG_INT64_MIN)
713+
ereport(ERROR,
714+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
715+
errmsg("bigint out of range")));
716+
717+
/*
718+
* Some machines throw a floating-point exception for INT64_MIN % -1,
719+
* which is a bit silly since the correct answer is perfectly
720+
* well-defined, namely zero. Guard against this and just return the
721+
* result, gcd(INT64_MIN, -1) = 1.
722+
*/
723+
if (arg2 == -1)
724+
return 1;
725+
}
726+
727+
/* Use the Euclidean algorithm to find the GCD */
728+
while (arg2 != 0)
729+
{
730+
swap = arg2;
731+
arg2 = arg1 % arg2;
732+
arg1 = swap;
733+
}
734+
735+
/*
736+
* Make sure the result is positive. (We know we don't have INT64_MIN
737+
* anymore).
738+
*/
739+
if (arg1 < 0)
740+
arg1 = -arg1;
741+
742+
return arg1;
743+
}
744+
745+
Datum
746+
int8gcd(PG_FUNCTION_ARGS)
747+
{
748+
int64 arg1 = PG_GETARG_INT64(0);
749+
int64 arg2 = PG_GETARG_INT64(1);
750+
int64 result;
751+
752+
result = int8gcd_internal(arg1, arg2);
753+
754+
PG_RETURN_INT64(result);
755+
}
756+
757+
/*
758+
* Least Common Multiple
759+
*/
760+
Datum
761+
int8lcm(PG_FUNCTION_ARGS)
762+
{
763+
int64 arg1 = PG_GETARG_INT64(0);
764+
int64 arg2 = PG_GETARG_INT64(1);
765+
int64 gcd;
766+
int64 result;
767+
768+
/*
769+
* Handle lcm(x, 0) = lcm(0, x) = 0 as a special case. This prevents a
770+
* division-by-zero error below when x is zero, and an overflow error from
771+
* the GCD computation when x = INT64_MIN.
772+
*/
773+
if (arg1 == 0 || arg2 == 0)
774+
PG_RETURN_INT64(0);
775+
776+
/* lcm(x, y) = abs(x / gcd(x, y) * y) */
777+
gcd = int8gcd_internal(arg1, arg2);
778+
arg1 = arg1 / gcd;
779+
780+
if (unlikely(pg_mul_s64_overflow(arg1, arg2, &result)))
781+
ereport(ERROR,
782+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
783+
errmsg("bigint out of range")));
784+
785+
/* If the result is INT64_MIN, it cannot be represented. */
786+
if (unlikely(result == PG_INT64_MIN))
787+
ereport(ERROR,
788+
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
789+
errmsg("bigint out of range")));
790+
791+
if (result < 0)
792+
result = -result;
793+
794+
PG_RETURN_INT64(result);
795+
}
670796

671797
Datum
672798
int8inc(PG_FUNCTION_ARGS)

0 commit comments

Comments
 (0)