How do I compute the first decimal digit of a division between large numbers?

Question

I have two unsigned long longs X and Y, where X < Y, but both may be very large. I want to compute the first digit past the decimal point of X / Y. For example, if X is 11 and Y is 14, then 11 / 14 is .785…, so the result should be 7.

(X * 10) / Y would work, except it produces the wrong result if X * 10 overflows. Conversion to double would work if I had reason to believe it is precise enough to compute the right result.

This is in C. Thanks for any help!

Answer 1

Conversion to double would work if I had reason to believe it is precise enough to compute the right result.

If you only need the first digit, then surely double is precise enough.

Edit: wrang-wrang's counter-example in the comments proves me wrong.

Answer 2

I'm unwilling to guarantee an accurate result via floating point cast unless the mantissa has enough bits to exactly represent all the integral values. For instance, tonsider y=9223372036854775807 and x = (y div 10) - 1 = 922337203685477579. where "div" is integer division. x/y is 0.09999999999999999981568563067746..., but using doubles gets you >= 0.1. This is a result of doubles having only 52 digits of precision in the significand (while y needs 61 bits and x about 58)

You may be able to use 80bit or 128bit FP precision, in which case you'd get the right answer because the mantissa would be >=64bit (ULL are 64bit, right?) and so you'd losslessly represent the numbers.

I'd start with an approximation (either using integer or FP arithmetic) then trial multiplication to see if the answer should be 1 less or more. The key insight is that you can still compare two possibly overflowed ints as long as you know that the difference between the two quantities is less than half the max unsigned int. This sort of comparison technique is necessary when e.g. TCP sequence numbers overflow.

In case you want to use only integer arithmetic, the below function "fdd(x,y)" works. I've included a main() to show some results:

#include <iostream>
using namespace std;

typedef unsigned char ull; // change char to any integral type e.g. long long

const ull maxull=(ull)-1;
const ull halfull = maxull/2;
typedef unsigned long long asint;

// x = X mod (maxull+1), y= Y mod (maxull+1).  we only know x and y
// if we assume |X-Y|<halfull, then we return X<Y:
inline bool less_mod_near(ull x, ull y) {

    return (x<=halfull == y<=halfull) ? x<y : y>x;
}

// assuming x<y, return first decimal digit of 10x/y (return is in [0..9])
inline int fdd(ull x, ull y) { 
// assert(x<y);
 if (x<=maxull/10) return (10*x)/y; 
  // for speed, and to ensure that y>10 to avoid division by 0 later
 ull r=y/10;
 if (r*10==y) return x/r;
 ull ub=x/(r+1); // ub >= 10x div y (without overflow)
 ull x10=x*10; // allow overflow
 cout<<"ub="<<(asint)ub<<" x10="<<(asint)x10<<" r="<<(asint)r<<" ";
 return less_mod_near(x10,ub) ? ub-1 : ub; 
  // we already handled the 10 evenly divides y case
}

int pdd(ull x, ull y,ull mustbe)
{
    ull d=fdd(x,y);
    cout << (asint)x << '/' << (asint)y << " = ." << (asint)d << "...";
    if (d!=mustbe) cout << " (should be "<<(asint)mustbe<<")";
    cout<<endl;
//    assert(a==d);
}

int main() {
    pdd(0,1,0);
    pdd(1,2,5);
    pdd(11,101,1);
    pdd(10,101,0);
    pdd(49,69,7);
    pdd(50,69,7);
    pdd(48,69,6);
    pdd(160,200,8);
    pdd(161,200,8);
    pdd(159,200,7);
    pdd(254,255,9);
}

Output:

0/1 = .0...
1/2 = .5...
11/101 = .1...
10/101 = .0...
ub=7 x10=234 r=6 49/69 = .7...
ub=7 x10=244 r=6 50/69 = .7...
ub=6 x10=224 r=6 48/69 = .6...
160/200 = .8...
161/200 = .8...
159/200 = .7...
ub=9 x10=236 r=25 254/255 = .9...

Answer 3

Converting to Double is going to give up 12 bits of precision on both the dividend and the donominator out of 64. The majority of the time it will give you the correct answer, however it will occationally have rounding errors unless it's stored usng 80bit floating format.

Edit: Unless I'm just too sleepy to see the error, I think wrang-wrang's answer will work. I got work in morning, 6 hour drive to customer's site. Ugg. night.

Edit: One last thing. x86 uses 80bit internal representation. I think there are opcodes that will garentee int64 to float80 conversion, if you want to toss a couple asm instructions around. That would be a more elegant and certainly faster solution than a pure C implementation, although it wouldn't be portable.

Answer 4

If you are working on PPC or x86-64, you have access to native 64-bit integer division operations. PPC does it always, x86 supports it with the divq opcode.

Your compiler should be clever enough to emit those instructions if you simply do

unsigned long long X, Y; // assuming unsigned long long == uint64 
unsigned long long result = X/Y; // 64 bit divison.

but you should verify this in the emitted machine code (/FAcs in Visual Studio, -S in GCC) to be sure; the gap between "the compiler should" and "the compiler does" is huge. With GCC you may need to force it by using the __uint128_t intrinsic type; with Visual Studio, if your version doesn't emit divq, you may have to resort to inline assembly just for that one op.

Once you've got your result, sounds like you've got a handle on finding the top base-10 sigfig already. 64-bit mod will work in just the same way as 64-bit division, because it's the same hardware op.

See also prior questions.

Answer 5

X % Y gives the remainder R

You can then compute the fractional part of your answer from R and Y

R / Y

To get the 1st digit use integer division:

(long)((long)10*R/Y)

This should round down the numbers and drop any extra decimals.

Edit:

To match your question (for those who want to be picky) its

Y % X = R

and

R / X

Answer 6

If you are prpared to accept the range limitation you could do it all in integer arithmatic:-

(X * 10) / Y

In your example:

(11 * 10) / 14

=> 110 / 14

=> 7

The limitation being you have reduced the maximum value of X by a factor of 10.

Answer 7

How about

x / ((y + 9) / 10)

The y + 9 is for rounding up the quotient y/10 in the denominator, so the overall result is rounded down.

For large x and y it should be almost always correct, but it's not perfect, it produces 5 for your example 11 / 14.

The problem is that you loose information due to the division. Since multiplication by ten already overflows, you can't get around that except by using a bigger number data type.