Оценка умножения с экспоненциальной функцией

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

Поскольку значение примера c=0.95 близко к 1, любой корень будет еще ближе к 1 и, таким образом, потеряет точность. Избегайте этого, вычисляя разницу в 1 напрямую, z=1-c^(1/n) через

z = -expm1(log(c)/N).

Теперь вам нужно вычислить многочлен

sum of x[i] * (1-z)^i

что можно сделать путем тщательной модификации схемы Горнера. Вместо

for(i=N; i-->0; ) {
  res = res*(1-z)+x[i]
}

использовать

for(i=N; i-->0; ) {
  res = (res+x[i])-res*z
}

что математически эквивалентно, но потеря цифр в 1-z происходит как можно позже без использования более сложного метода, такого как сложение с двойной точностью.

В тестах эти два метода, вопреки замыслу, дали почти одинаковые результаты, существенное улучшение можно было наблюдать, разделив результат на его значение в c=1, z=0 и кратное z, как в

double res0 = 0, resz=0;
int i;
for(i=N; i-->0; ) { 
    /* res0+z*resz = (res0+z*resz)*(1-z)+x[i]; */
    resz = resz - res0 -z*resz;
    res0 = res0 + x[i];
}

Тестовый пример, который показал это улучшение, был для последовательности коэффициентов

f(u) = (1-u/N)^(N-2)*(1-u)

где для N=1000 оценки приводят к

   c               z=1-c^(1/N)             f(1-z)         diff for 1st proc     diff for 3rd proc

  0.950000     0.000051291978909      0.000018898570629  1.33289104579937e-17  4.43845264361253e-19
  0.951000     0.000050239954368      0.000018510931892  1.23765066121009e-16  -9.24959978401696e-19
  0.952000     0.000049189034371      0.000018123700958  1.67678642238461e-17  -5.38712954453735e-19
  0.953000     0.000048139216599      0.000017736876972  -2.86635949350855e-17  -2.37169225231204e-19
...
  0.994000     0.000006018054217      0.000002217256601  1.31645860662263e-17  1.15619997300212e-19
  0.995000     0.000005012529261      0.000001846785028  -4.15668713370839e-17  -3.5363625547867e-20
  0.996000     0.000004008013365      0.000001476685973  8.48811716443534e-17    8.470329472543e-22
  0.997000     0.000003004504507      0.000001106958687  1.44711343873661e-17  -2.92226366802734e-20
  0.998000     0.000002002000667      0.000000737602425   5.6734266807093e-18  -6.56450534122083e-21
  0.999000     0.000001000499833      0.000000368616443  -3.72557383333555e-17  1.47701370177469e-20

Ответ Ива вдохновил меня.

Кажется, что лучший подход — вычислять pow(c, 1.0/N) не напрямую, а косвенно:

cc[0]=c; cc[1]=sqrt(cc[0]), cc[2]=sqrt(cc[1]),... cc[logN] = sqrt(cc[logN-1])

Или в двоичном формате,

cc[0]=c, cc[1]=c^0.1, cc[2]=c^0.01, cc[3]=c^0.001, ....

Теперь, если нам нужно x[0b100100] * c^0.100100, мы можем рассчитать это как x[0b100100]* c^0.1 * c^0.0001. Мне не нужно предварительно вычислять таблицу размера N, как предложил geza. Таблицы размера log(N), вероятно, достаточно, и ее можно создать, многократно извлекая квадратные корни.

[править] Как указано в ветке комментариев к другому ответу, попарное суммирование очень эффективно для контроля над ошибками. И это очень хорошо сочетается с этим ответом.

Начнем с наблюдения, что мы суммируем

x[0] * c^0.0000000
x[1] * c^0.0000001
x[2] * c^0.0000010
x[3] * c^0.0000011
...

Итак, мы запускаем log(N) итераций. На итерации 1 мы добавляем N/2 пар x[i]+x[i+1]*c^0.000001 и сохраняем результат в x[i/2]. На итерации 2 мы добавляем пары x[i]+x[i+1]*c^0.000010 и так далее. Основное отличие от обычного попарного суммирования состоит в том, что это умножение и сложение на каждом шаге.

Теперь мы видим, что на каждой итерации мы используем один и тот же множитель pow(c, 2^i/N), что означает, что нам нужно вычислить только множители log(N). Это также довольно эффективно с точки зрения кэширования, поскольку мы выполняем доступ только к непрерывной памяти. Это также позволяет легко распараллеливать SIMD, особенно когда у вас есть инструкции FMA.

Если N является степенью 2, вы можете заменить оценки степеней геометрическими средствами, используя

a^(i+j)/2 = √(a^i.a^j)

и рекурсивно разделить от c^N/N.c^0/N. С рекурсией предварительного заказа вы можете накапливать, увеличивая веса.

В любом случае, ускорение sqrt по сравнению с pow может быть незначительным.

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

Вы можете смешать многократное умножение на pow(c, 1./N) с некоторыми явными вызовами pow. т.е. каждая 16-я итерация или около того делает реальное pow и в противном случае продвигается вперед с умножением. Это должно дать большие преимущества в производительности при незначительной стоимости точности.

В зависимости от того, насколько сильно варьируется c, вы даже можете предварительно вычислить и заменить все вызовы pow поиском или только те, которые необходимы в приведенном выше методе (= меньшая таблица поиска = лучшее кэширование).