Repository: numpy/numpy
Pull request: BUG, ENH: fix array2string rounding bug by adding min_digits option
Fixes #18609

This PR adds a new `min_digits` argument to the dragon4 float printing methods, such as `np.format_float_positional`. This kwd guarantees that at least the given number of digits will be printed, when printing in `unique=True` mode, even if the extra digits are unnecessary to uniquely specify the value. It is the counterpart to the `precision` argument which sets the maximum number of digits to be printed. (When `unique=False` fixed-precision mode, it has no effect, and the `precision` arg fixes the number of digits)

Using this new keyword, I then fix a bug in arrayprint ( #18609). The bug in the old code is that, in order to guarantee that all elements of an array print using the same number of digits, we turned off `unique=True` because we needed to use fixed-precision mode to fix the number of digits. However, unique mode rounds differently from fixed-precision mode, so the rounding was wrong. This PR instead now stays in unique mode to print the elements, but uses the `min_digits` and `precision` arguments to fix the number of output digits.

It looks like a lof of chages, but a lot of it is just introducing the `min_digits` arg into all the function signatures in dragon4 and in the user api.

Problem statement:
floatmode=unique sometimes giving wrong rounding
Let `a` be a special float64:
```python
In [26]: a = np.float64.fromhex('-1p-97')

In [27]: a
Out[27]: -6.310887241768095e-30
```

It has `uniquestra`, a shortest decimal string that uniquely identifies it (thus can be eval-ed to reproduce it exactly):
```python
In [33]: uniquestra = np.format_float_scientific(a, unique=True)

In [34]: uniquestra
Out[34]: '-6.310887241768095e-30'

In [35]: eval(uniquestra)==a
Out[35]: True
```

Because `a` is a special float, unlike most other floats, it has a correctly-rounded precision=15 decimal string _that is different from `uniquestra`_:
```python
In [36]: roundedstra = np.format_float_scientific(a, unique=False, precision=15)

In [37]: roundedstra
Out[37]: '-6.310887241768094e-30'

In [38]: eval(roundedstra)==a
Out[38]: False
```
So far, so good. This kind of thing is why `format_float_scientific` has the `unique` flag, to choose between these two cases depending on which is needed in a given situation.

The problem comes with `array2string`, who takes a `floatmode=unique` option that should give `uniquestr` and instead gives `roundedstr`. The problem is that internally `array2string` uses `FloatingFormat` which in `floatmode='unique'` mode effectively uses format_float_scientific twice: Once with `unique=True, precision=None` to figure out how many digits of precision are enough to uniquely represent the float, and then a second time with `unique=False, precision=foo` to actually do the formatting. Most of the time this correctly gives a unique representation of the float, but for a special float like `a` it fails and gives `roundedstra` instead of `uniquestra`. 

### Reproducing code example:


```python
import numpy as np
a = np.float64.fromhex('-1p-97')
assert eval(np.array2string(a, floatmode='unique'))==a
```


### NumPy/Python version information:
```
1.20.1 3.8.8 (default, Feb 24 2021, 13:46:16) 
[Clang 10.0.0 ]
```

Hints from issue discussion:
Looks like a bug, probably/possibly with `dragon4` itself?  Ping @ahaldane do you know where to fix this and if we should report it "upstream"?
I think @levbishop found the problem, it is in arrayprint/FloatingFormatter rather than in dragon4.

It happens because although we use unique mode to determine the element length, when actually printing the values we switch to non-unique mode to make sure all the floats printed in scientific mode in an array have the same length (i.e, same number of digits after the decimal point) , which evidently unintentionally causes a different kind of rounding to occur.

I remember thinking about this when we updated that code, but it seems the current code solution doesn't work for these cases.

Probably we need to add some kind of "digits_right" or "min_digits" option to `np.format_float_scientific` (and  np.format_float_positional for symmetry), which adds extra unnecessary digits when printing in unique mode. Easiest to implement would be to just add 0s to the end, but maybe more satisfying to try to add computed digits somehow. (Edit: or try to update the meaning of the "precision" argument so it can add digits in unique mode, not just remove them)

I'll take a deeper look later.
Ah, right. I thought it was weird that with identical number of digits the result depends on `unique`.  But it actually does, since `unique=True` has to find the representation actual closest, while without `unique` it has to do correct decimal rounding. Simple, but tricky :).  Probably repeating the initial report, but if anyone wonders, this sums up the fact nicely that ending in  ending in `94` is the correct decimal rounding, but is actually closer to the next smaller float:
```
In [18]: np.format_float_scientific(np.float64.fromhex('-1p-97'), 17, unique=False)
Out[18]: '-6.31088724176809444e-30'

In [19]: np.format_float_scientific(np.nextafter(np.float64.fromhex('-1p-97'), 0), 17, unique=False)
Out[19]: '-6.31088724176809374e-30'

In [20]: np.format_float_scientific(np.nextafter(np.float64.fromhex('-1p-97'), -1), 17, unique=False)
Out[20]: '-6.31088724176809584e-30'
```

Commit message: BUG/ENH: fix array2string rounding bug by adding min_digits option

Fixes #18609

Patch to review:
diff --git a/numpy/core/arrayprint.py b/numpy/core/arrayprint.py
index 5c1d6cb63c6c..523820c33c11 100644
--- a/numpy/core/arrayprint.py
+++ b/numpy/core/arrayprint.py
@@ -914,6 +914,7 @@ def fillFormat(self, data):
             self.trim = '.'
             self.exp_size = -1
             self.unique = True
+            self.min_digits = None
         elif self.exp_format:
             trim, unique = '.', True
             if self.floatmode == 'fixed' or self._legacy == '1.13':
@@ -927,6 +928,8 @@ def fillFormat(self, data):
 
             self.trim = 'k'
             self.precision = max(len(s) for s in frac_part)
+            self.min_digits = self.precision
+            self.unique = unique
 
             # for back-compat with np 1.13, use 2 spaces & sign and full prec
             if self._legacy == '1.13':
@@ -936,10 +939,7 @@ def fillFormat(self, data):
                 self.pad_left = max(len(s) for s in int_part)
             # pad_right is only needed for nan length calculation
             self.pad_right = self.exp_size + 2 + self.precision
-
-            self.unique = False
         else:
-            # first pass printing to determine sizes
             trim, unique = '.', True
             if self.floatmode == 'fixed':
                 trim, unique = 'k', False
@@ -955,14 +955,14 @@ def fillFormat(self, data):
                 self.pad_left = max(len(s) for s in int_part)
             self.pad_right = max(len(s) for s in frac_part)
             self.exp_size = -1
+            self.unique = unique
 
             if self.floatmode in ['fixed', 'maxprec_equal']:
-                self.precision = self.pad_right
-                self.unique = False
+                self.precision = self.min_digits = self.pad_right
                 self.trim = 'k'
             else:
-                self.unique = True
                 self.trim = '.'
+                self.min_digits = 0
 
         if self._legacy != '1.13':
             # account for sign = ' ' by adding one to pad_left
@@ -991,6 +991,7 @@ def __call__(self, x):
         if self.exp_format:
             return dragon4_scientific(x,
                                       precision=self.precision,
+                                      min_digits=self.min_digits,
                                       unique=self.unique,
                                       trim=self.trim,
                                       sign=self.sign == '+',
@@ -999,6 +1000,7 @@ def __call__(self, x):
         else:
             return dragon4_positional(x,
                                       precision=self.precision,
+                                      min_digits=self.min_digits,
                                       unique=self.unique,
                                       fractional=True,
                                       trim=self.trim,
@@ -1009,7 +1011,8 @@ def __call__(self, x):
 
 @set_module('numpy')
 def format_float_scientific(x, precision=None, unique=True, trim='k',
-                            sign=False, pad_left=None, exp_digits=None):
+                            sign=False, pad_left=None, exp_digits=None,
+                            min_digits=None):
     """
     Format a floating-point scalar as a decimal string in scientific notation.
 
@@ -1027,11 +1030,12 @@ def format_float_scientific(x, precision=None, unique=True, trim='k',
         If `True`, use a digit-generation strategy which gives the shortest
         representation which uniquely identifies the floating-point number from
         other values of the same type, by judicious rounding. If `precision`
-        was omitted, print all necessary digits, otherwise digit generation is
-        cut off after `precision` digits and the remaining value is rounded.
+        is given fewer digits than necessary can be printed, or if `min_digits`
+        is given more can be printed, in which cases the last digit is rounded
+        with unbiased rounding.
         If `False`, digits are generated as if printing an infinite-precision
         value and stopping after `precision` digits, rounding the remaining
-        value.
+        value with unbiased rounding
     trim : one of 'k', '.', '0', '-', optional
         Controls post-processing trimming of trailing digits, as follows:
 
@@ -1048,6 +1052,10 @@ def format_float_scientific(x, precision=None, unique=True, trim='k',
     exp_digits : non-negative integer, optional
         Pad the exponent with zeros until it contains at least this many digits.
         If omitted, the exponent will be at least 2 digits.
+    min_digits : non-negative integer or None, optional
+        Minimum number of digits to print. Only has an effect if `unique=True`
+        in which case additional digits past those necessary to uniquely
+        identify the value may be printed, rounding the last additional digit.
 
     Returns
     -------
@@ -1071,15 +1079,18 @@ def format_float_scientific(x, precision=None, unique=True, trim='k',
     precision = _none_or_positive_arg(precision, 'precision')
     pad_left = _none_or_positive_arg(pad_left, 'pad_left')
     exp_digits = _none_or_positive_arg(exp_digits, 'exp_digits')
+    min_digits = _none_or_positive_arg(min_digits, 'min_digits')
+    if min_digits > 0 and precision > 0 and min_digits > precision:
+        raise ValueError("min_digits must be less than or equal to precision")
     return dragon4_scientific(x, precision=precision, unique=unique,
                               trim=trim, sign=sign, pad_left=pad_left,
-                              exp_digits=exp_digits)
+                              exp_digits=exp_digits, min_digits=min_digits)
 
 
 @set_module('numpy')
 def format_float_positional(x, precision=None, unique=True,
                             fractional=True, trim='k', sign=False,
-                            pad_left=None, pad_right=None):
+                            pad_left=None, pad_right=None, min_digits=None):
     """
     Format a floating-point scalar as a decimal string in positional notation.
 
@@ -1097,16 +1108,19 @@ def format_float_positional(x, precision=None, unique=True,
         If `True`, use a digit-generation strategy which gives the shortest
         representation which uniquely identifies the floating-point number from
         other values of the same type, by judicious rounding. If `precision`
-        was omitted, print out all necessary digits, otherwise digit generation
-        is cut off after `precision` digits and the remaining value is rounded.
+        is given fewer digits than necessary can be printed, or if `min_digits`
+        is given more can be printed, in which cases the last digit is rounded
+        with unbiased rounding.
         If `False`, digits are generated as if printing an infinite-precision
         value and stopping after `precision` digits, rounding the remaining
-        value.
+        value with unbiased rounding
     fractional : boolean, optional
-        If `True`, the cutoff of `precision` digits refers to the total number
-        of digits after the decimal point, including leading zeros.
-        If `False`, `precision` refers to the total number of significant
-        digits, before or after the decimal point, ignoring leading zeros.
+        If `True`, the cutoffs of `precision` and `min_digits` refer to the
+        total number of digits after the decimal point, including leading
+        zeros.
+        If `False`, `precision` and `min_digits` refer to the total number of
+        significant digits, before or after the decimal point, ignoring leading
+        zeros.
     trim : one of 'k', '.', '0', '-', optional
         Controls post-processing trimming of trailing digits, as follows:
 
@@ -1123,6 +1137,10 @@ def format_float_positional(x, precision=None, unique=True,
     pad_right : non-negative integer, optional
         Pad the right side of the string with whitespace until at least that
         many characters are to the right of the decimal point.
+    min_digits : non-negative integer or None, optional
+        Minimum number of digits to print. Only has an effect if `unique=True`
+        in which case additional digits past those necessary to uniquely
+        identify the value may be printed, rounding the last additional digit.
 
     Returns
     -------
@@ -1147,10 +1165,16 @@ def format_float_positional(x, precision=None, unique=True,
     precision = _none_or_positive_arg(precision, 'precision')
     pad_left = _none_or_positive_arg(pad_left, 'pad_left')
     pad_right = _none_or_positive_arg(pad_right, 'pad_right')
+    min_digits = _none_or_positive_arg(min_digits, 'min_digits')
+    if not fractional and precision == 0:
+        raise ValueError("precision must be greater than 0 if "
+                         "fractional=False")
+    if min_digits > 0 and precision > 0 and min_digits > precision:
+        raise ValueError("min_digits must be less than or equal to precision")
     return dragon4_positional(x, precision=precision, unique=unique,
                               fractional=fractional, trim=trim,
                               sign=sign, pad_left=pad_left,
-                              pad_right=pad_right)
+                              pad_right=pad_right, min_digits=min_digits)
 
 
 class IntegerFormat:
diff --git a/numpy/core/arrayprint.pyi b/numpy/core/arrayprint.pyi
index d2a5fdef963f..ac2b6f5a8abb 100644
--- a/numpy/core/arrayprint.pyi
+++ b/numpy/core/arrayprint.pyi
@@ -103,6 +103,7 @@ def format_float_scientific(
     sign: bool = ...,
     pad_left: Optional[int] = ...,
     exp_digits: Optional[int] = ...,
+    min_digits: Optional[int] = ...,
 ) -> str: ...
 def format_float_positional(
     x: _FloatLike_co,
@@ -113,6 +114,7 @@ def format_float_positional(
     sign: bool = ...,
     pad_left: Optional[int] = ...,
     pad_right: Optional[int] = ...,
+    min_digits: Optional[int] = ...,
 ) -> str: ...
 def array_repr(
     arr: ndarray[Any, Any],
diff --git a/numpy/core/src/multiarray/dragon4.c b/numpy/core/src/multiarray/dragon4.c
index a7b252a77e7e..1d8c275700a2 100644
--- a/numpy/core/src/multiarray/dragon4.c
+++ b/numpy/core/src/multiarray/dragon4.c
@@ -1130,8 +1130,9 @@ BigInt_ShiftLeft(BigInt *result, npy_uint32 shift)
  *   * exponent - value exponent in base 2
  *   * mantissaBit - index of the highest set mantissa bit
  *   * hasUnequalMargins - is the high margin twice as large as the low margin
- *   * cutoffMode - how to interpret cutoffNumber: fractional or total digits?
- *   * cutoffNumber - cut off printing after this many digits. -1 for no cutoff
+ *   * cutoffMode - how to interpret cutoff_*: fractional or total digits?
+ *   * cutoff_max - cut off printing after this many digits. -1 for no cutoff
+ *   * cutoff_min - print at least this many digits. -1 for no cutoff
  *   * pOutBuffer - buffer to output into
  *   * bufferSize - maximum characters that can be printed to pOutBuffer
  *   * pOutExponent - the base 10 exponent of the first digit
@@ -1142,7 +1143,7 @@ static npy_uint32
 Dragon4(BigInt *bigints, const npy_int32 exponent,
         const npy_uint32 mantissaBit, const npy_bool hasUnequalMargins,
         const DigitMode digitMode, const CutoffMode cutoffMode,
-        npy_int32 cutoffNumber, char *pOutBuffer,
+        npy_int32 cutoff_max, npy_int32 cutoff_min, char *pOutBuffer,
         npy_uint32 bufferSize, npy_int32 *pOutExponent)
 {
     char *curDigit = pOutBuffer;
@@ -1169,7 +1170,8 @@ Dragon4(BigInt *bigints, const npy_int32 exponent,
     BigInt *temp2 = &bigints[6];
 
     const npy_float64 log10_2 = 0.30102999566398119521373889472449;
-    npy_int32 digitExponent, cutoffExponent, hiBlock;
+    npy_int32 digitExponent, hiBlock;
+    npy_int32 cutoff_max_Exponent, cutoff_min_Exponent;
     npy_uint32 outputDigit;    /* current digit being output */
     npy_uint32 outputLen;
     npy_bool isEven = BigInt_IsEven(mantissa);
@@ -1294,9 +1296,9 @@ Dragon4(BigInt *bigints, const npy_int32 exponent,
      * increases the number. This will either correct digitExponent to an
      * accurate value or it will clamp it above the accurate value.
      */
-    if (cutoffNumber >= 0 && cutoffMode == CutoffMode_FractionLength &&
-            digitExponent <= -cutoffNumber) {
-        digitExponent = -cutoffNumber + 1;
+    if (cutoff_max >= 0 && cutoffMode == CutoffMode_FractionLength &&
+            digitExponent <= -cutoff_max) {
+        digitExponent = -cutoff_max + 1;
     }
 
 
@@ -1347,26 +1349,44 @@ Dragon4(BigInt *bigints, const npy_int32 exponent,
     }
 
     /*
-     * Compute the cutoff exponent (the exponent of the final digit to print).
-     * Default to the maximum size of the output buffer.
+     * Compute the cutoff_max exponent (the exponent of the final digit to
+     * print).  Default to the maximum size of the output buffer.
      */
-    cutoffExponent = digitExponent - bufferSize;
-    if (cutoffNumber >= 0) {
+    cutoff_max_Exponent = digitExponent - bufferSize;
+    if (cutoff_max >= 0) {
         npy_int32 desiredCutoffExponent;
 
         if (cutoffMode == CutoffMode_TotalLength) {
-            desiredCutoffExponent = digitExponent - cutoffNumber;
-            if (desiredCutoffExponent > cutoffExponent) {
-                cutoffExponent = desiredCutoffExponent;
+            desiredCutoffExponent = digitExponent - cutoff_max;
+            if (desiredCutoffExponent > cutoff_max_Exponent) {
+                cutoff_max_Exponent = desiredCutoffExponent;
             }
         }
-        /* Otherwise it's CutoffMode_FractionLength. Print cutoffNumber digits
+        /* Otherwise it's CutoffMode_FractionLength. Print cutoff_max digits
          * past the decimal point or until we reach the buffer size
          */
         else {
-            desiredCutoffExponent = -cutoffNumber;
-            if (desiredCutoffExponent > cutoffExponent) {
-                cutoffExponent = desiredCutoffExponent;
+            desiredCutoffExponent = -cutoff_max;
+            if (desiredCutoffExponent > cutoff_max_Exponent) {
+                cutoff_max_Exponent = desiredCutoffExponent;
+            }
+        }
+    }
+    /* Also compute the cutoff_min exponent. */
+    cutoff_min_Exponent = digitExponent;
+    if (cutoff_min >= 0) {
+        npy_int32 desiredCutoffExponent;
+
+        if (cutoffMode == CutoffMode_TotalLength) {
+            desiredCutoffExponent = digitExponent - cutoff_min;
+            if (desiredCutoffExponent < cutoff_min_Exponent) {
+                cutoff_min_Exponent = desiredCutoffExponent;
+            }
+        }
+        else {
+            desiredCutoffExponent = -cutoff_min;
+            if (desiredCutoffExponent < cutoff_min_Exponent) {
+                cutoff_min_Exponent = desiredCutoffExponent;
             }
         }
     }
@@ -1432,14 +1452,17 @@ Dragon4(BigInt *bigints, const npy_int32 exponent,
 
             /*
              * stop looping if we are far enough away from our neighboring
-             * values or if we have reached the cutoff digit
+             * values (and we have printed at least the requested minimum
+             * digits) or if we have reached the cutoff digit
              */
             cmp = BigInt_Compare(scaledValue, scaledMarginLow);
             low = isEven ? (cmp <= 0) : (cmp < 0);
             cmp = BigInt_Compare(scaledValueHigh, scale);
             high = isEven ? (cmp >= 0) : (cmp > 0);
-            if (low | high | (digitExponent == cutoffExponent))
+            if (((low | high) & (digitExponent <= cutoff_min_Exponent)) |
+                    (digitExponent == cutoff_max_Exponent)) {
                 break;
+            }
 
             /* store the output digit */
             *curDigit = (char)('0' + outputDigit);
@@ -1471,7 +1494,7 @@ Dragon4(BigInt *bigints, const npy_int32 exponent,
             DEBUG_ASSERT(outputDigit < 10);
 
             if ((scaledValue->length == 0) |
-                    (digitExponent == cutoffExponent)) {
+                    (digitExponent == cutoff_max_Exponent)) {
                 break;
             }
 
@@ -1589,6 +1612,7 @@ typedef struct Dragon4_Options {
     DigitMode digit_mode;
     CutoffMode cutoff_mode;
     npy_int32 precision;
+    npy_int32 min_digits;
     npy_bool sign;
     TrimMode trim_mode;
     npy_int32 digits_left;
@@ -1617,11 +1641,12 @@ FormatPositional(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
                  npy_int32 exponent, char signbit, npy_uint32 mantissaBit,
                  npy_bool hasUnequalMargins, DigitMode digit_mode,
                  CutoffMode cutoff_mode, npy_int32 precision,
-                 TrimMode trim_mode, npy_int32 digits_left,
-                 npy_int32 digits_right)
+                 npy_int32 min_digits, TrimMode trim_mode,
+                 npy_int32 digits_left, npy_int32 digits_right)
 {
     npy_int32 printExponent;
     npy_int32 numDigits, numWholeDigits=0, has_sign=0;
+    npy_int32 add_digits;
 
     npy_int32 maxPrintLen = (npy_int32)bufferSize - 1, pos = 0;
 
@@ -1644,8 +1669,9 @@ FormatPositional(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
     }
 
     numDigits = Dragon4(mantissa, exponent, mantissaBit, hasUnequalMargins,
-                        digit_mode, cutoff_mode, precision, buffer + has_sign,
-                        maxPrintLen - has_sign, &printExponent);
+                        digit_mode, cutoff_mode, precision, min_digits,
+                        buffer + has_sign, maxPrintLen - has_sign,
+                        &printExponent);
 
     DEBUG_ASSERT(numDigits > 0);
     DEBUG_ASSERT(numDigits <= bufferSize);
@@ -1744,9 +1770,10 @@ FormatPositional(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
         buffer[pos++] = '.';
     }
 
-    desiredFractionalDigits = precision;
-    if (cutoff_mode == CutoffMode_TotalLength && precision >= 0) {
-        desiredFractionalDigits = precision - numWholeDigits;
+    add_digits = digit_mode == DigitMode_Unique ? min_digits : precision;
+    desiredFractionalDigits = add_digits < 0 ? 0 : add_digits;
+    if (cutoff_mode == CutoffMode_TotalLength) {
+        desiredFractionalDigits = add_digits - numWholeDigits;
     }
 
     if (trim_mode == TrimMode_LeaveOneZero) {
@@ -1757,10 +1784,9 @@ FormatPositional(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
         }
     }
     else if (trim_mode == TrimMode_None &&
-             digit_mode != DigitMode_Unique &&
              desiredFractionalDigits > numFractionDigits &&
              pos < maxPrintLen) {
-        /* add trailing zeros up to precision length */
+        /* add trailing zeros up to add_digits length */
         /* compute the number of trailing zeros needed */
         npy_int32 count = desiredFractionalDigits - numFractionDigits;
         if (pos + count > maxPrintLen) {
@@ -1778,7 +1804,7 @@ FormatPositional(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
      * when rounding, we may still end up with trailing zeros. Remove them
      * depending on trim settings.
      */
-    if (precision >= 0 && trim_mode != TrimMode_None && numFractionDigits > 0) {
+    if (trim_mode != TrimMode_None && numFractionDigits > 0) {
         while (buffer[pos-1] == '0') {
             pos--;
             numFractionDigits--;
@@ -1852,7 +1878,7 @@ static npy_uint32
 FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
                   npy_int32 exponent, char signbit, npy_uint32 mantissaBit,
                   npy_bool hasUnequalMargins, DigitMode digit_mode,
-                  npy_int32 precision, TrimMode trim_mode,
+                  npy_int32 precision, npy_int32 min_digits, TrimMode trim_mode,
                   npy_int32 digits_left, npy_int32 exp_digits)
 {
     npy_int32 printExponent;
@@ -1860,12 +1886,12 @@ FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
     char *pCurOut;
     npy_int32 numFractionDigits;
     npy_int32 leftchars;
+    npy_int32 add_digits;
 
     if (digit_mode != DigitMode_Unique) {
         DEBUG_ASSERT(precision >= 0);
     }
 
-
     DEBUG_ASSERT(bufferSize > 0);
 
     pCurOut = buffer;
@@ -1893,7 +1919,9 @@ FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
     }
 
     numDigits = Dragon4(mantissa, exponent, mantissaBit, hasUnequalMargins,
-                        digit_mode, CutoffMode_TotalLength, precision + 1,
+                        digit_mode, CutoffMode_TotalLength,
+                        precision < 0 ? -1 : precision + 1,
+                        min_digits < 0 ? -1 : min_digits + 1,
                         pCurOut, bufferSize, &printExponent);
 
     DEBUG_ASSERT(numDigits > 0);
@@ -1928,6 +1956,8 @@ FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
         --bufferSize;
     }
 
+    add_digits = digit_mode == DigitMode_Unique ? min_digits : precision;
+    add_digits = add_digits < 0 ? 0 : add_digits;
     if (trim_mode == TrimMode_LeaveOneZero) {
         /* if we didn't print any fractional digits, add the 0 */
         if (numFractionDigits == 0 && bufferSize > 1) {
@@ -1937,13 +1967,12 @@ FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
             ++numFractionDigits;
         }
     }
-    else if (trim_mode == TrimMode_None &&
-            digit_mode != DigitMode_Unique) {
-        /* add trailing zeros up to precision length */
-        if (precision > (npy_int32)numFractionDigits) {
+    else if (trim_mode == TrimMode_None) {
+        /* add trailing zeros up to add_digits length */
+        if (add_digits > (npy_int32)numFractionDigits) {
             char *pEnd;
             /* compute the number of trailing zeros needed */
-            npy_int32 numZeros = (precision - numFractionDigits);
+            npy_int32 numZeros = (add_digits - numFractionDigits);
 
             if (numZeros > (npy_int32)bufferSize - 1) {
                 numZeros = (npy_int32)bufferSize - 1;
@@ -1961,7 +1990,7 @@ FormatScientific (char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
      * when rounding, we may still end up with trailing zeros. Remove them
      * depending on trim settings.
      */
-    if (precision >= 0 && trim_mode != TrimMode_None && numFractionDigits > 0) {
+    if (trim_mode != TrimMode_None && numFractionDigits > 0) {
         --pCurOut;
         while (*pCurOut == '0') {
             --pCurOut;
@@ -2153,14 +2182,14 @@ Format_floatbits(char *buffer, npy_uint32 bufferSize, BigInt *mantissa,
         return FormatScientific(buffer, bufferSize, mantissa, exponent,
                                 signbit, mantissaBit, hasUnequalMargins,
                                 opt->digit_mode, opt->precision,
-                                opt->trim_mode, opt->digits_left,
-                                opt->exp_digits);
+                                opt->min_digits, opt->trim_mode,
+                                opt->digits_left, opt->exp_digits);
     }
     else {
         return FormatPositional(buffer, bufferSize, mantissa, exponent,
                                 signbit, mantissaBit, hasUnequalMargins,
                                 opt->digit_mode, opt->cutoff_mode,
-                                opt->precision, opt->trim_mode,
+                                opt->precision, opt->min_digits, opt->trim_mode,
                                 opt->digits_left, opt->digits_right);
     }
 }
@@ -3100,7 +3129,7 @@ Dragon4_Positional_##Type##_opt(npy_type *val, Dragon4_Options *opt)\
 \
 PyObject *\
 Dragon4_Positional_##Type(npy_type *val, DigitMode digit_mode,\
-                   CutoffMode cutoff_mode, int precision,\
+                   CutoffMode cutoff_mode, int precision, int min_digits, \
                    int sign, TrimMode trim, int pad_left, int pad_right)\
 {\
     Dragon4_Options opt;\
@@ -3109,6 +3138,7 @@ Dragon4_Positional_##Type(npy_type *val, DigitMode digit_mode,\
     opt.digit_mode = digit_mode;\
     opt.cutoff_mode = cutoff_mode;\
     opt.precision = precision;\
+    opt.min_digits = min_digits;\
     opt.sign = sign;\
     opt.trim_mode = trim;\
     opt.digits_left = pad_left;\
@@ -3136,7 +3166,8 @@ Dragon4_Scientific_##Type##_opt(npy_type *val, Dragon4_Options *opt)\
 }\
 PyObject *\
 Dragon4_Scientific_##Type(npy_type *val, DigitMode digit_mode, int precision,\
-                   int sign, TrimMode trim, int pad_left, int exp_digits)\
+                   int min_digits, int sign, TrimMode trim, int pad_left, \
+                   int exp_digits)\
 {\
     Dragon4_Options opt;\
 \
@@ -3144,6 +3175,7 @@ Dragon4_Scientific_##Type(npy_type *val, DigitMode digit_mode, int precision,\
     opt.digit_mode = digit_mode;\
     opt.cutoff_mode = CutoffMode_TotalLength;\
     opt.precision = precision;\
+    opt.min_digits = min_digits;\
     opt.sign = sign;\
     opt.trim_mode = trim;\
     opt.digits_left = pad_left;\
@@ -3166,8 +3198,8 @@ make_dragon4_typefuncs(LongDouble, npy_longdouble, NPY_LONGDOUBLE_BINFMT_NAME)
 
 PyObject *
 Dragon4_Positional(PyObject *obj, DigitMode digit_mode, CutoffMode cutoff_mode,
-                   int precision, int sign, TrimMode trim, int pad_left,
-                   int pad_right)
+                   int precision, int min_digits, int sign, TrimMode trim,
+                   int pad_left, int pad_right)
 {
     npy_double val;
     Dragon4_Options opt;
@@ -3176,6 +3208,7 @@ Dragon4_Positional(PyObject *obj, DigitMode digit_mode, CutoffMode cutoff_mode,
     opt.digit_mode = digit_mode;
     opt.cutoff_mode = cutoff_mode;
     opt.precision = precision;
+    opt.min_digits = min_digits;
     opt.sign = sign;
     opt.trim_mode = trim;
     opt.digits_left = pad_left;
@@ -3208,7 +3241,8 @@ Dragon4_Positional(PyObject *obj, DigitMode digit_mode, CutoffMode cutoff_mode,
 
 PyObject *
 Dragon4_Scientific(PyObject *obj, DigitMode digit_mode, int precision,
-                   int sign, TrimMode trim, int pad_left, int exp_digits)
+                   int min_digits, int sign, TrimMode trim, int pad_left,
+                   int exp_digits)
 {
     npy_double val;
     Dragon4_Options opt;
@@ -3217,6 +3251,7 @@ Dragon4_Scientific(PyObject *obj, DigitMode digit_mode, int precision,
     opt.digit_mode = digit_mode;
     opt.cutoff_mode = CutoffMode_TotalLength;
     opt.precision = precision;
+    opt.min_digits = min_digits;
     opt.sign = sign;
     opt.trim_mode = trim;
     opt.digits_left = pad_left;
diff --git a/numpy/core/src/multiarray/dragon4.h b/numpy/core/src/multiarray/dragon4.h
index 3a99bde6c53b..4b76bf9e582c 100644
--- a/numpy/core/src/multiarray/dragon4.h
+++ b/numpy/core/src/multiarray/dragon4.h
@@ -112,12 +112,12 @@ typedef enum TrimMode
     PyObject *\
     Dragon4_Positional_##Type(npy_type *val, DigitMode digit_mode,\
                               CutoffMode cutoff_mode, int precision,\
-                              int sign, TrimMode trim, int pad_left,\
-                              int pad_right);\
+                              int min_digits, int sign, TrimMode trim, \
+                              int pad_left, int pad_right);\
     PyObject *\
     Dragon4_Scientific_##Type(npy_type *val, DigitMode digit_mode,\
-                              int precision, int sign, TrimMode trim,\
-                              int pad_left, int exp_digits);
+                              int precision, int min_digits, int sign, \
+                              TrimMode trim, int pad_left, int exp_digits);
 
 make_dragon4_typedecl(Half, npy_half)
 make_dragon4_typedecl(Float, npy_float)
@@ -128,12 +128,13 @@ make_dragon4_typedecl(LongDouble, npy_longdouble)
 
 PyObject *
 Dragon4_Positional(PyObject *obj, DigitMode digit_mode, CutoffMode cutoff_mode,
-                   int precision, int sign, TrimMode trim, int pad_left,
-                   int pad_right);
+                   int precision, int min_digits, int sign, TrimMode trim,
+                   int pad_left, int pad_right);
 
 PyObject *
 Dragon4_Scientific(PyObject *obj, DigitMode digit_mode, int precision,
-                   int sign, TrimMode trim, int pad_left, int exp_digits);
+                   int min_digits, int sign, TrimMode trim, int pad_left,
+                   int exp_digits);
 
 #endif
 
diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c
index a0f7afeb5938..415aa7dacbbd 100644
--- a/numpy/core/src/multiarray/multiarraymodule.c
+++ b/numpy/core/src/multiarray/multiarraymodule.c
@@ -3481,7 +3481,7 @@ dragon4_scientific(PyObject *NPY_UNUSED(dummy),
         PyObject *const *args, Py_ssize_t len_args, PyObject *kwnames)
 {
     PyObject *obj;
-    int precision=-1, pad_left=-1, exp_digits=-1;
+    int precision=-1, pad_left=-1, exp_digits=-1, min_digits=-1;
     DigitMode digit_mode;
     TrimMode trim = TrimMode_None;
     int sign=0, unique=1;
@@ -3495,6 +3495,7 @@ dragon4_scientific(PyObject *NPY_UNUSED(dummy),
             "|trim", &trimmode_converter, &trim,
             "|pad_left", &PyArray_PythonPyIntFromInt, &pad_left,
             "|exp_digits", &PyArray_PythonPyIntFromInt, &exp_digits,
+            "|min_digits", &PyArray_PythonPyIntFromInt, &min_digits,
             NULL, NULL, NULL) < 0) {
         return NULL;
     }
@@ -3507,7 +3508,7 @@ dragon4_scientific(PyObject *NPY_UNUSED(dummy),
         return NULL;
     }
 
-    return Dragon4_Scientific(obj, digit_mode, precision, sign, trim,
+    return Dragon4_Scientific(obj, digit_mode, precision, min_digits, sign, trim,
                               pad_left, exp_digits);
 }
 
@@ -3522,7 +3523,7 @@ dragon4_positional(PyObject *NPY_UNUSED(dummy),
         PyObject *const *args, Py_ssize_t len_args, PyObject *kwnames)
 {
     PyObject *obj;
-    int precision=-1, pad_left=-1, pad_right=-1;
+    int precision=-1, pad_left=-1, pad_right=-1, min_digits=-1;
     CutoffMode cutoff_mode;
     DigitMode digit_mode;
     TrimMode trim = TrimMode_None;
@@ -3538,6 +3539,7 @@ dragon4_positional(PyObject *NPY_UNUSED(dummy),
             "|trim", &trimmode_converter, &trim,
             "|pad_left", &PyArray_PythonPyIntFromInt, &pad_left,
             "|pad_right", &PyArray_PythonPyIntFromInt, &pad_right,
+            "|min_digits", &PyArray_PythonPyIntFromInt, &min_digits,
             NULL, NULL, NULL) < 0) {
         return NULL;
     }
@@ -3552,8 +3554,8 @@ dragon4_positional(PyObject *NPY_UNUSED(dummy),
         return NULL;
     }
 
-    return Dragon4_Positional(obj, digit_mode, cutoff_mode, precision, sign,
-                              trim, pad_left, pad_right);
+    return Dragon4_Positional(obj, digit_mode, cutoff_mode, precision,
+                              min_digits, sign, trim, pad_left, pad_right);
 }
 
 static PyObject *
@@ -3572,7 +3574,7 @@ format_longfloat(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *kwds)
                 "not a longfloat");
         return NULL;
     }
-    return Dragon4_Scientific(obj, DigitMode_Unique, precision, 0,
+    return Dragon4_Scientific(obj, DigitMode_Unique, precision, -1, 0,
                               TrimMode_LeaveOneZero, -1, -1);
 }
 
diff --git a/numpy/core/src/multiarray/scalartypes.c.src b/numpy/core/src/multiarray/scalartypes.c.src
index 10f304fe7133..a001500b0a97 100644
--- a/numpy/core/src/multiarray/scalartypes.c.src
+++ b/numpy/core/src/multiarray/scalartypes.c.src
@@ -331,13 +331,13 @@ format_@name@(@type@ val, npy_bool scientific,
 {
     if (scientific) {
         return Dragon4_Scientific_@Name@(&val,
-                        DigitMode_Unique, precision,
+                        DigitMode_Unique, precision, -1,
                         sign, trim, pad_left, exp_digits);
     }
     else {
         return Dragon4_Positional_@Name@(&val,
                         DigitMode_Unique, CutoffMode_TotalLength, precision,
-                        sign, trim, pad_left, pad_right);
+                        -1, sign, trim, pad_left, pad_right);
     }
 }
 
diff --git a/numpy/core/tests/test_arrayprint.py b/numpy/core/tests/test_arrayprint.py
index 2c5f1577ddab..03eca2f7438e 100644
--- a/numpy/core/tests/test_arrayprint.py
+++ b/numpy/core/tests/test_arrayprint.py
@@ -759,6 +759,10 @@ def test_floatmode(self):
         assert_equal(repr(c),
             "array([1.00000000+1.00000000j, 1.12345679+1.12345679j])")
 
+        # test unique special case (gh-18609)
+        a = np.float64.fromhex('-1p-97')
+        assert_equal(np.float64(np.array2string(a, floatmode='unique')), a)
+
     def test_legacy_mode_scalars(self):
         # in legacy mode, str of floats get truncated, and complex scalars
         # use * for non-finite imaginary part
diff --git a/numpy/core/tests/test_scalarprint.py b/numpy/core/tests/test_scalarprint.py
index 6502ec4c1486..620472683391 100644
--- a/numpy/core/tests/test_scalarprint.py
+++ b/numpy/core/tests/test_scalarprint.py
@@ -9,7 +9,7 @@
 
 from tempfile import TemporaryFile
 import numpy as np
-from numpy.testing import assert_, assert_equal
+from numpy.testing import assert_, assert_equal, assert_raises
 
 class TestRealScalars:
     def test_str(self):
@@ -176,7 +176,8 @@ def test_dragon4(self):
                     "87538682506419718265533447265625")
 
         # largest numbers
-        assert_equal(fpos32(np.finfo(np.float32).max, **preckwd(0)),
+        f32x = np.finfo(np.float32).max
+        assert_equal(fpos32(f32x, **preckwd(0)),
                     "340282346638528859811704183484516925440.")
         assert_equal(fpos64(np.finfo(np.float64).max, **preckwd(0)),
                     "1797693134862315708145274237317043567980705675258449965989"
@@ -186,10 +187,66 @@ def test_dragon4(self):
                     "6580855933212334827479782620414472316873817718091929988125"
                     "0404026184124858368.")
         # Warning: In unique mode only the integer digits necessary for
-        # uniqueness are computed, the rest are 0. Should we change this?
-        assert_equal(fpos32(np.finfo(np.float32).max, precision=0),
+        # uniqueness are computed, the rest are 0.
+        assert_equal(fpos32(f32x),
                     "340282350000000000000000000000000000000.")
 
+        # Further tests of zero-padding vs rounding in different combinations
+        # of unique, fractional, precision, min_digits
+        # precision can only reduce digits, not add them.
+        # min_digits can only extend digits, not reduce them.
+        assert_equal(fpos32(f32x, unique=True, fractional=True, precision=0),
+                    "340282350000000000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=True, precision=4),
+                    "340282350000000000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=True, min_digits=0),
+                    "340282346638528859811704183484516925440.")
+        assert_equal(fpos32(f32x, unique=True, fractional=True, min_digits=4),
+                    "340282346638528859811704183484516925440.0000")
+        assert_equal(fpos32(f32x, unique=True, fractional=True,
+                                    min_digits=4, precision=4),
+                    "340282346638528859811704183484516925440.0000")
+        assert_raises(ValueError, fpos32, f32x, unique=True, fractional=False,
+                                          precision=0)
+        assert_equal(fpos32(f32x, unique=True, fractional=False, precision=4),
+                    "340300000000000000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=False, precision=20),
+                    "340282350000000000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=False, min_digits=4),
+                    "340282350000000000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=False,
+                                  min_digits=20),
+                    "340282346638528859810000000000000000000.")
+        assert_equal(fpos32(f32x, unique=True, fractional=False,
+                                  min_digits=15),
+                    "340282346638529000000000000000000000000.")
+        assert_equal(fpos32(f32x, unique=False, fractional=False, precision=4),
+                    "340300000000000000000000000000000000000.")
+        # test that unique rounding is preserved when precision is supplied
+        # but no extra digits need to be printed (gh-18609)
+        a = np.float64.fromhex('-1p-97')
+        assert_equal(fsci64(a, unique=True), '-6.310887241768095e-30')
+        assert_equal(fsci64(a, unique=False, precision=15),
+                     '-6.310887241768094e-30')
+        assert_equal(fsci64(a, unique=True, precision=15),
+                     '-6.310887241768095e-30')
+        assert_equal(fsci64(a, unique=True, min_digits=15),
+                     '-6.310887241768095e-30')
+        assert_equal(fsci64(a, unique=True, precision=15, min_digits=15),
+                     '-6.310887241768095e-30')
+        # adds/remove digits in unique mode with unbiased rnding
+        assert_equal(fsci64(a, unique=True, precision=14),
+                     '-6.31088724176809e-30')
+        assert_equal(fsci64(a, unique=True, min_digits=16),
+                     '-6.3108872417680944e-30')
+        assert_equal(fsci64(a, unique=True, precision=16),
+                     '-6.310887241768095e-30')
+        assert_equal(fsci64(a, unique=True, min_digits=14),
+                     '-6.310887241768095e-30')
+        # test min_digits in unique mode with different rounding cases
+        assert_equal(fsci64('1e120', min_digits=3), '1.000e+120')
+        assert_equal(fsci64('1e100', min_digits=3), '1.000e+100')
+
         # test trailing zeros
         assert_equal(fpos32('1.0', unique=False, precision=3), "1.000")
         assert_equal(fpos64('1.0', unique=False, precision=3), "1.000")
@@ -200,7 +257,9 @@ def test_dragon4(self):
         assert_equal(fsci32('1.5', unique=False, precision=3), "1.500e+00")
         assert_equal(fsci64('1.5', unique=False, precision=3), "1.500e+00")
         # gh-10713
-        assert_equal(fpos64('324', unique=False, precision=5, fractional=False), "324.00")
+        assert_equal(fpos64('324', unique=False, precision=5,
+                                   fractional=False), "324.00")
+
 
     def test_dragon4_interface(self):
         tps = [np.float16, np.float32, np.float64]


Code review comments:
diff_hunk: @@ -1048,6 +1052,10 @@ def format_float_scientific(x, precision=None, unique=True, trim='k',
     exp_digits : non-negative integer, optional
         Pad the exponent with zeros until it contains at least this many digits.
         If omitted, the exponent will be at least 2 digits.
+    min_digits : non-negative integer or None, optional
+        Minimum number of digits to print. Only has an effect if `unique=True`
+        in which case additional digits past those necessary to uniquely
+        identify the value may be printed, rounding the last additional digit.
path: numpy/core/arrayprint.py
text: @user1:
Hmm, an intended comment got lost. Was thinking we could try to simplify this a bit:
```suggestion
        Minimum number of digits to print. This only has an effect for `unique=True`.
        In that case more digits than necessary to uniquely identify the value may
        be printed. The last additional digit is rounded unbiased.
```
but happy to just leave it.

diff_hunk: @@ -1027,11 +1030,12 @@ def format_float_scientific(x, precision=None, unique=True, trim='k',
         If `True`, use a digit-generation strategy which gives the shortest
         representation which uniquely identifies the floating-point number from
         other values of the same type, by judicious rounding. If `precision`
-        was omitted, print all necessary digits, otherwise digit generation is
-        cut off after `precision` digits and the remaining value is rounded.
+        is given fewer digits than necessary can be printed, or if `min_digits`
path: numpy/core/arrayprint.py
text: @user1:
```suggestion
        is given fewer digits than necessary can be printed. If `min_digits`
```