Блог

Предположим, у меня есть некоторый код, такой как:

 float a, b = ...; // both positive
int s1 = ceil(sqrt(a/b));
int s2 = ceil(sqrt(a/b))   0.1;
  

Возможно ли это когда-нибудь s1 != s2 ? Меня беспокоит, когда a/b это идеальный квадрат. Например, возможно a=100.0 , и b=4.0 , тогда вывод ceil должен быть 5.00000 , но что, если вместо этого это 4.99999 так?

Аналогичный вопрос: есть ли вероятность, что 100.0/4.0 вычисляется значение say 5.00001 , а затем ceil округляется до 6.00000 ?

Я бы предпочел сделать это в целочисленной математике, но sqrt это своего рода план.

РЕДАКТИРОВАТЬ: предложения о том, как лучше реализовать это, также будут оценены! Значения a и b являются целочисленными значениями, поэтому фактический код больше похож на: ceil(sqrt(float(a)/b))

РЕДАКТИРОВАТЬ: основываясь на ответе levis501, я думаю, что сделаю это:

 float a, b = ...; // both positive
int s = sqrt(a/b);
while (s*s*b < a)   s;
  

Спасибо вам всем!

Я не думаю, что это возможно. Независимо от значения sqrt(a/b) , оно выдает некоторое значение N, которое мы используем как:

 int s1 = ceil(N);
int s2 = ceil(N)   0.1;
  

Поскольку ceil всегда выдает целочисленное значение (хотя и представленное как double ), у нас всегда будет некоторое значение X, для которого первое производит X.0 и второе X.1 . Преобразование в int всегда будет усекать это .1 , поэтому оба приведут к X .

Может показаться, что было бы исключение, если бы X было настолько большим, что X.1 превысил диапазон double . Я не вижу, где это может быть возможно. За исключением случаев, близких к 0 (где переполнение не вызывает беспокойства), квадратный корень из числа всегда будет меньше входного числа. Следовательно, до того, как ceil(N) 0.1 может переполниться, значение, a/b используемое в качестве входных sqrt(a/b) данных, должно было бы уже переполниться.

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

 /* return the smallest positive integer whose square is at least x */
int isqrt(double x) {
  int y1 = ceil(sqrt(x));
  int y2 = y1 - 1;
  if ((y2 * y2) >= x) return y2;
  return y1;
}
  

Это будет обрабатывать нечетный случай, когда квадратный корень из вашего отношения a/b находится в пределах точности double .

Равенство чисел с плавающей запятой действительно является проблемой, но ИМХО нет, если мы имеем дело с целыми числами.

Если у вас есть случай 100.0/4.0 , он должен идеально вычисляться 25.0 , поскольку 25.0 точно представим как значение с плавающей точкой, в отличие от, например 25.1 .

Да, вполне возможно, что s1 != s2 . Но почему это проблема? Кажется достаточно естественным, что s1 != (s1 0.1) .

Кстати, если вы предпочитаете 5.00001 округлять до 5.00000 вместо 6.00000 , используйте rint вместо ceil .

И чтобы ответить на актуальный вопрос (в вашем комментарии) — вы можете использовать sqrt , чтобы получить отправную точку, а затем просто найти правильный квадрат, используя целочисленную арифметику.

 int min_dimension_greater_than(int items, int buckets)
{
    double target = double(items) / buckets;
    int min_square = ceil(target);
    int dim = floor(sqrt(target));
    int square = dim * dim;
    while (square < min_square) {
        seed  = 1;
        square = dim * dim;
    }
    return dim;
}
  

И да, это можно значительно улучшить, это всего лишь краткий набросок.

s1 всегда будет равно s2 .

Стандарты C и C мало что говорят о точности математических процедур. В буквальном смысле невозможно реализовать стандарт, поскольку стандарт C гласит, что sqrt(x) возвращает квадратный корень из x, но квадратный корень из двух не может быть точно представлен с плавающей запятой.

Реализация подпрограмм с хорошей производительностью, которые всегда возвращают правильно округленный результат (в режиме округления до ближайшего это означает, что результатом является представимое число с плавающей запятой, ближайшее к точному результату, с разрешенными связями в пользу младшего нулевого бита), является сложной исследовательской проблемой. Хорошие математические библиотеки нацелены на точность менее 1 ULP (поэтому возвращается одно из двух ближайших представимых чисел), возможно, что-то чуть больше .5 ULP. (ULP — это единица наименьшей точности, значение младшего бита, заданное конкретным значением в поле экспоненты.) Некоторые математические библиотеки могут быть значительно хуже этого. Вам нужно будет спросить своего поставщика или проверить документацию для получения дополнительной информации.

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

Аналогично, если библиотека гарантирует, что ошибки меньше 1 ULP, то ceil должен вернуть точный результат, опять же, потому что точный результат представим, а любой другой результат будет удален как минимум на 1 ULP. Кроме того, природа ceil такова, что я ожидаю, что любая разумная математическая библиотека всегда будет возвращать целое число, даже если остальная часть библиотеки не была высокого качества.

Что касается случаев переполнения, если ceil(x) находится за пределами диапазона, в котором все целые числа точно представимы, то ceil(x) .1 ближе к ceil(x), чем к любому другому представимому числу, поэтому округленный результат добавления .1 к ceil(x) должен быть ceil(x) в любой системе, реализующей стандарт с плавающей запятой (IEEE 754). Это при условии, что вы находитесь в режиме округления по умолчанию, который является округлением до ближайшего. Можно изменить режим округления на что-то вроде округления до бесконечности, что может привести к тому, что ceil(x) .1 будет целым числом выше, чем ceil(x) .