2

I am managing some big (128~256bits) integers with gmp. It has come a point were I would like to multiply them for a double close to 1 (0.1 < double < 10), the result being still an approximated integer. A good example of the operation I need to do is the following:

int i = 1000000000000000000 * 1.23456789

I searched in the gmp documentation but I didn't find a function for this, so I ended up writing this code which seems to work well:

mpz_mult_d(mpz_class & r, const mpz_class & i, double d, int prec=10) {
  if (prec > 15) prec=15; //avoids overflows
  uint_fast64_t m = (uint_fast64_t) floor(d);
  r = i * m;
  uint_fast64_t pos=1;
  for (uint_fast8_t j=0; j<prec; j++) {
    const double posd = (double) pos;
    m = ((uint_fast64_t) floor(d * posd * 10.)) -
        ((uint_fast64_t) floor(d * posd)) * 10;
    pos*=10;
    r += (i * m) /pos;
  }
}

Can you please tell me what do you think? Do you have any suggestion to make it more robust or faster?

DarioP
  • 5,377
  • 1
  • 33
  • 52
  • 1
    That's a question for [Code Review](http://codereview.stackexchange.com/), not for StackOverflow :) – Morwenn May 27 '13 at 12:40
  • Oh sorry I didn't know about that branch. However that's just my solution to a very precise problem. Please, consider to answer the general question and only eventually comment on the code. Thanks! – DarioP May 27 '13 at 12:45
  • if you need an approximation, why don't you convert the big integer to double? – Karoly Horvath May 27 '13 at 12:47
  • @KarolyHorvath: because doing that I am compressing the 256bits of the integer into the 52bits of the double mantissa. That totally waste the use of ints! A possibility would be to cast the big int to a big float, but I think it's problematic in terms of size, speed and accuracy. – DarioP May 27 '13 at 12:54
  • so what? The double you are multiplying with doesn't have bigger precision either.. if you need precision, use big float, or do the multiplication with 123456789 (store it as decimal) – Karoly Horvath May 27 '13 at 12:56
  • @KarolyHorvath that's more or less what I'm doing in the code but trying to keep a reasonable memory usage. Double shouldn't have problems to handle values like 1.23456789 if I just look at the first 10 or something figures. PS. I need a double because that comes from some nasty physical quantities. – DarioP May 27 '13 at 13:05
  • 1
    GMP supports arbitrary precision rational and floating point numbers. Convert your two values to one of these formats and let GMP do the multiply. – brian beuning May 27 '13 at 13:33
  • If it comes from "physical quantities" (do you mean measurements? something else?) isn't it an approximation already? I mean, even if the value 1.23456789 can be represented exactly, it isn't the exact quantity, but merely an approximation with as much accuracy and precision as provided by the measuring device. – R. Martinho Fernandes May 27 '13 at 13:52
  • 1
    You should look to MPFR (http://www.mpfr.org/), which is the equivalent of GMP for floating-point numbers. It is as easy/simple to use as GMP. – Vincent Jul 12 '13 at 17:24

3 Answers3

0

this is what you wanted:

// BYTE lint[_N]   ... lint[0]=MSB, lint[_N-1]=LSB
void mul(BYTE *c,BYTE *a,double b)  // c[_N]=a[_N]*b
    {
    int i; DWORD cc;
    double q[_N+1],aa,bb;
    for (q[0]=0.0,i=0;i<_N;)        // mul,carry down
        {
        bb=double(a[i])*b; aa=floor(bb); bb-=aa;
        q[i]+=aa; i++;
        q[i]=bb*256.0;
        }
    cc=0; if (q[_N]>127.0) cc=1.0;  // round
    for (i=_N-1;i>=0;i--)           // carry up
        {
        double aa,bb;
        cc+=q[i];
        c[i]=cc&255;
        cc>>=8;
        }
    }

_N is number of bits/8 per large int, large int is array of _N BYTEs where first byte is MSB (most significant BYTE) and last BYTE is LSB (least significant BYTE) function is not handling signum, but it is only one if and some xor/inc to add.

trouble is that double has low precision even for your number 1.23456789 !!! due to precision loss the result is not exact what it should be (1234387129122386944 instead of 1234567890000000000) I think my code is mutch quicker and even more precise than yours because i do not need to mul/mod/div numbers by 10, instead i use bit shifting where is possible and not by 10-digit but by 256-digit (8bit). if you need more precision than use long arithmetic. you can speed up this code by using larger digits (16,32, ... bit)

My long arithmetics for precise astro computations are usually fixed point 256.256 bits numbers consist of 2*8 DWORDs + signum, but of course is much slower and some goniometric functions are realy tricky to implement, but if you want just basic functions than code your own lon arithmetics is not that hard.

also if you want to have numbers often in readable form is good to compromise between speed/size and consider not to use binary coded numbers but BCD coded numbers

Spektre
  • 49,595
  • 11
  • 110
  • 380
0

I am not so familiar with either C++ or GMP what I could suggest source code without syntax errors, but what you are doing is more complicated than it should and can introduce unnecessary approximation.

Instead, I suggest you write function mpz_mult_d() like this:

mpz_mult_d(mpz_class & r, const mpz_class & i, double d) {
  d = ldexp(d, 52); /* exact, no overflow because 1 <= d <= 10 */
  unsigned long long l = d; /* exact because d is an integer */
  p = l * i; /* exact, in GMP */
  (quotient, remainder) = p / 2^52; /* in GMP */

And now the next step depends on the kind of rounding you wish. If you wish the multiplication of d by i to give a result rounded toward -inf, just return quotient as result of the function. If you wish a result rounded to the nearest integer, you must look at remainder:

 assert(0 <= remainder);  /* proper Euclidean division */
 assert(remainder < 2^52);
 if (remainder < 2^51) return quotient;
 if (remainder > 2^51) return quotient + 1; /* in GMP */
 if (remainder == 2^51) return quotient + (quotient & 1); /* in GMP, round to “even” */

PS: I found your question by random browsing but if you had tagged it “floating-point”, people more competent than me could have answered it quickly.

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
0

Try this strategy:

  1. Convert integer value to big float
  2. Convert double value to big float
  3. Make product

  4. Convert result to integer

    mpf_set_z(...)
mpf_set_d(...)
mpf_mul(...)
mpz_set_f(...)
GeScript
  • 1