Пытается умножить число от 0..4096 на 0,3, используя только целую математику

Ответ 1

Умножение на 0,3 такое же, как умножение на (0.3*2^n), затем деление на 2^n. Второй этап эквивалентен праву сдвига n.

Но какое лучшее значение n?

Чтобы найти это, возьмите наибольшее целое число и найдите наибольшее значение n, которое можно умножить на (0.3*2^n) без переполнения. Для 64-битных целых чисел без знака и 4096 в качестве максимального значения вам нужно

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

или

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

Это неравенство имеет максимум n, когда RHS равно 0,5, поэтому

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

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

Итак, ответ умножается на 2^53*0.3, а затем сдвигается вправо на 53, т.е.

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

Чтобы проверить, что не переполняется, и у нас есть лучший n, от bc:

2702159776422297*4096 = 11068046444225728512
2^64                  = 18446744073709551616

IE он не будет переполняться, но если мы снова умножим его на 2, это будет.

Для 32-битных целых чисел ответ умножается на 2^21*0.3, а затем сдвигается вправо на 21, т.е.

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

И, наконец, вы можете разложить любое умножение на несколько добавлений, посмотрев на двоичный 1 в множителе. Поэтому вы разрешаете "масштабирование", и я предположил, что это означает умножение. Если нет, вот 32-битная версия (64-разрядная версия оставлена ​​в качестве упражнения для читателя), использующей тот факт, что 629146 - 10011001100110011010 (аккуратный шаблон из-за повторяющейся двоичной дроби). Мы обойдемсь в другую сторону и вместо этого используем 10011001100110011001.

/* 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;
}

Ответ 2

Если у вас есть быстрое целочисленное умножение, но у вас нет целочисленного деления, вы можете получить разумное приближение, умножив его на 1229, а затем сдвигая на 12 бит. Например:

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

Это работает, потому что 1229 составляет около 0.3 * 1 << 12. "Реальное" значение составляет 1228,8, поэтому в некоторых случаях оценка будет высокой на 1 (68 из 4097 значений). Это никогда не будет выключено более чем на 1, однако.

Ответ 3

Использование классического взлома для divu10() ниже. Я предпочитаю подход * n и shift, но думал, что предлагаю еще один 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);  
}

Ответ 4

Я знаю, что это не пример кодирования, извините за это, но если набор настолько мал (range [0; 4096]), почему бы не создать блок результатов и не использовать указатель для извлечения значений? Это значительно сократит циклы GPU, так как ограничений памяти нет.

Ответ 5

Просто попробовал много комбинаций формы (A*x + B) >> n с приведенным ниже тестовым кодом и придумал:

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

Тестовый код

#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;
}