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

Question

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

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

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

У Кнута и исходного кода 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.

Answer 1

Число с плавающей запятой - не что иное, как произведение трех терминов. Три термина - это знак, значимое (иногда называемое мантиссой) и показатель степени. Значение этих трех терминов вычисляется как

(-1) ^знак * значимое и * базовый ^{показатель}

База обычно равна 2, хотя стандарт C++ не гарантирует этого. Соответственно, ваш расчет становится

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

== (M * (-1) ^знак * значимое и * базовый ^{показатель}) мод М

== ((-1) ^{знак (M) + знак} * abs (M) * значение и * базовый ^{показатель}) mod M

Вычисление знака результата должно быть довольно тривиальным. Вычислить базовый ^{показатель} X * довольно просто: это либо подходящий сдвиг битов, если основание равно 2, либо умножение на / деление на степень основания (сдвиг влево или умножение для положительного показателя, сдвиг вправо или деление для отрицательного экспонент). Предполагая, что ваше большое целочисленное представление уже поддерживает операцию модуля, единственным интересным термином является умножение abs (M) * значимый, но это обычное целочисленное умножение, хотя и для вашего большого целочисленного представления. Я не проверял слишком близко, но я думаю, что это то, что делает первый ответ, на который вы ссылаетесь (тот, который вы назвали "слишком экзотическим").

Оставшийся бит состоит в том, чтобы извлечь знак, значение и показатель из double, Знак можно легко определить путем сравнения с 0.0 и значение, и показатель степени могут быть получены с использованием frexp() (см. этот ответ, например). Значение и возвращается как double вы, вероятно, хотите умножить его на 2 std::numeric_limits<double>::digits и соответствующим образом скорректировать показатель (я не делал это некоторое время, т. е. я не совсем уверен в предварительном договоре frexp()).

Чтобы ответить на ваш вопрос: я не знаком с операциями GMP, которые вы используете, но я думаю, что операции, которые вы действительно выполняете, выполняют вычисления, описанные выше.