Trying to multiply a number from 0..4096 by 0.3 using only integer math

Question

I'm trying to wrap my head around how I'd multiply a number from 0..4096 by 0.3 using just integers with shift operations and scaling, without dividing, in C. I'm new to this stuff and any input or step by step suggestions would be extremely helpful.

Answer 1

Multiplying by 0.3 is the same as multiplying by (0.3*2^n), then dividing by 2^n. The second stage is equivalent to a shift right of n.

But what's the best value of n?

To find this, take your largest integer and find the largest value of n such that you can multiply by (0.3*2^n) without an overflow. With 64 bit unsigned integers and 4096 as the maximum value, you need

0.3*2^n <= 2^(64-12)

or

0.3 <= 2^(64-12-n)

That inequality has maximum n when the RHS is equal to 0.5, so

2^-1 = 2^(64-12-n)

so -1 = 64-12-n, n = 64-12+1 = 53.

So the answer is multiply by 2^53*0.3 then shift right by 53, i.e.

/* designed to work with input values 0 .. 4096 only */
uint64_t
multiplyby0_3 (uint64_t x)
{
    return (x * 2702159776422297ULL) >> 53;
}

To check that doesn't overflow and we've got the best n, from bc:

2702159776422297*4096 = 11068046444225728512
2^64                  = 18446744073709551616

IE it won't overflow, but if we multiplied it by 2 again, it would.

For 32 bit integers, the answer is multiply by 2^21*0.3 then shift right by 21, i.e.

/* designed to work with input values 0 .. 4096 only */
uint32_t
multiplyby0_3 (uint32_t x)
{
    return (x * 629146U) >> 21;
}

And finally, you can decompose any multiply into a number of additions by looking at the binary 1's in the multiplier. So you permit 'scaling' and I assumed that meant multiply. If not, here is the 32 bit version (64 bit version left as an exercise for the reader) exploiting the fact that 629146 is 10011001100110011010 (a neat pattern due to the recurring binary fraction). We'll round the other way and use 10011001100110011001 instead.

/* designed to work with input values 0 .. 4096 only */
uint32_t
multiplyby0_3 (uint32_t x)
{
    uint32_t y;
    x += x<<3;    /* * 1001 */
    y = x<<4 + x; /* * 10011001 */
    y += y<<8;    /* * 1001100110011001 */
    y += x<<16;   /* * 10011001100110011001 */
    return y >> 21;
}

Answer 2

If you have fast integer multiplication, but you don't have integer division, you can get a reasonable approximation by multiplying by 1229, then shifting right by 12 bits. For instance:

>> 100 * 1229
122900
>> 122900 >> 12
30

This works because 1229 is about 0.3 * 1 << 12. The "real" value is 1228.8, so the estimate will be high by 1 in some cases (68 out of 4097 values). It will never be off by more than 1, though.

Answer 3

Using a classic hack for divu10() below. I prefer the *n and shift approach, but thought I offer another POV.

unsigned mult3tenths(unsigned x) {
  return divu10(x*3);
}

unsigned divu10(unsigned n) {
  unsigned q, rem;
  q = (n >> 1) + (n >> 2);
  q = q + (q >> 4);
  q = q + (q >> 8);
  q = q + (q >> 16);
  q = q >> 3;
  rem = n - q*10;
  return q + ((rem + 6) >> 4);  
}

Answer 4

I know that this is not a coding example, sorry about that, but if the set is so small (range [0;4096]), why not create a block of results and use a pointer to extract values? It will greatly reduce the GPU cycles, since there are not memory restrictions.

Answer 5

Simply tried a lot of combinations of the form (A*x + B) >> n with the below test code and came up with:

// Scale by 4915, then shift 14.
int Mult3Tenths(int x) {
  return (x*4915 + 0) >> 14;  // Use 4915L if `int` is 16 bit.
}

Test code

#define N3 (4096)
int main(void) {
  int target[N3 + 1];
  unsigned i;
  for (i = 0; i <= N3; i++) {
    target[i] = 0.3 * i;
  }

  // form (A*x + B) >> n
  int A, B, n;
  int besti = 0;
  for (n = 0; n < 31; n++) {
    int Amin = ((N3 * 0.3 - 1) * (1 << n) - (1 << n)) / N3 - 1;
    int Amax = ((N3 * 0.3 + 1) * (1 << n) + (1 << n)) / N3 + 1;
    for (A = Amin; A <= Amax; A++) {
      int Bmax = 1 << n;
      for (B = -Bmax; B <= Bmax; B++) {
        for (i = 0; i <= N3; i++) {
          int y = (A * i + B) >> n;
          if (y != target[i])
            break;
          if (i > besti) {
            besti = i;
            if (i == N3) {
              printf("i:%i A:%d B:%d n:%d\n", i, A, B, n);
              printf("!!!\n");
              exit(0);
            }
          }
        }
      }
    }
  }
  printf("???\n");
  return 0;
}