Инверсия fmod

fmod(x,10) = 0

Это учитывается для каждого целого числа, заканчивающегося нулем, например 0, -10, 20, -30,..., 100, -110, 120, -130,..., 200, -210, 220, -230. , ..., 300, -310, 320, -330, ..., 1000, -1010, 1020, -1030, ..., 1100, -1110, 1120, -1130, ... (вы также встречаются ±123456789876543210 и ±98765432123456890 и многие-многие другие) :-)

Вы понимаете, почему на ваш вопрос не может быть ответа? :-)

Из этой ссылки на std::fmod:

Остаток операции деления x / y с плавающей запятой, вычисленный этой функцией, точно равен значению x - rem * y, где rem — это x / y с усеченной дробной частью.

_{[Emphasis mine]}

Поскольку дробная часть усекается, невозможно вернуть исходное значение.

И это до вовлечения этого.

Далее я предполагаю вычисление double, где double отображается в формате IEEE-754 binary64.

Любое общее деление дает две порции информации: частное и остаток. fmod (a, b) возвращает остаток a / b, когда неявное частное преобразуется в целое число с округлением до нуля (усечение). Исходное делимое может быть вычислено из делителя b, усеченного частного и остатка, возвращаемого fmod(), при условии, что частное точно представимо в виде double. Без частного недостаточно информации для повторного вычисления дивиденда.

Поскольку операнд double обеспечивает 52 сохраненных значащих (мантиссы) бита, частное в общем случае точно представимо только в том случае, если разница двоичных показателей делимого и делителя не превышает 52. Очевидно, нам нужно избегать делимых и делителей, которые равны бесконечности. или NaN, и делитель не может быть равен нулю. Пересчет дивиденда как (частное * делитель + остаток) должен выполняться с однократным округлением, поэтому необходимо использовать плавное умножение-сложение (FMA).

Следующая программа демонстрирует успешный пересчет дивиденда при указанных ограничениях. Обратите внимание, что для достижения предполагаемой функциональности fdiv_rz() требуется настроить компилятор для самого строгого соблюдения IEEE-754. В противном случае компилятор может рассматривать операцию деления как независимую от вызовов fesetround(), в результате чего порядок операций будет изменен таким образом, что деление больше не будет зажато между вызовами fesetround() по мере необходимости. Из того, что я вижу, это особенно влияет на gcc, тогда как clang по умолчанию ведет себя лучше. Я использовал компилятор Intel для создания этого кода, который позволяет легко управлять поведением с плавающей запятой. В частности, я скомпилировал с icl /Qstd=c++11 /Ox /QxHOST /W4 /fp:strict fdiv_fmod.cpp.

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <fenv.h>

#define N (2000000000)  // number of random test cases

/* divide a by b, rounding result towards zero */
double fdiv_rz (double a, double b) 
{
    double quot;
    int curr_mode = fegetround ();
    fesetround (FE_TOWARDZERO);
    quot = a / b;
    fesetround (curr_mode);
    return quot;
}

/* reinterpret 64-bit unsigned integer as an IEEE-754 binary64 */
double uint64_as_double (uint64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

int main (void)
{
    double a, b, t, quot, rem;
    int expa, expb;

    for (int i = 0; i < N; i++) {
        do {
            a = uint64_as_double (KISS64);
        } while (isinf (a) || isnan (a));
        (void)frexp (a, &expa);
        do {
            b = uint64_as_double (KISS64);
            (void)frexp (b, &expb);
        } while ((b == 0) || isinf (b) || isnan (b) || (abs (expa - expb) > 52));
        
        quot = trunc (fdiv_rz (a, b));
        rem  = fmod (a, b);
        t = (fabs (quot) < 1.0) ? rem : fma (quot, b, rem);

        if (t != a) {
            printf ("a=% 23.13a  b=% 23.13a  quot=% 23.13a  rem=%23.13a  t=% 23.13a\n",
                    a, b, quot, rem, t);
            return EXIT_FAILURE;
        }
    }
    printf ("test passed\n");
    return EXIT_SUCCESS;
}