Fortran: speed up decimal output of integers

libgfortran/ChangeLog:

	PR libfortran/98076
	* runtime/string.c (itoa64, itoa64_pad19): New helper functions.
	(gfc_itoa): On targets with 128-bit integers, call fast
	64-bit functions to avoid many slow divisions.

gcc/testsuite/ChangeLog:

	PR libfortran/98076
	* gfortran.dg/pr98076.f90: New test.

1 files changed, 60 insertions, 5 deletions

diff --git a/libgfortran/runtime/string.c b/libgfortran/runtime/string.c
index 835027a..0ccd731 100644
--- a/libgfortran/runtime/string.c
+++ b/libgfortran/runtime/string.c
@@ -23,6 +23,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 <http://www.gnu.org/licenses/>.  */
 
 #include "libgfortran.h"
+#include <assert.h>
 #include <string.h>
 #include <strings.h>
 
@@ -169,6 +170,38 @@ find_option (st_parameter_common *cmp, const char *s1, gfc_charlen_type s1_len,
 }
 
 
+/* Fast helper function for a positive value that fits in uint64_t.  */
+
+static inline char *
+itoa64 (uint64_t n, char *p)
+{
+  while (n != 0)
+    {
+      *--p = '0' + (n % 10);
+      n /= 10;
+    }
+  return p;
+}
+
+
+#if defined(HAVE_GFC_INTEGER_16)
+# define TEN19 ((GFC_UINTEGER_LARGEST) 1000000 * (GFC_UINTEGER_LARGEST) 1000000 * (GFC_UINTEGER_LARGEST) 10000000)
+
+/* Same as itoa64(), with zero padding of 19 digits.  */
+
+static inline char *
+itoa64_pad19 (uint64_t n, char *p)
+{
+  for (int k = 0; k < 19; k++)
+    {
+      *--p = '0' + (n % 10);
+      n /= 10;
+    }
+  return p;
+}
+#endif
+
+
 /* Integer to decimal conversion.
 
    This function is much more restricted than the widespread (but
@@ -195,11 +228,33 @@ gfc_itoa (GFC_UINTEGER_LARGEST n, char *buffer, size_t len)
   p = buffer + GFC_ITOA_BUF_SIZE - 1;
   *p = '\0';
 
-  while (n != 0)
+#if defined(HAVE_GFC_INTEGER_16)
+  /* On targets that have a 128-bit integer type, division in that type
+     is slow, because it occurs through a function call. We avoid that.  */
+
+  if (n <= UINT64_MAX)
+    /* If the value fits in uint64_t, use the fast function. */
+    return itoa64 (n, p);
+  else
     {
-      *--p = '0' + (n % 10);
-      n /= 10;
+      /* Otherwise, break down into smaller bits by division. Two calls to
+	 the uint64_t function are not sufficient for all 128-bit unsigned
+	 integers (we would need three calls), but they do suffice for all
+	 values up to 2^127, which is the largest that Fortran can produce
+	 (-HUGE(0_16)-1) with its signed integer types.  */
+      static_assert(sizeof(GFC_UINTEGER_LARGEST) <= 2 * sizeof(uint64_t));
+
+      GFC_UINTEGER_LARGEST r;
+      r = n % TEN19;
+      n = n / TEN19;
+      assert (r <= UINT64_MAX);
+      p = itoa64_pad19 (r, p);
+
+      assert(n <= UINT64_MAX);
+      return itoa64 (n, p);
     }
-
-  return p;
+#else
+  /* On targets where the largest integer is 64-bit, just use that.  */
+  return itoa64 (n, p);
+#endif
 }