Задача с погрешностью и переполнением

Хочу с вами поделиться небольшой алгоритмической задачей, которую мне не так давно подкинули коллеги по цеху. Как мне кажется, решение получилось довольно элегантное.

Суть задачи заключалась в следующем. Есть некое устройство (микроконтроллер), которое умеет обращаться только с 32-битными целыми значениями и не умеет работать с плавающей запятой.
На таком аппарате есть таймер, который в секунду генерирует 32768 тиков. Необходимо написать функцию, переводящую тики в миллисекунды без потери точности (желательно с округлением).
Примечание: все примеры в статье сделаны на Java, но это не должно помешать пониманию алгоритмов.

Очевидно, что решение «в лоб» — m = (ticks * 1000) / (2^15) не подходит, так как может легко привести к переполнению.

Эврика?

Первым вариантом, о котором говорили сотрудники — тупо делить на 32, но здесь получается большая погрешность.
Другим вариантов решения предполагалось «умное» деление с умножением, в несколько этапов с серией проверок, стараясь на каждом шаге не переполнить 32 бита.

Достаточно быстро мне в голову пришла следующая цепочка преобразований (помня, что 1000 = 1024 - 24 ):
ticks * 1000 * 2^(-15) =
ticks * 2^(-5) - ticks * 24 * 2^(-15) = ticks * 2^(-5) - 3 * ticks * 2^(-12)

Итого имеем:

	private int getMillis(int ticks) {
		return (ticks >> 5) - 3*(ticks >> 12);
	}

Супер! Жаль, что это не работает.

Первое приближение

Подумал, проверил лишний раз — очевидных ошибок в рассуждении нет. Проверил на нескольких тестовых примерах — такой алгоритм считает, но не всегда корректно. В ряде случаев есть отклонения от правильного результата на +-3. Судя по всему, в математике проблем нет (иначе ответы бы не были похожи на правду), но имеется погрешность.

Чтобы понять, в чем проблема, стоит вспомнить, как работает битовый сдвиг. Операция битового сдвига вправо на n соответствует делению нацело на 2 в степени n. Это целочисленная операция, поэтому биты, оказавшиеся «правее» нулевого бита просто отбрасываются.

Если присмотреться к предложенному алгоритму, мы имеем 2 операции сдвига, умножение одного из результатов и разность. Схематично эта операция изображена ниже.

Очевидно, что погрешность может накапливаться именно в битах, оказавшихся правее нулевого бита.

Первоначально у меня пало подозрение на последний отброшенный (минус-первый) бит, возможно именно его не хватало в операциях умножения и разности. В результате код функции стал несколько сложнее:

	private int getMillis(int ticks) {
		int err = ((ticks & 16) >> 4) - 3 * ((ticks & 2048) >> 11);
		if (err < 0) {
			return (ticks >> 5) - 3 * (ticks >> 12) - (Math.abs(err) >> 1);
		} else {
			return (ticks >> 5) - 3 * (ticks >> 12) + (err >> 1);
		}
	}

В переменную err мы помещаем тот бит, который был отброшен при операции сдвига, и над ним отдельно проводится та же операция вида x - 3 * y. Для этого выделяются соответственно 4-й и 11-й биты, сдвигаются в нулевую позицию, вычисляется операция и результат еще раз сдвигается вправо на 1 бит.
В результате все, что осталось после сдвига — это потерянные биты, которые мы используем для корректировки результата.
Результат улучшился: теперь абсолютная погрешность составляла не более 2.

Решение

Теперь стало очевидно, что направление мысли было правильное, но участвовать должен не один бит правее запятой, а все 5 пересекающихся (то есть которые оказываются в позиции от -1 до -5 на приведенной схеме).

В результате я получил окончательный вариант функции, погрешность которой составляла менее 0.5.
Общая идея была та же, что в предыдущем примере, но в качестве источника погрешности выделяется не по 1 биту из каждого операнда, а по 5.

Кроме этого, реализовано округление, которое основано на идее, что i-й бит равен половине i+1-го.

Все остальное должно быть понятно из кода :)

	private int getMillis(int ticks) {
		int l = (ticks & 0x1F) << 7;
		int r = ticks & 0xFFF; // 2^12 - 1

		int err = l - 3 * r;
		int fraction = (Math.abs(err) & 0xFFF);

		int fractionRound = fraction >> 11; //== 1 <=> there is 0.5xxx
		fraction = fraction & 0x7FF; // 0.5xxx => 0.0xxx

		if (err < 0) {
			fraction = fraction > 0 ? 1 : 0;
			return (ticks >> 5) - 3 * (ticks >> 12) - (Math.abs(err) >> 12)
					- (fractionRound & fraction);
		} else {
			return (ticks >> 5) - 3 * (ticks >> 12) + (err >> 12)
					+ fractionRound;
		}
	}

Будьте хорошими программистами и приятных вам выходных!