Расчет Fast Exp: можно повысить точность, не теряя при этом слишком высокой производительности?

Ответ 1

В случае, если кто-то хочет воспроизвести функцию относительной ошибки, показанную в вопросе, здесь способ использования Matlab ( "быстрый" показатель не очень быстр в Matlab, но он точен):

t = 1072632447+[0:ceil(1512775*pi)];
x = (t - 1072632447)/1512775;
ex = exp(x);
t = uint64(t);
import java.lang.Double;
et = arrayfun( @(n) java.lang.Double.longBitsToDouble(bitshift(n,32)), t );
plot(x, et./ex);

Теперь период ошибки точно совпадает с тем, когда двоичное значение tmp переполняется из мантиссы в экспоненту. Позвольте разбить наши данные на бины, отбросив биты, которые становятся экспонентами (делая их периодическими), и сохраняя только восемь оставшихся бит (чтобы сделать нашу таблицу поиска разумным):

index = bitshift(bitand(t,uint64(2^20-2^12)),-12) + 1;

Теперь вычислим среднюю требуемую настройку:

relerrfix = ex./et;
adjust = NaN(1,256);
for i=1:256; adjust(i) = mean(relerrfix(index == i)); end;
et2 = et .* adjust(index);

Относительная ошибка уменьшается до +/-.0006. Конечно, возможны и другие размеры таблиц (например, 6-битная таблица с 64 записями дает +/-.0025), и ошибка почти линейна в размере таблицы. Линейная интерполяция между табличными элементами еще более улучшила бы ошибку, но за счет производительности. Поскольку мы уже достигли цели точности, позвольте избежать дальнейших поражений производительности.

На этом этапе некоторые тривиальные навыки редактора принимают значения, вычисленные MatLab, и создают таблицу поиска на С#. Для каждого вычисления мы добавляем битмаску, поиск таблицы и умножение с двойной точностью.

static double FastExp(double x)
{
    var tmp = (long)(1512775 * x + 1072632447);
    int index = (int)(tmp >> 12) & 0xFF;
    return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}

Ускорение очень похоже на исходный код - для моего компьютера это примерно на 30% быстрее компилируется как x86 и примерно в 3 раза быстрее для x64. С моно на идеоне это существенный чистый убыток (но это оригинал).

Полный исходный код и тестовый файл: http://ideone.com/UwNgx

using System;
using System.Diagnostics;

namespace fastexponent
{
    class Program
    {
        static double[] ExpAdjustment = new double[256] {
            1.040389835,
            1.039159306,
            1.037945888,
            1.036749401,
            1.035569671,
            1.034406528,
            1.033259801,
            1.032129324,
            1.031014933,
            1.029916467,
            1.028833767,
            1.027766676,
            1.02671504,
            1.025678708,
            1.02465753,
            1.023651359,
            1.022660049,
            1.021683458,
            1.020721446,
            1.019773873,
            1.018840604,
            1.017921503,
            1.017016438,
            1.016125279,
            1.015247897,
            1.014384165,
            1.013533958,
            1.012697153,
            1.011873629,
            1.011063266,
            1.010265947,
            1.009481555,
            1.008709975,
            1.007951096,
            1.007204805,
            1.006470993,
            1.005749552,
            1.005040376,
            1.004343358,
            1.003658397,
            1.002985389,
            1.002324233,
            1.001674831,
            1.001037085,
            1.000410897,
            0.999796173,
            0.999192819,
            0.998600742,
            0.998019851,
            0.997450055,
            0.996891266,
            0.996343396,
            0.995806358,
            0.995280068,
            0.99476444,
            0.994259393,
            0.993764844,
            0.993280711,
            0.992806917,
            0.992343381,
            0.991890026,
            0.991446776,
            0.991013555,
            0.990590289,
            0.990176903,
            0.989773325,
            0.989379484,
            0.988995309,
            0.988620729,
            0.988255677,
            0.987900083,
            0.987553882,
            0.987217006,
            0.98688939,
            0.98657097,
            0.986261682,
            0.985961463,
            0.985670251,
            0.985387985,
            0.985114604,
            0.984850048,
            0.984594259,
            0.984347178,
            0.984108748,
            0.983878911,
            0.983657613,
            0.983444797,
            0.983240409,
            0.983044394,
            0.982856701,
            0.982677276,
            0.982506066,
            0.982343022,
            0.982188091,
            0.982041225,
            0.981902373,
            0.981771487,
            0.981648519,
            0.981533421,
            0.981426146,
            0.981326648,
            0.98123488,
            0.981150798,
            0.981074356,
            0.981005511,
            0.980944219,
            0.980890437,
            0.980844122,
            0.980805232,
            0.980773726,
            0.980749562,
            0.9807327,
            0.9807231,
            0.980720722,
            0.980725528,
            0.980737478,
            0.980756534,
            0.98078266,
            0.980815817,
            0.980855968,
            0.980903079,
            0.980955475,
            0.981017942,
            0.981085714,
            0.981160303,
            0.981241675,
            0.981329796,
            0.981424634,
            0.981526154,
            0.981634325,
            0.981749114,
            0.981870489,
            0.981998419,
            0.982132873,
            0.98227382,
            0.982421229,
            0.982575072,
            0.982735318,
            0.982901937,
            0.983074902,
            0.983254183,
            0.983439752,
            0.983631582,
            0.983829644,
            0.984033912,
            0.984244358,
            0.984460956,
            0.984683681,
            0.984912505,
            0.985147403,
            0.985388349,
            0.98563532,
            0.98588829,
            0.986147234,
            0.986412128,
            0.986682949,
            0.986959673,
            0.987242277,
            0.987530737,
            0.987825031,
            0.988125136,
            0.98843103,
            0.988742691,
            0.989060098,
            0.989383229,
            0.989712063,
            0.990046579,
            0.990386756,
            0.990732574,
            0.991084012,
            0.991441052,
            0.991803672,
            0.992171854,
            0.992545578,
            0.992924825,
            0.993309578,
            0.993699816,
            0.994095522,
            0.994496677,
            0.994903265,
            0.995315266,
            0.995732665,
            0.996155442,
            0.996583582,
            0.997017068,
            0.997455883,
            0.99790001,
            0.998349434,
            0.998804138,
            0.999264107,
            0.999729325,
            1.000199776,
            1.000675446,
            1.001156319,
            1.001642381,
            1.002133617,
            1.002630011,
            1.003131551,
            1.003638222,
            1.00415001,
            1.004666901,
            1.005188881,
            1.005715938,
            1.006248058,
            1.006785227,
            1.007327434,
            1.007874665,
            1.008426907,
            1.008984149,
            1.009546377,
            1.010113581,
            1.010685747,
            1.011262865,
            1.011844922,
            1.012431907,
            1.013023808,
            1.013620615,
            1.014222317,
            1.014828902,
            1.01544036,
            1.016056681,
            1.016677853,
            1.017303866,
            1.017934711,
            1.018570378,
            1.019210855,
            1.019856135,
            1.020506206,
            1.02116106,
            1.021820687,
            1.022485078,
            1.023154224,
            1.023828116,
            1.024506745,
            1.025190103,
            1.02587818,
            1.026570969,
            1.027268461,
            1.027970647,
            1.02867752,
            1.029389072,
            1.030114973,
            1.030826088,
            1.03155163,
            1.032281819,
            1.03301665,
            1.033756114,
            1.034500204,
            1.035248913,
            1.036002235,
            1.036760162,
            1.037522688,
            1.038289806,
            1.039061509,
            1.039837792,
            1.040618648
        };

        static double FastExp(double x)
        {
            var tmp = (long)(1512775 * x + 1072632447);
            int index = (int)(tmp >> 12) & 0xFF;
            return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
        }

        static void Main(string[] args)
        {
            double[] x = new double[1000000];
            double[] ex = new double[x.Length];
            double[] fx = new double[x.Length];
            Random r = new Random();
            for (int i = 0; i < x.Length; ++i)
                x[i] = r.NextDouble() * 40;

            Stopwatch sw = new Stopwatch();
            sw.Start();
            for (int j = 0; j < x.Length; ++j)
                ex[j] = Math.Exp(x[j]);
            sw.Stop();
            double builtin = sw.Elapsed.TotalMilliseconds;
            sw.Reset();
            sw.Start();
            for (int k = 0; k < x.Length; ++k)
                fx[k] = FastExp(x[k]);
            sw.Stop();
            double custom = sw.Elapsed.TotalMilliseconds;

            double min = 1, max = 1;
            for (int m = 0; m < x.Length; ++m) {
                double ratio = fx[m] / ex[m];
                if (min > ratio) min = ratio;
                if (max < ratio) max = ratio;
            }

            Console.WriteLine("minimum ratio = " + min.ToString() + ", maximum ratio = " + max.ToString() + ", speedup = " + (builtin / custom).ToString());
         }
    }
}

Ответ 2

Попробуйте следующие альтернативы (exp1 быстрее, exp7 более точно).

код

public static double exp1(double x) { 
    return (6+x*(6+x*(3+x)))*0.16666666f; 
}

public static double exp2(double x) {
    return (24+x*(24+x*(12+x*(4+x))))*0.041666666f;
}

public static double exp3(double x) {
    return (120+x*(120+x*(60+x*(20+x*(5+x)))))*0.0083333333f;
}

public static double exp4(double x) {
    return 720+x*(720+x*(360+x*(120+x*(30+x*(6+x))))))*0.0013888888f;
}

public static double exp5(double x) {
    return (5040+x*(5040+x*(2520+x*(840+x*(210+x*(42+x*(7+x)))))))*0.00019841269f;
}

public static double exp6(double x) {
    return (40320+x*(40320+x*(20160+x*(6720+x*(1680+x*(336+x*(56+x*(8+x))))))))*2.4801587301e-5;
}

public static double exp7(double x) {
  return (362880+x*(362880+x*(181440+x*(60480+x*(15120+x*(3024+x*(504+x*(72+x*(9+x)))))))))*2.75573192e-6;
}

Точность

Function     Error in [-1...1]              Error in [3.14...3.14]

exp1         0.05           1.8%            8.8742         38.40%
exp2         0.01           0.36%           4.8237         20.80%
exp3         0.0016152      0.59%           2.28            9.80%
exp4         0.0002263      0.0083%         0.9488          4.10%
exp5         0.0000279      0.001%          0.3516          1.50%
exp6         0.0000031      0.00011%        0.1172          0.50%
exp7         0.0000003      0.000011%       0.0355          0.15%

Кредиты
Эти реализации exp() были рассчитаны с помощью "scoofy" с использованием серии Тейлора из реализации tanh() "fuzzpilz" (кем бы они ни были, у меня были только эти ссылки на мой код).

Ответ 3

Приближения серии Тейлора (такие как функции expX() в Adriano ответить) наиболее точны около нуля и могут иметь огромные ошибки при -20 или даже -5. Если вход имеет известный диапазон, например от -20 до 0, как и исходный вопрос, вы можете использовать небольшую таблицу поиска и один дополнительный умножить, чтобы значительно повысить точность.

Фокус в том, чтобы распознать, что exp() можно разделить на целые и дробные части. Например:

exp(-2.345) = exp(-2.0) * exp(-0.345)

Дробная часть всегда будет находиться между -1 и 1, поэтому приближение рядов Тейлора будет довольно точным. Целочисленная часть имеет только 21 возможное значение для exp (-20) до exp (0), поэтому они могут храниться в небольшой таблице поиска.

Ответ 4

Я изучил статью Николь Швахольф, где теперь была более подробно описана первоначальная реализация C вышеупомянутой функции. Похоже, что, возможно, невозможно существенно утвердить точность вычисления exp, не оказывая серьезного влияния на производительность. С другой стороны, приближение справедливо и для больших величин х, до +/- 700, что, конечно, выгодно.

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

"exp" >= exp for all x                      1072693248 -  (-1) = 1072693249
"exp" <= exp for all x                                 - 90253 = 1072602995
"exp" symmetric around exp                             - 45799 = 1072647449
Mimimum possible mean deviation                        - 68243 = 1072625005
Minimum possible root-mean-square deviation            - 60801 = 1072632447

Он также указывает, что на "микроскопическом" уровне примерная "exp" функция демонстрирует поведение лестничной клетки, поскольку 32 бита отбрасываются при преобразовании из long в double. Это означает, что функция является кусочно-постоянной в очень малом масштабе, но функция, по крайней мере, никогда не уменьшается с ростом x.

Ответ 5

В следующем коде должны быть указаны требования к точности, так как для входов в [-87,88] результаты имеют относительную погрешность <= 1,73e-3. Я не знаю С#, так что это код C, но преобразование должно быть неудачным.

Я предполагаю, что, поскольку требование точности является низким, использование вычисления с одной точностью является прекрасным. Используется классический алгоритм, в котором вычисление exp() сопоставляется с вычислением exp2(). После преобразования аргумента посредством умножения на log2 (e) экспонента по дробной части обрабатывается с использованием минимаксного многочлена степени 2, тогда как экспоненциальная составляющая аргумента выполняется путем непосредственного манипулирования экспоненциальной частью одиночного элемента IEEE-754 -резервный номер.

Волатильный союз разрешает повторную интерпретацию битового шаблона как целого числа или числа с плавающей запятой с одной точностью, необходимого для манипуляции экспонентами. Похоже, что С# предлагает для этой функции функции повторной интерпретации, которая намного чище.

Две потенциальные проблемы производительности - это функция floor() и float- > int conversion. Традиционно оба были медленными на x86 из-за необходимости обрабатывать динамическое состояние процессора. Но SSE (в частности SSE 4.1) содержит инструкции, которые позволяют быстро выполнять эти операции. Я не знаю, может ли С# использовать эти инструкции.

 /* max. rel. error <= 1.73e-3 on [-87,88] */
 float fast_exp (float x)
 {
   volatile union {
     float f;
     unsigned int i;
   } cvt;

   /* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */
   float t = x * 1.442695041f;
   float fi = floorf (t);
   float f = t - fi;
   int i = (int)fi;
   cvt.f = (0.3371894346f * f + 0.657636276f) * f + 1.00172476f; /* compute 2^f */
   cvt.i += (i << 23);                                          /* scale by 2^i */
   return cvt.f;
 }