Переполнение при умножении

Перед выполнением умножения C++ приводит множители к одному типу не короче int, а разрядность результата совпадает с разрядностью приведенных множителей. Для того, чтобы не потерять точность, иногда требуется для умножения выполнять дополнительные операции.
Рассмотрим задачу. Система определяет время с момента запуска программы в тиках микропроцессора вызовом функции:
unsigned long long getTickCount();

Длина unsigned long long — 64 бита. Для преобразования в физические единицы времени в системе существует константа:
const unsigned long long TICKS_PER_SECOND = 1999000001ULL;

Требуется определить функцию перевода тиков в наносекунды getNsec(unsigned long long ticks) с семантикой:
unsigned long long getNsecNaive(unsigned long long ticks) {
    static const unsigned long long NSEC_PER_SECOND = 1000000000ULL;
    unsigned long long nsec =  NSEC_PER_SECOND * ticks / TICKS_PER_SECOND;
    return nsec;
}

Для функции getNsec() необходимо обеспечить максимально возможную точность. Параметр ticks может быть как большим, так и маленьким. Для маленьких ticks (скажем, до 2^34) нужно вычисления выполнять в следующем порядке:
(NSEC_PER_SECOND * ticks) / TICKS_PER_SECOND

То есть сначала умножать, потом делить. Поскольку NSEC_PER_SECOND < 2^30, при умножении не возникнет переполнения, так как его результат будет меньше 2^64.
Для больших ticks, поскольку умножение может вызвать переполнение, лучше выполнять операции в следующем порядке:
NSEC_PER_SECOND * (ticks / TICKS_PER_SECOND)

Проблема в данном случае состоит в том, что второй множитель — всегда целый, результат в десятичной системе счисления всегда будет оканчиваться девятью нулями, то есть фактически обеспечивается секундная точность и функцию getNsec() следует переименовать в getSec(). С другой стороны, исходные данные позволяют обеспечить большую точность.
С учетом неассоциативности машинных операций умножения и деления всего возможны 4 порядка вычислений, то есть еще 2 помимо двух рассмотренных:
(NSEC_PER_SECOND / TICKS_PER_SECOND) * ticks

и
ticks / (TICKS_PER_SECOND / NSEC_PER_SECOND)

Первый всегда дает ноль, а второй — ticks, то есть почти 50% ошибку (которую в данном случае можно свести до 0,1%, но эта погрешность все равно не наименьшая возможная).
Итак, в лучшем случае получаем секундную точность. Решить проблему увеличения точности можно следующими способами (в порядке увеличения предпочтения):
  1. Производить вычисления с вещественными (double) типами
  2. Приводить аргументы к 128-разрядным целым
  3. Использовать функцию MultDiv()
Данные подходы могут быть не применимы по причинам ограничений платформы (процессора), среды программирования и библиотек, требований проекта.
Воспользуемся следующим подходом. Пусть в результате деления ticks на TICKS_PER_SECOND получаем неполное частное seconds и остаток r:
unsigned long long seconds = ticks / TICKS_PER_SECOND; 
unsigned long long r = ticks % TICKS_PER_SECOND;

Тогда ticks = seconds*TICKS_PER_SECOND + r. Подставляем в формулу для nsec:
nsec = NSEC_PER_SECOND * (seconds*TICKS_PER_SECOND + r) / TICKS_PER_SECOND = NSEC_PER_SECOND*seconds + (NSEC_PER_SECOND*r) / TICKS_PER_SECOND. Поскольку r < TICKS_PER_SECOND, (NSEC_PER_SECOND*r) никогда не даст переполнения. Итоговая функция:
unsigned long long getNsec(unsigned long long ticks) {
    static const unsigned long long NSEC_PER_SECOND = 1000000000ULL;
    return NSEC_PER_SECOND * (ticks / TICKS_PER_SECOND)
        + ((NSEC_PER_SECOND)*(ticks % TICKS_PER_SECOND))/TICKS_PER_SECOND;
}

Результат получается ровно таким, как если бы мы посчитали NSEC_PER_SECOND*ticks как 128-битное значение, а потом разделили на TICKS_PER_SECOND, то есть обеспечивается максимальную точность для данных исходных значений и данной разрядности представления результата. Формула в операторе return никак не упрощается: ни NSEC_PER_SECOND, ни /TICKS_PER_SECOND за скобки не выносится.
Фактически решение задачи сведено к реализации функции MultDiv(a, b, c), вычисляющей a*b/c, где b и c — большие константы, а дробь b/c не сократима.
Поделиться публикацией

Комментарии 23

    +10
    Есть классная книжка "Алгоритмические трюки для программистов" — в ней ещё много полезных приёмов по целочисленному умножению и делению.
      0
      Есть у меня эта книга, только не всю прочитал. Там есть конкретно такой пример? Я не претендую на право первооткрывателя, но приведенный код придуман самостоятельно, хотя он и не сложный, а предварительное гугление результатов не дало
        +2
        Олег, вы неправильно меня поняли. Я коммент написал не чтобы уличить Вас в переизобретении велосипеда, а к тому, что есть ещё и другие интересные приёмы, про которые полезно почитать, и которые могут часто пригодиться.
      +2
      результат умножения имеет ту же разрядность, что и множители.

      Да ну? Какого типа тогда будет результат умножения unsigned char на unsigned char?
        +1
        Все арифметические типы Уже int перед арифметикой приводятся к типу int. Т.е. unsigned char * unsigned char = int.
          0
          Согласен. Нужно было написать что-то вроде:
          «результат умножения целых типа не короче int имеет ту же разрядность, что и множители»
            0
            Еще точнее: «C++ не умеет умножать целые короче int, а для целых не короче int разрядность результата совпадает с разрядностью умножения». Для long тоже можно изобрести что-то вроде:
            long a, b;
            //...
            long long c = a*static_cast<long long>(b);

            Но это уже выходит за рамки условий поста. Там задача — остаться в рамках исходной разрядности.
          +1
          > В C/C++, в отличие от большинства ассемблеров,
          > результат умножения имеет ту же разрядность, что и множители.

          unsigned short a = 60000;
          unsigned short b = 60000;
          unsigned int i = a*b;

          А так?
          • НЛО прилетело и опубликовало эту надпись здесь
              0
              4.5 Integral promotions — это часть «usual arithmetic conversions».
              0
              Вы правы, смотрите мой ответ на предыдущий комментарий
              –7
              Модельный пример несколько высосан из пальца, ибо гораздо лучшим решением будет использование функции QueryPerformanceCounter.
              • НЛО прилетело и опубликовало эту надпись здесь
                  0
                  LARGE_INTEGER c;
                  QueryPerformanceCounter(&c);
                  return c.QuadPart;

                  а в чем должна была быть проблема? этот фрагмент, конечно, забывает про секунды, но если тип результата позволит их вернуть, то можно просто прибавить их, умножив на 1000000000. разница в том, что здесь нет деления и оно не нужно вообще.
                    0
                    нет, все же наврал. делить нужно на соответствующее значение QueryPerformanceFrequency. отличие разве что в том, что делитель здесь небольшой.
                  +1
                  … особенно, где-нить на embedded системе вообще без ОС :)
                    +1
                    Нет, моедльный пример — в реально работающем проекте. Вот только операционка там — ядро линукс, а все драйвера — свои. Из сторонних библиотек только libc
                    –1
                    Не проще ли один раз посчитать TICKS_PER_NSEC?
                      +1
                      И чему оно будет равно? Обратите внимание, что в соответствии с условием вещественные типы использовать нельзя
                        0
                        unsigned long long getNsec(unsigned long long ticks) {
                            static const unsigned long long _GCD_TPS_NSPS = gcd(NSEC_PER_SECOND, TICKS_PER_SECOND);
                            return ticks * (NSEC_PER_SECOND / _GCD_TPS_NSPS) / (TICKS_PER_SECOND / _GCD_TPS_NSPS);
                        };
                        


                          0
                          Я для упрощения округлил до 1999000000ULL, полагая, что значение роли не играет. Считайте gcd = 1, такое будет, например, для TICKS_PER_SECOND = 1999000001ULL;
                            0
                            Сделал соответствующее исправление в тексте, чтобы не возникало дальнейших вопросов с GCD. Стоит также заметить, что в исходной задаче константы были взаимно просты, и лишь при подготовке начального варианта поста я произвел округление, которое по моему недогляду привело к еще одному решению через GCD.
                      0
                      Все корректно. Кое где можно было бы привести умножение вместо второй операции '%' — может съэкономить пару тактов, если часто вызываема будет, поскольку операция деления и взятия остатка не сможет быть даже гипотетически объединена в одну операцию конвеером процессора (существуют в разрыв им другие операции). Нахождение общего делителя интересный момент, но гораздо хуже, если коэффициенты будут иметь единицу в качестве онного.

                      Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                      Самое читаемое