Log2 (_mm_log2_pd intel - ламеры :))

Хотелось бы реализовать полноценный метод вычисления логарифма от double. Толи я не те запросы пишу, то ли поисковик обнищал. Первоначальная идея вычислит результат по битам как и в случае с целочисленным логарифмом:

Сомнительна возможность добиться скорости большей, чем у FPU. Синус/косинус вычисляется по таблицам медленнее, чем рядами, во всяком случае, а log - несколько более дешевая операция для сопроцессора, для него еще менее выгодно.

ну 60-50 тактов на FPU log, по спецификации интел, что собственно у меня и получается при подсчете тактов.
синус и косинус 120 тактов. В радах вторая приведенная мной формула за 10-20 тактов для float справляется. Синус с косинусом на SSE в виде ряда в разы быстрее 20-30.

foxes
> ну 60-50 тактов на FPU log, по спецификации интел, что собственно у меня и получается при подсчете тактов.
Советую все-таки перестать жить в прошлом веке и почитать про branch misprediction.
Вариант с кучей if-ов просто ужасен и более чем полностью состоит из ошибок предсказания ветвлений.
Если в случае целых чисел компилятор может попробовать реализовать все с помощью cmove, то в случае вещественных это уже на порядок сложнее.

Вариант с отбрасыванием мантиссы и коррекции полиномом — это наиболее очевидный и достаточно эффективный вариант.
Коэффициенты для полиномов могут быть получены разложением по полиномам Лежандра или просто методом наименьших квадратов.
Еще можно заморочиться и использовать минимаксный алгоритм.

}:+()___ [Smile]
> Вариант с отбрасыванием мантиссы и коррекции полиномом
Вот собственно по нему и вопрос.
}:+()___ [Smile]
> разложением по полиномам Лежандра
Забавно, но меня больше интересует метод Ньютона. Натуральный логарифм единственное что в нем не понятно чем по мнению автора является арифметико-геометрическое среднее.
}:+()___ [Smile]
> Советую все-таки перестать жить в прошлом веке и почитать про branch
> misprediction.
Если тебя смущает ход моих рассуждений, то живи лучше там сам. Меня прежде всего интересует точность вычислений.

Ghost2
http://www.chokkan.org/software/dist/fastexp.c.html

}:+()___ [Smile]
Тебе надо в торговлю идти работать, и обязательно продавать парики волосатым, а руки безногим, у тебя в этом талант.

foxes
>Собственно вопрос может кто в курсе как можно посчитать эти коэффициенты для полинома 6-й,7-й и тд степени.
В приведённом тобой посте вроде открытым текстом написано:

И даются ссылки таки на алгоритм Ремеза:
http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/ht… /minimax.html
http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/ht… rs/remez.html
на который уже раньше давалась ссылка в этой теме:
https://ru.wikipedia.org/wiki/Алгоритм_Ремеза

Написал следующий (убогий) код на Maxima:

poly(n):=sum(a[k]*x^k,k,0,n);
L2_squared(n):=integrate((poly(n)-log(x)/log(2))^2,x,1,2);
sol(n):=solve(makelist(diff(L2_squared(n),a[k]),k,0,n));
p(n):=ev(poly(n),sol(n),eval);

grind(ev(horner(p(5),x),simplify,numer));
sqrt(ev(L2_squared(5),sol(5),eval,simplify,numer));

реализующий упомянутый в #3 метод наименьших квадратов.
Стоит отметить, что он минимизирует L², а не L^∞, как алгоритм Ремеза.
Результаты:

>Забавно, но меня больше интересует метод Ньютона. Натуральный логарифм единственное что в нем не понятно чем по мнению автора является арифметико-геометрическое среднее.
http://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean
и
http://en.wikipedia.org/wiki/Logarithm

>Меня прежде всего интересует точность вычислений.

А вот здесь можно подробнее? Чего хочется? Если цель - добиться 1 ULP дабла... да ещё и быстрее процессора... благородная цель.

Что интелы говорят:

http://www.intel.com/content/www/us/en/processors/architectures-s… -manuals.html

и это за ~100 тактов.

они правда умеют и ошибаться.

Примечание: в этом посте я могу сплошняком ошибаться, ибо вычислительную математику помню поверхностно.

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

double _STDCALL _DLog2DEF(double _a)
{
double b=4/(_a*18446744073709551616.0);
  double a=1.0;
  while ((a-b)>0.0) {
    double fa=(a+b)*0.5;
    double fb=_DsqrtDEF(a*b);
    a=fa;
    b=fb;
  }

  return (4.5323601418271938096276829457166/(2.0*a))-64;
}

Ни чем не лучше моего первого примера. :)

Что же касается полинома, то у интелов в общем то все уже посчитано и реализовано в _mm_log2_ps и _mm_log2_pd - здесь мантисса результата считается полностью, для тех значений которые использовал я особо ошибок не было.
Ту ссылку что я приводил выше, в принципе это некая реализация _mm_log2_ps https://searchcode.com/codesearch/view/20747275/
Но это сурагат, потому что у интелов совсем другой результат.

FordPerfect
> И дальше точность ухудшается. Не уверен в причинах этого явления.
Ну тут надо смотреть на корень зла, хотя бы если рассмотреть ряд Тейлора(Меркатора) для LN, собственно понятно что поделив каждый коэффициент ряда на ln(2) получим тот же log2.

>Ограничение этого бесконечного ряда i-м членом порождает многочлены Тейлора i-го порядка, содержащие степени не выше i-й. На рисунке >справа приведены графики функции \ln (1+x)\, и некоторых многочленов Тейлора около x = 0. Аппроксимации сходятся к функции только в >области сходимости −1 < x ≤ 1, а за её пределами быстро отклоняются от точной функции, причем многочлены высших степеней дают >бо́льшую ошибку.
причем многочлены высших степеней дают бо́льшую ошибку.
В чем тут связь между полиномом и рядом, да по сути это те же яйца, с теми же свойствами, только с боку. Можно попробовать по выносить за скобки некоторые параметры, поиграться с конечным числом элементов ряда, и в конечном итоге получить аналог того же полинома.
Точность действительно ухудшается, насколько я успел покопаться вчера с полиномами более низкого порядка то, на некоторых диапазонах, достаточно n=3 n=4, дает большую точность, но только для быстрого вычисления float.

Кстати у интелов там на этот счет стоит вилка с условиями, и они не особо заморачивались с производительностью, а вставили определение процессора и всех его свойств (cpuid) прямо в реализацию _mm_log2_pd, так что там прям на AVX/SSE2/SSE4 все считается.

FordPerfect
> И даются ссылки таки на алгоритм Ремеза:
Вот моя вина ссылки не разглядел.

FordPerfect
> А вот здесь можно подробнее? Чего хочется? Если цель - добиться 1 ULP дабла...
> да ещё и быстрее процессора... благородная цель.
Да это основная утопия, поэтому начал с очевидных алгоритмов дающих правильный результат.

На Ньютона я засмотрелся потому что как то писал квадратный корень, и скрестил итерационный алгоритм Ньютона (уточняется методом касательных) с Героном и вычислением обратного квадратного корня. Получилась очень симпатичная весьма точная и быстрая реализация:

double _DsqrtDEF(double _a)
{
  if (_a==0) return 0;
  conv4b c;
  conv4b al;
  c.valfloat=_a;
  al.valuint=0x5f3504f3-(c.valuint>>1);
  float g=al.valfloat;
  float mr=g*(3.0f-c.valfloat*g*g);
  double r=2.0/mr+_a*mr*0.5;
  return r*0.25+_a/r;
}

Поэтому экспериментирую в этом направлении.

каждое деление - 40 тиков. Корень квадратный не быстрее вычисляется, чем одно деление?

Zab
> каждое деление - 40 тиков.
Смотря что и на чем считать. Я бы так сказал, весь алгоритм приведенный выше считается за 40 тактов для double. Для float одно деление, и полностью на алгоритм тратиться 28 тактов. Но это полностью с вызовом функции и передачей параметров, а если считать только инструкции внутри метода то и того меньше. Конечно это не идет в сравнение с SSE инструкциями там на квадратный корень уходит 10 тактов для float и 18 тактов для double. FPU считает квадратный корень за 20 тактов.
Но результат гораздо быстрее итерационного метода Герона (r*0.25+_a/r). Конкретно эта реализация подходит для процессоров без команд вычисления квадратного корня. А иногда позволяет распараллелить вычисление квадратного корня с другими вычислениями. Потому что чистый sqrt блокирует конвейеры, и если считать честные нормали на AVX (не для графики) то получается очень медленно.

В сравнении даже могу привести пример расчета обычной нормали на SSE и с быстрым обратным квадратным корнем.
doulbe fast 14-20
float fast 10-14
double 40 - при том что 3-4 элемента вектора делятся, а не умножаются на обратное значение корня.
float 10-24
То есть тут конкретно особенность вычислений с двойной точностью на SSE, получается гораздо дешевле вычислить 3 деления чем один раз обратить значение корня и сделать еще 3 умножения.
На FPU гораздо проще посчитать обратное значение и сделать 3 умножения, что дает 70-45 тактов на нормаль.

foxes
>В принципе с Ньютоном я вчера успел разобраться, за счет поиска среднего между средне арифметическим и средне геометрическим алгоритм получается достаточно весомым, но точност можно получить любую
>Получилась очень симпатичная весьма точная и быстрая реализация

А чиселки можешь привести? Ошибка в ulp и время в тактах. На твоей машине. Можно на глаз.

Оставлю тут пару версий.

Аппроксимация многочленом, порядок 7, точность ~8.65e-7, лучшее, что у меня получается полиномом:

double log2(double x)
{
    return (((((((
        +0.01443375071119004)*x
        -0.17637586388975270)*x
        +0.94298726381922360)*x
        -2.89539350282091500)*x
        +5.64554947963608300)*x
        -7.37179520425447900)*x
        +7.07342938148042600)*x
        -3.23283443968268600;
}

Аппроксимация Паде, порядок 4,4, точность ~9.2e-8, при повышении порядка точность падает:

double p0(double x)
{
    return ((((
        +.1273677649270379)*x
        +0.921530451918903)*x
        -7.222874751510376)*x
        +1.641481529588268)*x
        +4.53249439799058;
}

double p1(double x)
{
    return ((((
        -.01308052782623356)*x
        -.4778151982277505)*x
        +1.252106129978092)*x
        +4.844649683354856)*x
        +1.0;
}

double log2(double x)
{
    return -p0(x)/p1(x);
}

Код на Maxima для их нахождения (метод наименьших квадратов):

r0(n):=sum(a[k]*x^k,k,0,n);
r1(n):=1+sum(a[k+n]*x^k,k,1,n);
L2_squared(n):=integrate((r0(n)+r1(n)*log(x)/log(2))^2,x,1,2);
sol(n):=solve(makelist(diff(L2_squared(n),a[k]),k,0,2*n));
p0(n):=ev(r0(n),sol(n),eval);
p1(n):=ev(r1(n),sol(n),eval);

grind(ev(p0(4),numer));
grind(ev(p1(4),numer));

Обе функции приведены для x∈[1;2). Экспонента отбирается отдельно.

Более общие аппроксимации Паде-Эрмита пока не пробовал.

Битовая функция, использующая соотношение log(x^2)=2*log(x):

const int N=50;
const double cf=1.0/(double(1ULL<<(uint64_t(N))));

inline double log2(double x)
{
    uint64_t n,m;
    double y;
    memcpy(&n,&x,8);
    m=((n&0x7FF0000000000000)>>52ULL)-1023;
    n=(n&0x000FFFFFFFFFFFFF);
    n=(n|0x3FF0000000000000);
    memcpy(&y,&n,8);
    for(int i=0;i<N;++i)
    {
        memcpy(&n,&y,8);
        m=2ULL*m+(n>>62ULL);
        n=(n&0x000FFFFFFFFFFFFF);
        n=(n|0x3FF0000000000000);
        memcpy(&y,&n,8);
        y=y*y;
    }
    return 2.0*double(m)*cf;
}

Точность ~2.0e-15, и настраивается. Но медленная.
Интересно сколько она работает на какой-нибудь другой машине? Можете отписаться, кто читает?

2.0e-15 - это 9 ulp (если не ошибаюсь).
N=20 даёт 2.0e-6.
N=53 даёт 4.0е-16, т. е. 2 ulp.
N=60 даёт 2.5e-16, что несколько странно. Видимо вычисления проводятся на FPU в double extended-precision floating-point.
Соответственно в SIMD может быть не так хорошо.

>На Ньютона я засмотрелся потому что как то писал квадратный корень, и скрестил итерационный алгоритм Ньютона (уточняется методом касательных) с Героном и вычислением обратного квадратного корня.
Знакомо выглядит.
http://www.beyond3d.com/content/articles/8/
http://www.beyond3d.com/content/articles/15/

>Конечно это не идет в сравнение с SSE инструкциями там на квадратный корень уходит 10 тактов для float и 18 тактов для double.
Есть ещё sqrtss vs. rsqrtss и способы догнать точность rsqrtss с одной итерацией метода Ньютона.

P. S. А у тебя C или C++? В C++ использование union для доступа к битовому представлению вроде бы undefined behavior.
P. P. S. На всякий случай упомяну статьи о получении явных формул для аппроксимаций Паде-Эрмита для логарифма:

J.A.C. Weideman, Padé approximations to the logarithm I: Derivation via differential equations
Wenchang Chu, Harmonic number identitiesand Hermite–Padé approximations to the logarithm function

Закодил такое:

double f_exp(double x)
{
    double y,z;
    x*=0.03125;
    z=1.0+x;
    y=x*x;
    z+=0.5*y;
    y*=x;
    z+=0.1666666666666666574*y;
    y*=x;
    z+=0.04166666666666666435*y;
    y*=x;
    z+=0.008333333333333333218*y;
    y*=x;
    z+=0.001388888888888888942*y;
    y*=x;
    z+=0.0001984126984126984125*y;
//    y*=x;
//    z+=2.480158730158730157e-005*y;
//    y*=x;
//    z+=2.755731922398589251e-006*y;
//    y*=x;
//    z+=2.755731922398589357e-007*y;
    z*=z;
    z*=z;
    z*=z;
    z*=z;
    z*=z;
    return z;
}

double f_log(double x)
{
    double y,z,w;
    y=(((((
        +0.1433639818182963)*x
        -0.9991961322375756)*x
        +2.911969145635241)*x
        -4.673045914613294)*x
        +4.889507570787543)*x
        -2.272607233668817;
    for(int i=0;i<3;++i)
    {
        z=f_exp(y);
        w=y*y;
        y-=(z-x)*(1.0-y+
            w*(
                0.5
               -0.16666666666666666*y
               +0.04166666666666666435*w
               -0.008333333333333333218*w*y));
    }
    return y;
}

Обе функции вроде бы имеют погрешность 1e-16 на интервалах [-1/2,+1/2] и [e^(-1/2),e^(+1/2)] соответственно.
Стоимость:
f_exp - 7 сложений, 19 умножений,
f_log - 48 сложений, 83 умножения.

В тактах не замерял, неизвестно, насколько конкретно мои замеры - релевантны.

Точность - для FPU, на SIMD может быть хуже.

Соображения?

ребят, где у вас логарифм в time-critical месте мог всплыть? ._.

FordPerfect
> А чиселки можешь привести?
Дык, уже привел, для квадратного корня:
foxes
> Я бы так сказал, весь алгоритм приведенный выше считается за 40 тактов для
> double.
А для логарифма через Ньютона смысла нет, очень медленно (секунд 10-20), нужно че то допиливать.

FordPerfect
> Ошибка в ulp и время в тактах. На твоей машине. Можно на глаз.
а ошибки сам глянуть можешь код выложен. Я считаю сколько бит последних не совпадает с правильным ответом, стараюсь держать до 3 бит разницы или 0. Логарифм по Ньютону не калибровал.

FordPerfect
> Знакомо выглядит.
Ну конечно же стандартный быстрый обратный квадратный корень - считается экспонента, а потом доводится методом касательной, и я довернул точность Героном, а заодно перевернул дробь. Героном гораздо проще получить большую точность там сразу 20 бит точности прибавляется. Для мантиссы в 52 бита:2 бита на стадии вычисления экспоненты (разница логарифмов), 10 бит считается ньютоном, 20 и 20 Героном - Получается полный double sqrt.

Suslik
не переживай это так, для тренировки мозга. А вообще ты знаешь стандартный FPU FYL2X за 60 тактов меня не устраивает.
Второй пример в нулевом посте 26 тактов и разница на 3 бита (для float), вот это интересно.

FordPerfect
> Соображения?
я тут интелевцев _mm_log2_pd реинженерил, у них большая таблица логарифмов и они промежуточные результаты доводят чем то напоминающим полином. Чтоб такой гадостью воспользоваться без заморочек нужно intel compiller иметь.
печально как то.

FordPerfect
> Битовая функция, использующая соотношение log(x^2)=2*log(x):
Жесть 2400 тактов. Где то ты первый бит не там поставил:
для X=14.5678
FYL2X = 3.8647111158218332
у тебя = 6.8647111158218319

FordPerfect
> N=60 даёт 2.5e-16, что несколько странно. Видимо вычисления проводятся на FPU в
> double extended-precision floating-point.
> Соответственно в SIMD может быть не так хорошо.
Ну ты больше 52 не посчитаешь, дальше просто не влезет, поэтому за счет округления на fpu с 80 бит в 64 бит и получается 2.5e-16
N=52 самое оптимальное.
По сути альтернатива моему первому примеру (я там в обратном направлении пошел - квадратный корень вычисляю). Только он за 70 тактов считается на x32 и за 40 на x64 - в противовес твоим 2000-3000 тактов при N=52. Но это надо сравнивать при N=18 600-700 тактов.
FordPerfect
> Соответственно в SIMD может быть не так хорошо.
Наоборот будет нормально. Тут надо не забывать про особенность чисел с плавающей точкой, можно считать отдельные части ответа во FLOAT главное их сумму получить в DOUBLE.

FordPerfect
> P. S. А у тебя C или C++?
С++