Разделить диапазон на равные интервалы - PullRequest
Новая
       44

Разделить диапазон на равные интервалы

0 голосов
/ 21 октября 2018

Я хочу разбить диапазон с double границами на N>=2 равные или почти равные интервалы.

Я нашел подходящую функцию в Научная библиотека GNU : 

make_uniform (double range[], size_t n, double xmin, double xmax)
{
  size_t i;

  for (i = 0; i <= n; i++)
    {
      double f1 = ((double) (n-i) / (double) n);
      double f2 = ((double) i / (double) n);
      range[i] = f1 * xmin +  f2 * xmax;
    }
}

Однако, когда
xmin = 241141 (двоичный 0x410D6FA800000000)
xmax = 241141.0000000001 (двоичный 0x410D6FA800000003)
N = 3
, функция выдает 

[0x410D6FA800000000,
 0x410D6FA800000000,
 0x410D6FA800000002,
 0x410D6FA800000003]

вместо желаемого 

[0x410D6FA800000000,
 0x410D6FA800000001,
 0x410D6FA800000002,
 0x410D6FA800000003]

Как добиться однородности, не прибегая к длинной арифметике (у меня уже есть решение с длинной арифметикой, но оно уродливо и медленно)?Подходящие битовые операции и x86 (x86-64, поэтому нет повышенной точности) являются приемлемыми.

ОБНОВЛЕНИЕ:

Требуется общее решение, без предпосылки, что xmin, xmax имеют равныеэкспонента и знак:

  • xmin и xmax могут иметь любое значение, кроме бесконечности и NaN (возможно, также исключая денормализованные значения для простоты).
  • xmin < xmax
  • (1<<11)-1>=N>=2
  • Я готов к значительным (в 2-3 заказах) потерям производительности

Ответы [ 2 ]

0 голосов
/ 23 октября 2018

x87 все еще существует в x86-64, и 64-разрядные ядра для основных операционных систем правильно сохраняют / восстанавливают состояние x87 для 64-разрядных процессов.Несмотря на то, что вы, возможно, прочитали, x87 полностью применим в 64-битном коде.

Вне Windows (т. Е. X86-64 System V ABI, который используется везде), long double - это 80-битный собственный x87родной формат.Это, вероятно, решит вашу проблему точности только для x86 / x86-64, если вы не заботитесь о переносимости в ARM / PowerPC / что-то еще, что имеет только 64-битную точность в HW.

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

Я не уверен, что вам нужно делать в Windows, чтобы получить компилятор для генерации 80-битной расширенной математики FP.Это, конечно, возможно в asm и поддерживается ядром, но набор инструментов и ABI делают его неудобным для использования.

x87 только несколько медленнее, чем скалярная математика SSE на современных процессорах.Однако 80-битная загрузка / хранение выполняется слишком медленно, например, 4 мопа на Skylake вместо 1 (https://agner.org/optimize/) и несколько циклов дополнительной задержки для fld m80.

Для цикла, в котором необходимо преобразоватьint в FP, сохраняя и используя x87 fild, это может быть примерно в 2 раза медленнее, чем то, что хороший компилятор может сделать с SSE2 для 64-битного двойного.

И, конечно, long double предотвратит авто-векторизацию.

0 голосов
/ 21 октября 2018

Я вижу два варианта: переупорядочить операции как xmin + (i * (xmax - xmin)) / n или иметь дело непосредственно с двоичными представлениями.Вот пример для обоих. 

#include <iostream>
#include <iomanip>

int main() {
    double xmin = 241141;
    double xmax = 241141.0000000001;
    size_t n = 3, i;
    double range[4];

    std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::fixed;

    for (i = 0; i <= n; i++) {
        range[i] = xmin + (i * (xmax - xmin)) / n;

        std::cout << range[i] << "\n";
    }
    std::cout << "\n";

    auto uxmin = reinterpret_cast<unsigned long long&>(xmin);
    auto uxmax = reinterpret_cast<unsigned long long&>(xmax);

    for (i = 0; i <= n; i++) {
        auto rangei = ((n-i) * uxmin + i * uxmax) / n;
        range[i] = reinterpret_cast<double&>(rangei);

        std::cout << range[i] << "\n";
    }
}

Live на Coliru

...