Subject: Re: Greatest Common Divisor
To: Fabien COELHO <coelho@cri.ensmp.fr>
Cc: PostgreSQL Hackers <pgsql-hackers@lists.postgresql.org>
References: <adbd3e0b-e3f1-5bbc-21db-03caf1cef0f7@2ndquadrant.com>
 <alpine.DEB.2.21.1912281848540.24861@pseudo>
 <dbf86808-91fa-3262-f939-2ef3d38ec242@2ndquadrant.com>
 <alpine.DEB.2.21.1912290815520.889@pseudo>
From: Vik Fearing <vik.fearing@2ndquadrant.com>
Message-ID: <9890a02a-94d4-537c-9f4a-ff02a5414720@2ndquadrant.com>
Date: Sun, 29 Dec 2019 14:58:49 +0100
User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:68.0) Gecko/20100101
 Thunderbird/68.2.2
MIME-Version: 1.0
In-Reply-To: <alpine.DEB.2.21.1912290815520.889@pseudo>
Content-Type: multipart/mixed;
 boundary="------------A07F80366BA3A32A5E612C67"
Content-Language: en-US
Precedence: bulk

This is a multi-part message in MIME format.
--------------A07F80366BA3A32A5E612C67
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: 8bit

On 29/12/2019 08:30, Fabien COELHO wrote:
>
>>> I'm wondering what it should do on N, 0 and 0, 0. Raise an error?
>>> Return 0? Return 1? return N? There should be some logic and comments
>>> explaining it.
>>
>> Well, gcd(N, 0) is N, and gcd(0, 0) is 0, so I don't see an issue here?
>
> I think that there should be a comment.


Done.


>>> It does not seem that fixing the sign afterwards is the right thing to
>>> do. I'd rather turn both arguments positive before doing the descent.
>>
>> Why isn't it the right thing to do?
>
> Because I do not trust C modulo as I had a lot of problems with it? :-)
>
> If it works, but it should deserve a clear comment explaining why.


Surely such a comment should be on the mod functions and not in this patch.


>
>>> Which raises an issue with INT_MIN by the way, which has no positive:-(
>>
>> That's an argument against abs-ing the input values.  It's only an issue
>> with gcd(INT_MIN, INT_MIN) which currently returns INT_MIN.
>
> That should be an error instead, because -INT_MIN cannot be represented?


Why should it error?  Is INT_MIN not a valid divisor of INT_MIN?  I
added a comment instead.


>>> Also, the usual approach is to exchange args so that the largest is
>>> first, because there may be a software emulation if the hardware does
>>> not implement modulo. At least it was the case with some sparc
>>> processors 25 years ago:-)
>>
>> The args will exchange themselves.
>
> Yep, but after a possibly expensive software-emulated modulo operation?
>

I'll just trust you on this.  Swap added.


Thanks!

-- 

Vik


--------------A07F80366BA3A32A5E612C67
Content-Type: text/x-patch; charset=UTF-8;
 name="gcd.0003.patch"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
 filename="gcd.0003.patch"

diff --git a/doc/src/sgml/func.sgml b/doc/src/sgml/func.sgml
index 57a1539506..b2c663cd92 100644
--- a/doc/src/sgml/func.sgml
+++ b/doc/src/sgml/func.sgml
@@ -870,6 +870,19 @@
        <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>
diff --git a/src/backend/utils/adt/int.c b/src/backend/utils/adt/int.c
index 04825fc77d..61b50db71c 100644
--- a/src/backend/utils/adt/int.c
+++ b/src/backend/utils/adt/int.c
@@ -1196,6 +1196,76 @@ int2abs(PG_FUNCTION_ARGS)
 	PG_RETURN_INT16(result);
 }
 
+/* int[24]gcd()
+ * Greatest Common Divisor
+ *
+ * By definition, gcd(0, 0) = 0 and gcd(x, 0) = x.
+ */
+Datum
+int4gcd(PG_FUNCTION_ARGS)
+{
+	int32		arg1 = PG_GETARG_INT32(0);
+	int32		arg2 = PG_GETARG_INT32(1);
+	int32		swap;
+
+	/*
+	 * Put the greater value in arg1.
+	 * This would happen automatically in the loop below, but avoids an
+	 * expensive modulo simulation on some architectures.
+	 */
+	if (arg1 < arg2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	/* Find GCD using the basic Euclidean algorithm */
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	/*
+	 * Make sure the result is positive. (INT_MIN will still be negative
+	 * because there is no positive representation of it).
+	 */
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	PG_RETURN_INT32(arg1);
+}
+
+Datum
+int2gcd(PG_FUNCTION_ARGS)
+{
+	/* See int4gcd for commented version. */
+	int16		arg1 = PG_GETARG_INT16(0);
+	int16		arg2 = PG_GETARG_INT16(1);
+	int16		swap;
+
+	if (arg1 < arg2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	PG_RETURN_INT16(arg1);
+}
+
 Datum
 int2larger(PG_FUNCTION_ARGS)
 {
diff --git a/src/backend/utils/adt/int8.c b/src/backend/utils/adt/int8.c
index a40ae40dff..59be1534a4 100644
--- a/src/backend/utils/adt/int8.c
+++ b/src/backend/utils/adt/int8.c
@@ -667,6 +667,38 @@ int8mod(PG_FUNCTION_ARGS)
 	PG_RETURN_INT64(arg1 % arg2);
 }
 
+/* int8gcd()
+ * Greatest Common Divisor
+ *
+ * See int4gcd for commented version.
+ */
+Datum
+int8gcd(PG_FUNCTION_ARGS)
+{
+	int64		arg1 = PG_GETARG_INT64(0);
+	int64		arg2 = PG_GETARG_INT64(1);
+	int64		swap;
+
+	if (arg1 < arg2)
+	{
+		swap = arg1;
+		arg1 = arg2;
+		arg2 = swap;
+	}
+
+	while (arg2 != 0)
+	{
+		swap = arg2;
+		arg2 = arg1 % arg2;
+		arg1 = swap;
+	}
+
+	if (arg1 < 0)
+		arg1 = -arg1;
+
+	PG_RETURN_INT64(arg1);
+}
+
 
 Datum
 int8inc(PG_FUNCTION_ARGS)
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
index ac8f64b219..b3057f2cdd 100644
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -10729,4 +10729,15 @@
   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' },
+
 ]
diff --git a/src/test/regress/expected/int2.out b/src/test/regress/expected/int2.out
index 8c255b9e4d..cca47af7fb 100644
--- a/src/test/regress/expected/int2.out
+++ b/src/test/regress/expected/int2.out
@@ -306,3 +306,23 @@ 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, (-32768)::int2);
+  gcd   
+--------
+ -32768
+(1 row)
+
diff --git a/src/test/regress/expected/int4.out b/src/test/regress/expected/int4.out
index bda7a8daef..3f0f37ffc7 100644
--- a/src/test/regress/expected/int4.out
+++ b/src/test/regress/expected/int4.out
@@ -403,3 +403,23 @@ 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, (-2147483648)::int4);
+     gcd     
+-------------
+ -2147483648
+(1 row)
+
diff --git a/src/test/regress/expected/int8.out b/src/test/regress/expected/int8.out
index 8447a28c3d..d6d1a76b29 100644
--- a/src/test/regress/expected/int8.out
+++ b/src/test/regress/expected/int8.out
@@ -886,3 +886,23 @@ 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, (-9223372036854775808)::int8);
+         gcd          
+----------------------
+ -9223372036854775808
+(1 row)
+
diff --git a/src/test/regress/sql/int2.sql b/src/test/regress/sql/int2.sql
index 7dbafb6dac..656b200e5e 100644
--- a/src/test/regress/sql/int2.sql
+++ b/src/test/regress/sql/int2.sql
@@ -112,3 +112,10 @@ 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, (-32768)::int2);
diff --git a/src/test/regress/sql/int4.sql b/src/test/regress/sql/int4.sql
index f014cb2d32..1faea33945 100644
--- a/src/test/regress/sql/int4.sql
+++ b/src/test/regress/sql/int4.sql
@@ -155,3 +155,10 @@ 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, (-2147483648)::int4);
diff --git a/src/test/regress/sql/int8.sql b/src/test/regress/sql/int8.sql
index e890452236..af9501c415 100644
--- a/src/test/regress/sql/int8.sql
+++ b/src/test/regress/sql/int8.sql
@@ -225,3 +225,10 @@ 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, (-9223372036854775808)::int8);

--------------A07F80366BA3A32A5E612C67--