Самый точный способ вычисления asinhf() из log1pf()?

Question

Обратная гиперболическая функция asinh() тесно связана с натуральным логарифмом. Я пытаюсь определить наиболее точный способ вычисления asinh() из стандартной математической функции C99 log1p(). Для удобства экспериментов я ограничиваю себя вычислением одноточечной обработки IEEE-754 прямо сейчас, то есть смотрю asinhf() и log1pf(). Я намерен повторно использовать тот же самый алгоритм для вычисления двойной точности, то есть asinh() и log1p(), позже.

Моя основная цель - минимизировать ошибку ulp, вторичная цель - свести к минимуму количество неверно округленных результатов, при условии, что улучшенный код будет минимально медленнее, чем версии, указанные ниже. Любое постепенное улучшение точности, скажем 0,2 ulp, было бы очень желанным. Добавление нескольких FMA (fused multiply-add) было бы неплохо, с другой стороны, я надеюсь, кто-то сможет определить решение, в котором используется быстрый rsqrtf() (обратный квадратный корень).

Полученный код C99 должен поддаваться векторизации, возможно, некоторыми незначительными прямыми преобразованиями. Все промежуточные вычисления должны соответствовать с точностью до аргумента функции и результата, так как любой переключатель с более высокой точностью может оказать серьезное отрицательное влияние на производительность. Код должен корректно работать как с поддержкой денормальной поддержки IEEE-754, так и в режиме FTZ (с нуля до нуля).

До сих пор я идентифицировал следующие две реализации кандидата. Обратите внимание, что код может быть легко преобразован в безветровую векторную версию с одним вызовом до log1pf(), но я не сделал этого на этом этапе, чтобы избежать ненужной запутывания.

/* for a >= 0, asinh(a) = log (a + sqrt (a*a+1)) 
         = log1p (a + (sqrt (a*a+1) - 1)) 
         = log1p (a + sqrt1pm1 (a*a)) 
         = log1p (a + (a*a/(1 + sqrt(a*a + 1)))) 
         = log1p (a + a * (a/(1 + sqrt(a*a + 1)))) 
         = log1p (fma (a/(1 + sqrt(a*a + 1)), a, a) 
         = log1p (fma (1/(1/a + sqrt(1/a*a + 1)), a, a) 
*/ 
float my_asinhf (float a) 
{ 
    float fa, t; 
    fa = fabsf (a); 
#if !USE_RECIPROCAL 
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation 
     t = log1pf (fa) + 0x1.62e430p-1f; // log(2) 
    } else { 
     t = fmaf (fa/(1.0f + sqrtf (fmaf (fa, fa, 1.0f))), fa, fa); 
     t = log1pf (t); 
    } 
#else // USE_RECIPROCAL 
    if (fa > 0x1.0p126f) { // prevent underflow in intermediate computation 
     t = log1pf (fa) + 0x1.62e430p-1f; // log(2) 
    } else { 
     t = 1.0f/fa; 
     t = fmaf (1.0f/(t + sqrtf (fmaf (t, t, 1.0f))), fa, fa); 
     t = log1pf (t); 
    } 
#endif // USE_RECIPROCAL 
    return copysignf (t, a); // restore sign 
} 

С конкретной log1pf() реализации, который с точностью до 0,6 < ulps, я наблюдаю следующую статистику ошибок при тестировании исчерпывающе во всех 2 возможных IEEE-754 с одинарной точностью входов. Когда USE_RECIPROCAL = 0, максимальная ошибка равна 1.49486 ulp, и есть 353 587 822 неверно округленные результаты. С USE_RECIPROCAL = 1 максимальная ошибка 1.50805 ulp, и только 77 569 390 неверно округленных результатов.

С точки зрения производительности, вариант USE_RECIPROCAL = 0 будет быстрее, если обратные и полные деления будут занимать примерно одинаковое количество времени, но вариант USE_RECIPROCAL = 1 может быть быстрее, если доступна очень быстрая взаимная поддержка. Ответы могут предполагать, что все основные арифметические операции, в том числе FMA (плавное умножение), правильно округлены в соответствии с режимом «круглый-самый-четный» IEEE-754. Кроме того,, более быстрая, почти правильно закругленная, версии взаимных и rsqrtf()могут быть доступны, где «почти правильно округленная» означает, что максимальная ошибка ulp будет ограничена примерно 0,53 ulps и подавляющим большинством результатов, скажем, 95%, правильно округлены. Базовая арифметика с направленными закруглениями может быть доступной без каких-либо дополнительных затрат.

Answer 1

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

Поскольку мой вопрос о минимизации ошибки для преобразования аргумента, которые понесены в дополнение к ошибке в log1pf() себя, самый простой подход использовать для экспериментов является использование правильно закругленную реализацию этой функции логарифма. Обратите внимание на то, что реализация с корректным округлением вряд ли будет существовать в контексте высокопроизводительной среды. Согласно работам Ж.-М. Muller et. . Аль, для получения точных результатов с одинарной точностью, x86 расширенная точностью вычисления должны быть достаточными, например

float accurate_log1pf (float a) 
{ 
    float res; 
    __asm fldln2; 
    __asm fld  dword ptr [a]; 
    __asm fyl2xp1; 
    __asm fst  dword ptr [res]; 
    __asm fcompp; 
    return res; 
} 

Реализация asinhf() с использованием первого варианта от моего вопроса, то выглядит следующим образом:

float my_asinhf (float a) 
{ 
    float fa, s, t; 
    fa = fabsf (a); 
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation 
     t = log1pf (fa) + 0x1.62e430p-1f; // log(2) 
    } else { 
     t = fmaf (fa/(1.0f + sqrtf (fmaf (fa, fa, 1.0f))), fa, fa); 
     t = accurate_log1pf (t); 
    } 
    return copysignf (t, a); // restore sign 
} 

Тестирование со всеми 2 Одиночные операции с одиночной точностью IEEE-754 показывают, что максимальная погрешность 1,49486070 ulp происходит при ± 0x1.ff5022p-9 и имеется 353 521 140 неправильных округленных результатов. Что произойдет, если преобразование всего аргумента использует арифметику с двойной точностью? Код изменен на

float my_asinhf (float a) 
{ 
    float fa, s, t; 
    fa = fabsf (a); 
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation 
     t = log1pf (fa) + 0x1.62e430p-1f; // log(2) 
    } else { 
     double tt = fa; 
     tt = fma (tt/(1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt); 
     t = (float)tt; 
     t = accurate_log1pf (t); 
    } 
    return copysignf (t, a); // restore sign 
} 

Однако с этой заменой ошибка не улучшилась! Максимальная погрешность 1,49486070 ульпа по-прежнему наблюдается при ± 0x1.ff5022p-9, и сейчас имеется 350 971 046 неправильно округленных результатов, немного меньше, чем раньше. Проблема заключается в том, что операнд float не может передавать достаточную информацию до log1pf() для получения более точных результатов. Аналогичная проблема возникает при вычислении sinf() и cosf(). Если приведенный аргумент, представленный как правильно округленный операнд float, передается на основные полиномы, результирующая ошибка в sinf() и cosf() составляет всего чуть меньше 1,5 ulp, как мы здесь наблюдаем с my_asinhf().

Одним из решений является вычисление преобразованного аргумента более высокой, чем одинарная точность, например, в виде пары операндов с двойным поплавком (полезный краткий обзор методов с двойным поплавком можно найти в this paper by Andrew Thall). В этом случае мы можем использовать дополнительную информацию для выполнения линейной интерполяции по результату на основе знания о том, что производная от логарифма является обратной. Это дает нам:

float my_asinhf (float a) 
{ 
    float fa, s, t; 
    fa = fabsf (a); 
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation 
     t = log1pf (fa) + 0x1.62e430p-1f; // log(2) 
    } else { 
     double tt = fa; 
     tt = fma (tt/(1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt); 
     t = (float)tt;    // "head" of double-float 
     s = (float)(tt - (double)t); // "tail" of double-float 
     t = fmaf (s, 1.0f/(1.0f + t), accurate_log1pf (t)); // interpolate 
    } 
    return copysignf (t, a); // restore sign 
} 

Исчерпывающий тест этой версии указывает на то, что максимальная ошибка была уменьшена до 0.99999948 ULP, это происходит при ± 0x1.deeea0p-22. Есть 349 653 534 неверно округленных результатов. Достигнута добросовестная реализация asinhf().

К сожалению, практическая полезность этого результата ограничена. В зависимости от платформы HW пропускная способность арифметических операций на double может составлять от 1/2 до 1/32 от пропускной способности float. Вычисление с двойной точностью может быть заменено вычислением с двойным поплавком, но это также потребует значительных затрат. Наконец, мой подход здесь заключался в том, чтобы использовать реализацию с одной точностью как доказательство для последующей работы с двойной точностью, и многие аппаратные платформы (конечно, все, что меня интересует) не предлагают аппаратную поддержку для числового формата с более высокой точностью чем IEEE-75 (двойная точность). Поэтому любое решение не должно требовать арифметики более высокой точности при промежуточном вычислении.

Поскольку все неприятные аргументы в случае asinhf() малы по величине, можно было [частично?] Задать вопрос о точности, используя полиномиальное минимаксное приближение для области вокруг начала координат. Поскольку это создало бы другую ветвь кода, это, скорее всего, сделало бы векторизацию более сложной.

Answer 2

Во-первых, вы можете взглянуть на точность и скорость вашей функции log1pf: они могут немного различаться между libms (я обнаружил, что математические функции OS X бывают быстрыми, glibc - медленнее, но обычно правильно округлены).

Openlibm, на основе BSD libm, которая, в свою очередь, базируется на fdlibm Sun, использовать несколько подходов по диапазону, но основной бит the relation:

t = x*x; 
w = log1pf(fabsf(x)+t/(one+sqrtf(one+t))); 

Вы также можете попытаться скомпилировать с -fno-math-errno, который отключает старые коды ошибок System V для sqrt (исключения IEEE-754 будут по-прежнему работать).