Schubfach algorithm

README.md

@@ -12,9 +12,9 @@
 
 **dconv** is a header-only C++20 library for printing and parsing floating point numbers.
 
-Double to string conversion uses the **Grisu2** algorithm by default, described by **Florian Loitsch** in its publication [Printing Floating-Point Numbers Quickly and Accurately with Integers](https://florian.loitsch.com/publications).
+Double to string conversion uses the **Schubfach** algorithm by default, described by **Raffaello Giulietti** in its publication *The Schubfach way to render doubles*. It always prints the shortest representation that reads back to the same double, and the closest one when several are possible.
 
-The **Grisu** and **Grisu3** algorithms from the same publication are also available. Grisu prints up to 18 digits without looking for the shortest output, Grisu2 prints the shortest output for about 99.9% of the doubles, and Grisu3 detects the remaining cases and falls back to `std::to_chars`.
+The **Grisu**, **Grisu2** and **Grisu3** algorithms, described by **Florian Loitsch** in [Printing Floating-Point Numbers Quickly and Accurately with Integers](https://florian.loitsch.com/publications), are also available. Grisu prints up to 18 digits without looking for the shortest output, Grisu2 prints the shortest output for about 99.9% of the doubles, and Grisu3 detects the remaining cases and falls back to `std::to_chars`.
 
 String to double conversion uses the **Eisel-Lemire** algorithm on the first 19 significant digits. Longer mantissas are bracketed between two Eisel-Lemire results; only when they disagree is the number rebuilt from the digits already read and passed to strtod. A stream is read once and never rewound; a buffer is only read again for mantissas longer than 19 digits.
 
@@ -85,6 +85,7 @@ char* end = dconv::grisu3::dtoa (value, -2.22507e-308);
 
 | Header | Function |
 | :--- | :--- |
+| `dconv/schubfach.hpp` | `dconv::schubfach::dtoa` |
 | `dconv/grisu3.hpp` | `dconv::grisu3::dtoa` |
 | `dconv/grisu2.hpp` | `dconv::grisu2::dtoa` |
 | `dconv/grisu.hpp` | `dconv::grisu::dtoa` |

VERSION

@@ -1 +1 @@
-1.0.5
+1.1.0

dconv/CMakeLists.txt

@@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.22.1)
 set(PUBLIC_HEADERS
     include/dconv/error.hpp
     include/dconv/probe.hpp
+    include/dconv/power.hpp
     include/dconv/eisellemirepow.hpp
     include/dconv/eisellemire.hpp
     include/dconv/atod.hpp
@@ -13,6 +14,8 @@ set(PUBLIC_HEADERS
     include/dconv/grisu.hpp
     include/dconv/grisu2.hpp
     include/dconv/grisu3.hpp
+    include/dconv/schubfachpow.hpp
+    include/dconv/schubfach.hpp
     include/dconv/dtoa.hpp
 )
 

dconv/include/dconv/digits.hpp

@@ -40,6 +40,22 @@ namespace dconv
             "0001020304050607080910111213141516171819202122232425262728293031323334353637383940414243444546474849505152"
             "5354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899";
 
+        /// trailing zeros of the two digits decimal representations of 0 to 99.
+        alignas (64) constexpr uint8_t pairZeros[] = {
+            2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
+            0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
+            0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+        /**
+         * @brief count the trailing zeros of a four digits group.
+         * @param group four digits group.
+         * @return number of trailing zeros.
+         */
+        inline int countTrailingZeros (uint32_t group) noexcept
+        {
+            return pairZeros[group % 100] + ((group % 100 == 0) ? pairZeros[group / 100] : 0);
+        }
+
         /**
          * @brief count the digits of a decimal significand.
          * @param digits decimal significand.
@@ -82,6 +98,86 @@ namespace dconv
             return (digits >= 1000000000) ? 10 : 9;
         }
 
+        /**
+         * @brief write the digits of a decimal significand.
+         * @param buffer buffer to write the digits to.
+         * @param digits decimal significand.
+         * @param length number of digits.
+         * @return number of trailing zeros written.
+         */
+        inline int writeDigits (char* buffer, uint64_t digits, int length) noexcept
+        {
+            uint32_t high = static_cast<uint32_t> (digits);
+            uint32_t low = 0;
+            if (digits >= 100000000)
+            {
+                high = static_cast<uint32_t> (digits / 100000000);
+                low = static_cast<uint32_t> (digits % 100000000);
+            }
+
+            char* const end = buffer + length;
+            char* cur = end;
+            int zeros = 0;
+            if (digits >= 100000000)
+            {
+                cur -= 8;
+                if (low == 0)
+                {
+                    __builtin_memcpy (cur, "00000000", 8);
+                    zeros = 8;
+                }
+                else
+                {
+                    const uint32_t upper = low / 10000;
+                    const uint32_t lower = low % 10000;
+                    __builtin_memcpy (cur, &digitPairs[(upper / 100) * 2], 2);
+                    __builtin_memcpy (cur + 2, &digitPairs[(upper % 100) * 2], 2);
+                    __builtin_memcpy (cur + 4, &digitPairs[(lower / 100) * 2], 2);
+                    __builtin_memcpy (cur + 6, &digitPairs[(lower % 100) * 2], 2);
+                    zeros = (lower != 0) ? countTrailingZeros (lower) : 4 + countTrailingZeros (upper);
+                }
+            }
+
+            if (high >= 10000)
+            {
+                const uint32_t lower = high % 10000;
+                high /= 10000;
+                cur -= 4;
+                __builtin_memcpy (cur, &digitPairs[(lower / 100) * 2], 2);
+                __builtin_memcpy (cur + 2, &digitPairs[(lower % 100) * 2], 2);
+                if (zeros == end - cur - 4)
+                {
+                    zeros += countTrailingZeros (lower);
+                }
+            }
+
+            while (high >= 100)
+            {
+                cur -= 2;
+                __builtin_memcpy (cur, &digitPairs[(high % 100) * 2], 2);
+                if (zeros == end - cur - 2)
+                {
+                    zeros += pairZeros[high % 100];
+                }
+                high /= 100;
+            }
+
+            if (high >= 10)
+            {
+                cur -= 2;
+                __builtin_memcpy (cur, &digitPairs[high * 2], 2);
+                if (zeros == end - cur - 2)
+                {
+                    zeros += pairZeros[high];
+                }
+            }
+            else
+            {
+                *--cur = static_cast<char> ('0' + high);
+            }
+
+            return zeros;
+        }
     }
 }
 

dconv/include/dconv/dtoa.hpp

@@ -26,6 +26,7 @@
 #define __DCONV_DTOA_HPP__
 
 // dconv.
+#include <dconv/schubfach.hpp>
 #include <dconv/grisu3.hpp>
 
 namespace dconv
@@ -38,7 +39,7 @@ namespace dconv
      */
     inline char* dtoa (char* buffer, double value) noexcept
     {
-        return grisu2::dtoa (buffer, value);
+        return schubfach::dtoa (buffer, value);
     }
 }
 

dconv/include/dconv/eisellemire.hpp

@@ -40,34 +40,6 @@ namespace dconv
 {
     namespace eisellemire
     {
-        /**
-         * @brief multiply two 64 bits unsigned integers.
-         * @param a multiplicand.
-         * @param b multiplier.
-         * @param high high 64 bits of the 128-bit result.
-         * @param low low 64 bits of the 128-bit result.
-         */
-        inline void umul128 (uint64_t a, uint64_t b, uint64_t& high, uint64_t& low) noexcept
-        {
-#if defined(__SIZEOF_INT128__)
-            __uint128_t p = static_cast<__uint128_t> (a) * b;
-
-            high = static_cast<uint64_t> (p >> 64);
-            low = static_cast<uint64_t> (p);
-#else
-            uint64_t b_lo = static_cast<uint32_t> (b);
-            uint64_t b_hi = b >> 32;
-            uint64_t p0 = (a & 0xFFFFFFFF) * b_lo;
-            uint64_t p1 = (a >> 32) * b_lo;
-            uint64_t p2 = (a & 0xFFFFFFFF) * b_hi;
-            uint64_t p3 = (a >> 32) * b_hi;
-            uint64_t carry = (p0 >> 32) + (p1 & 0xFFFFFFFF) + (p2 & 0xFFFFFFFF);
-
-            low = (carry << 32) | (p0 & 0xFFFFFFFF);
-            high = (carry >> 32) + (p1 >> 32) + (p2 >> 32) + p3;
-#endif
-        }
-
         /**
          * @brief convert double using fast path.
          * @param significand significand digits (at most 19).
@@ -110,11 +82,11 @@ namespace dconv
 
             const details::Power& power = details::eisellemirepow[exponent + 342];
             uint64_t high, low;
-            umul128 (power.hi, significand, high, low);
+            details::umul128 (power.hi, significand, high, low);
             if ((high & 0x1FF) == 0x1FF)
             {
                 uint64_t secondHigh, secondLow;
-                umul128 (power.lo, significand, secondHigh, secondLow);
+                details::umul128 (power.lo, significand, secondHigh, secondLow);
                 low += secondHigh;
                 high += (secondHigh > low);
             }

dconv/include/dconv/eisellemirepow.hpp

@@ -25,18 +25,13 @@
 #ifndef __DCONV_EISELLEMIREPOW_HPP__
 #define __DCONV_EISELLEMIREPOW_HPP__
 
-#include <cstdint>
+// dconv.
+#include <dconv/power.hpp>
 
 namespace dconv
 {
     namespace details
     {
-        struct Power
-        {
-            uint64_t hi;
-            uint64_t lo;
-        };
-
         constexpr Power eisellemirepow[] = {
             {0xeef453d6923bd65a, 0x113faa2906a13b3f}, {0x9558b4661b6565f8, 0x4ac7ca59a424c507},
             {0xbaaee17fa23ebf76, 0x5d79bcf00d2df649}, {0xe95a99df8ace6f53, 0xf4d82c2c107973dc},

dconv/include/dconv/format.hpp

@@ -131,6 +131,63 @@ namespace dconv
                 return writeExponent (&buffer[length + 2], kk - 1);
             }
         }
+
+        /**
+         * @brief write a decimal significand with the same layout as prettify.
+         * @param buffer buffer to write the string representation to.
+         * @param digits decimal significand.
+         * @param length number of digits.
+         * @param k decimal exponent.
+         * @return end position.
+         */
+        inline char* format (char* buffer, uint64_t digits, int length, int k) noexcept
+        {
+            const int kk = length + k;
+
+            if ((kk > 21) || (kk <= -6))
+            {
+                length -= writeDigits (buffer + 1, digits, length);
+                buffer[0] = buffer[1];
+                buffer[1] = '.';
+                const int position = (length > 1) ? length + 1 : 1;
+                buffer[position] = 'e';
+                return writeExponent (&buffer[position + 1], kk - 1);
+            }
+
+            if (length <= kk)
+            {
+                writeDigits (buffer, digits, length);
+                for (int i = length; i < kk; ++i)
+                {
+                    buffer[i] = '0';
+                }
+                buffer[kk] = '.';
+                buffer[kk + 1] = '0';
+                return &buffer[kk + 2];
+            }
+            else if (0 < kk)
+            {
+                length -= writeDigits (buffer + 1, digits, length);
+                for (int i = 0; i < kk; ++i)
+                {
+                    buffer[i] = buffer[i + 1];
+                }
+                buffer[kk] = '.';
+                return &buffer[length + 1];
+            }
+            else
+            {
+                const int offset = 2 - kk;
+                buffer[0] = '0';
+                buffer[1] = '.';
+                for (int i = 2; i < offset; ++i)
+                {
+                    buffer[i] = '0';
+                }
+                length -= writeDigits (buffer + offset, digits, length);
+                return &buffer[length + offset];
+            }
+        }
     }
 }
 

dconv/include/dconv/power.hpp

@@ -0,0 +1,75 @@
+/**
+ * MIT License
+ *
+ * Copyright (c) 2026 Mathieu Rabine
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#ifndef __DCONV_POWER_HPP__
+#define __DCONV_POWER_HPP__
+
+// C.
+#include <cstdint>
+
+namespace dconv
+{
+    namespace details
+    {
+        /// 128 bits power.
+        struct Power
+        {
+            /// high 64 bits.
+            uint64_t hi;
+
+            /// low 64 bits.
+            uint64_t lo;
+        };
+
+        /**
+         * @brief multiply two 64 bits unsigned integers.
+         * @param a multiplicand.
+         * @param b multiplier.
+         * @param high high 64 bits of the 128-bit result.
+         * @param low low 64 bits of the 128-bit result.
+         */
+        inline void umul128 (uint64_t a, uint64_t b, uint64_t& high, uint64_t& low) noexcept
+        {
+#if defined(__SIZEOF_INT128__)
+            __uint128_t p = static_cast<__uint128_t> (a) * b;
+
+            high = static_cast<uint64_t> (p >> 64);
+            low = static_cast<uint64_t> (p);
+#else
+            uint64_t b_lo = static_cast<uint32_t> (b);
+            uint64_t b_hi = b >> 32;
+            uint64_t p0 = (a & 0xFFFFFFFF) * b_lo;
+            uint64_t p1 = (a >> 32) * b_lo;
+            uint64_t p2 = (a & 0xFFFFFFFF) * b_hi;
+            uint64_t p3 = (a >> 32) * b_hi;
+            uint64_t carry = (p0 >> 32) + (p1 & 0xFFFFFFFF) + (p2 & 0xFFFFFFFF);
+
+            low = (carry << 32) | (p0 & 0xFFFFFFFF);
+            high = (carry >> 32) + (p1 >> 32) + (p2 >> 32) + p3;
+#endif
+        }
+    }
+}
+
+#endif