Question about mpf_t comparison behavior for x - tiny

[Marc Glisse]

Hello,

I am not a specialist of MPF, but I think, as you suspect, that this is a consequence of how MPF works, and in particular the way it does not try to be very strict about the number of digits it keeps (the exponents don't need to be extreme to reproduce this phenomenon). In your example, if I set the precision to 64, I see

(gdb) p wu
$1 = {{_mp_prec = 2, _mp_size = -3, _mp_exp = 18014398509481984, _mp_d = 0x555555559330}}
(gdb) p d_minus_pivmin
$2 = {{_mp_prec = 2, _mp_size = -2, _mp_exp = 18014398509481984, _mp_d = 0x555555559370}}
(gdb) p/x *(mp_limb_t[3]*)wu._mp_d
$3 = {0xb1d0d65a96f5c9b, 0x1cd345dcc8169fef, 0x6ba27e656b4eb57a}
(gdb) p/x *(mp_limb_t[2]*)d_minus_pivmin._mp_d
$4 = {0x1cd345dcc8169fef, 0x6ba27e656b4eb57a}

That is, the initialization from a string gives more digits than we ask for, while later operations give a bit less, so an operation that should have been the identity may "randomly" end up truncating instead. Maybe mpf_set_str should be fixed to use fewer limbs, but the ultimate fix is to use MPFR instead.

----- Mail original -----
> De: "Maho NAKATA" <maho.nakata at gmail.com>
> À: gmp-discuss at gmplib.org
> Envoyé: Vendredi 24 Avril 2026 13:02:16
> Objet: Question about mpf_t comparison behavior for x - tiny

> Hello,
> 
> I hope this is an appropriate place to ask this question.
> 
> While debugging a numerical issue in MPLAPACK, I encountered an `mpf_t`
> behavior that I found a little surprising, and I would be very grateful
> if someone could tell me whether this is expected GMP behavior.
> 
> The pattern is roughly the following:
> 
> - `wu == d`
> - `pivmin > 0`
> - but `wu < d - pivmin`
> 
> So mathematically,
> 
> `wu - (d - pivmin) = (wu - d) + pivmin`
> 
> seems as though it should remain positive if `wu == d` and `pivmin > 0`.
> 
> If this were behaving in the most straightforward way, I would expect
> either
> 
> - `d - pivmin` to compare equal to `d`, in which case both comparisons
>  would agree, or
> - `d - pivmin` to be strictly smaller than `d`, in which case `wu == d`
>  should still imply `wu > d - pivmin`.
> 
> In either case, I would not expect to see
> 
> - `mpf_cmp(wu, d) == 0`
> - but at the same time `mpf_cmp(wu, d - pivmin) < 0`.
> 
> However, with `mpf_t` I observe exactly that, which is why I thought it
> might be worth asking whether this is simply a consequence of `mpf_t`
> semantics at extreme exponents, or whether I may be misunderstanding
> something about how these values are represented and compared.
> 
> I reduced the issue to the following small example:
> 
> ```cpp
> #include <cstdio>
> #include <cstdlib>
> #include <gmp.h>
> 
> static void set_or_die(mpf_t x, const char *s) {
>    if (mpf_set_str(x, s, 10) != 0) {
>        std::fprintf(stderr, "mpf_set_str failed for: %s\n", s);
>        std::exit(1);
>    }
> }
> 
> int main() {
>    const mp_bitcnt_t prec = 512;
> 
>    mpf_t d, wu, pivmin, d_minus_pivmin, diff_d, diff_p;
>    mpf_inits(d, wu, pivmin, d_minus_pivmin, diff_d, diff_p, nullptr);
> 
>    mpf_set_prec(d, prec);
>    mpf_set_prec(wu, prec);
>    mpf_set_prec(pivmin, prec);
>    mpf_set_prec(d_minus_pivmin, prec);
>    mpf_set_prec(diff_d, prec);
>    mpf_set_prec(diff_p, prec);
> 
>    mpf_set_str(d,
> "-2.4616939245935578262292233824060520172893599354879699820500368499e+347063955532709820",
> 10);
>    mpf_set_str(wu,
> "-2.4616939245935578262292233824060520172893599354879699820500368499e+347063955532709820",
> 10);
>    mpf_set_str(pivmin,
> "1.7019382623481672278259575819240965611355119939659392498165297957e-1388255822130839283",
> 10);
> 
>    mpf_sub(d_minus_pivmin, d, pivmin);
>    mpf_sub(diff_d, wu, d);
>    mpf_sub(diff_p, wu, d_minus_pivmin);
> 
>    gmp_printf("d                    = %.40Fe\n", d);
>    gmp_printf("wu                   = %.40Fe\n", wu);
>    gmp_printf("pivmin               = %.40Fe\n", pivmin);
>    gmp_printf("d_minus_pivmin       = %.40Fe\n", d_minus_pivmin);
>    gmp_printf("wu_minus_d           = %.40Fe\n", diff_d);
>    gmp_printf("wu_minus_d_minus_piv = %.40Fe\n", diff_p);
> 
>    std::printf("mpf_cmp(wu, d) = %d\n", mpf_cmp(wu, d));
>    std::printf("mpf_cmp(wu, d_minus_pivmin) = %d\n", mpf_cmp(wu,
> d_minus_pivmin));
> 
>    mpf_clears(d, wu, pivmin, d_minus_pivmin, diff_d, diff_p, nullptr);
>    return 0;
> }
> ```
> 
> On my system, this produces:
> 
> ```text
> d = -2.4616939245935578262292233824060520172894e+347063955532709820
> wu = -2.4616939245935578262292233824060520172894e+347063955532709820
> pivmin = 1.7019382623481672278259575819240965611355e-1388255822130839283
> d_minus_pivmin =
> -2.4616939245935578262292233824060520172894e+347063955532709820
> wu_minus_d = 0.0000000000000000000000000000000000000000e+00
> wu_minus_d_minus_pivmin =
> -1.6856831808738259700179295042359981612965e+347063955532709647
> mpf_cmp(wu, d) = 0
> mpf_cmp(wu, d_minus_pivmin) = -1
> ```
> 
> My environment is:
> 
> - GMP 6.3.0
> - GCC 13.3.0
> - `gcc (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0`
> - Ubuntu 24.04 amd64
> 
> I understand that `mpf_t` is not intended to provide MPFR-style
> guarantees, so this may well be normal behavior. I am mainly hoping to
> understand whether this is an expected consequence of `mpf_t` arithmetic
> at extreme exponents, or whether it suggests something more subtle.
> 
> If it would be helpful, I would be happy to provide:
> 
> - a version that prints internal fields / limbs
> - any additional environment details
> 
> Thank you very much for your time.
> 
> Best regards,
> Nakata Mho
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-discuss