Масштабирование большого целого числа двойным

Мне нужно масштабировать большие целые числа (несколько сотен бит) в два раза. В частности, мне нужно вычислить

(М * фактор) мод М

где M большое целое число и фактор это двойной. Я не буду использовать какие-либо библиотеки, если вы не хотите называть дюжину строк кода в заголовочном файле «библиотекой»; следовательно, большая математика с плавающей точкой здесь не вариант.

Кнут и GMPИсходный код / ​​MPIR не имеет ответов, и здесь я нашел только Умножение между большими целыми и двойными что не совсем применимо, так как второй ответ слишком экзотичен, а первый теряет слишком много точности.

Работая с первых принципов и симулируя большие целые числа с помощью uint64_t, я придумал это (работающий с 64-битным VC ++ или gcc / MinGW64):

#include <cassert>
#include <cfloat>
#include <climits>
#include <cmath>
#include <cstdint>
#include <cstdio>

#include <intrin.h>   // VC++, MinGW

#define PX(format,expression)  std::printf("\n%35s  ==  " format, #expression, expression);

typedef std::uint64_t limb_t;

// precision will be the lower of LIMB_BITS and DBL_MANT_DIG
enum {  LIMB_BITS = sizeof(limb_t) * CHAR_BIT  };

// simulate (M * factor) mod M with a 'big integer' M consisting of a single limb
void test_mod_mul (limb_t modulus, double factor)
{
assert( factor >= 0 );

// extract the fractional part of the factor and discard the integer portion

double ignored_integer_part;
double fraction = std::modf(factor, &ignored_integer_part);

// extract the significand (aligned at the upper end of the limb) and the exponent

int exponent;
limb_t significand = limb_t(std::ldexp(std::frexp(fraction, &exponent), LIMB_BITS));

// multiply modulus and single-limb significand; the product will have (n + 1) limbs

limb_t hi;
/* limb_t lo = */_umul128(modulus, significand, &hi);

// The result comprises at most n upper limbs of the product; the lowest limb will be
// discarded in any case, and potentially more. Factors >= 1 could be handled as well,
// by dropping the modf() and handling exponents > 0 via left shift.

limb_t result = hi;

if (exponent)
{
assert( exponent < 0 );

result >>= -exponent;
}

PX("%014llX", result);
PX("%014llX", limb_t(double(modulus) * fraction));
}

int main ()
{
limb_t const M = 0x123456789ABCDEull;  // <= 53 bits (for checking with doubles)

test_mod_mul(M, 1 - DBL_EPSILON);
test_mod_mul(M, 0.005615234375);
test_mod_mul(M, 9.005615234375);
test_mod_mul(M, std::ldexp(1, -16));
test_mod_mul(M, std::ldexp(1, -32));
test_mod_mul(M, std::ldexp(1, -52));
}

Умножение и сдвиг будут выполняться с помощью целочисленной математики в моем приложении, но принцип должен быть таким же.

Правильный ли базовый подход или тест работает только потому, что я здесь тестирую с целыми числами игрушек? Я не знаю ничего о математике с плавающей запятой, и я выбрал функции из C ++ ссылка.

Пояснение: все, начиная с умножения и далее, будет сделано с (частичной) большой целой математикой; здесь я использую только limb_t для этой цели, чтобы получить крошечную игрушечную программу, которая может быть размещена и действительно запускается. Окончательное приложение будет использовать моральный эквивалент mpn_mul_1 () и mpn_rshift () GMP.