c++ - Scaling a big integer by a double -

c++ - Scaling a big integer by a double -

i need scale big integers (several hundred bits) double. in particular, need compute

(m * factor) mod m

where m big integer , factor double. won't using libraries unless want phone call dozen lines of code in header file 'library'; hence big float math not alternative here.

knuth , gmp/mpir source code had no answers, , here found multiplication between big integers , doubles not applicable, since sec reply exotic , first loses much precision.

working first principles , simulating big integers uint64_t came (runnable 64-bit vc++ or 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 lower of limb_bits , dbl_mant_dig enum { limb_bits = sizeof(limb_t) * char_bit }; // simulate (m * factor) mod m 'big integer' m consisting of single limb void test_mod_mul (limb_t modulus, double factor) { assert( factor >= 0 ); // extract fractional part of factor , discard integer portion double ignored_integer_part; double fraction = std::modf(factor, &ignored_integer_part); // extract significand (aligned @ upper end of limb) , exponent int exponent; limb_t significand = limb_t(std::ldexp(std::frexp(fraction, &exponent), limb_bits)); // multiply modulus , single-limb significand; product have (n + 1) limbs limb_t hi; /* limb_t lo = */_umul128(modulus, significand, &hi); // result comprises @ n upper limbs of product; lowest limb // discarded in case, , potentially more. factors >= 1 handled well, // dropping modf() , 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 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)); }

the multiplication , shift done big integer math in app, principle should same.

is basic approach right or test working because i'm testing toy integers here? don't know first thing floating point math, , picked functions c++ reference.

clarification: multiplication onward done (partial) big integer math; here i'm using limb_t purpose in order tiny toy programme can posted , runs. final app using moral equivalent of gmp's mpn_mul_1() , mpn_rshift().

a floating point number nil product of 3 terms. 3 terms sign, significand (sometimes called mantissa), , exponent. value of these 3 terms computed as

(-1)sign * significand * baseexponent

the base of operations 2 although c++ standard doesn't guarantee that. correspondingly, computation becomes

(m * factor) mod m

== (m * (-1)sign * significand * baseexponent) mod m

== ((-1)sign(m) + sign * abs(m) * significand * baseexponent) mod m

computing sign of result should rather trivial. computing x * baseexponent rather straight forwards to: either suitable shift of bits if base of operations 2 or multiplication with/division powerfulness of base of operations (left shift or multiplication positive exponent, right shift or partition negative exponent). assuming big integer representation supports modulus operation already, interesting term multiplication of abs(m) * significand normal integer multiplication, although big integer representation. haven't checked closely think first reply linked (the 1 described "too exotic").

the remaining bit extract sign, significand, , exponent double. sign can determine comparing 0.0 , significand , exponent can obtained using frexp() (see this answer example). significand returned double, i.e., want multiply 2std::numeric_limits<double>::digits , adjust exponent appropriately (i haven't done in while, i.e., i'm not exclusively sure extact contract of frexp()).

to reply question: i'm not familiar gmp operations using think operations indeed execute computation described above.

c++ floating-point biginteger arbitrary-precision