Недавно я написал код, который пытался вычислить без использования библиотечных функций. Вчера я просматривал старый код и пытался сделать его как можно быстрее (и исправить). Вот моя попытка:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Здесь я пытаюсь найти так , что е чуть больше п, а затем добавить значение логарифм п
Недавно я заинтересовался численным анализом, и поэтому я не могу не задать вопрос, насколько быстрее этот сегмент кода может быть выполнен на практике, будучи достаточно правильным? Нужно ли мне переключаться на некоторые другие методы, например, используя продолжение дроби, как это ?
ОБНОВЛЕНИЕ 1 : Используя ряд гиперболических арктанов, упомянутый в Википедии , вычисление, кажется, почти в 2,2 раза медленнее, чем функция журнала стандартной библиотеки С. Хотя я не проверил детально производительность, и для больших чисел моя текущая реализация кажется ДЕЙСТВИТЕЛЬНО медленной. Я хочу проверить мою реализацию на наличие ошибок и среднее время для широкого диапазона чисел, если я могу управлять. Вот мое второе усилие.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Любое предложение или критика приветствуется.
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
источник
frexpОтвет Кирилла уже затронул большое количество актуальных вопросов. Я хотел бы остановиться на некоторых из них, основываясь на практическом опыте проектирования математических библиотек. Заранее: разработчики математических библиотек, как правило, используют каждую опубликованную алгоритмическую оптимизацию, а также множество машинно-специфических оптимизаций, не все из которых будут опубликованы. Код часто будет написан на ассемблере, а не скомпилированным кодом. Поэтому маловероятно, что простая и скомпилированная реализация достигнет более 75% производительности существующей высококачественной математической реализации библиотеки, предполагая идентичные наборы функций (точность, обработка особых случаев, отчеты об ошибках, поддержка режима округления).
Точность обычно оценивается путем сравнения с (сторонним) эталоном с более высокой точностью. Функции с одинарной точностью и одним аргументом легко тестируются исчерпывающе, другие функции требуют тестирования с помощью (направленных) случайных тестовых векторов. Очевидно, что невозможно вычислить бесконечно точные эталонные результаты, но исследование дилеммы Table-Maker предполагает, что для многих простых функций достаточно вычислить эталон с точностью, примерно в три раза превышающей точность цели. Смотрите, например:
Винсент Лефевр, Жан-Мишель Мюллер, «Худшие случаи для правильного округления элементарных функций в двойной точности». В материалах 15-го симпозиума IEEE по компьютерной арифметике , 2001, 111-118). (препринт онлайн)
С точки зрения производительности следует различать оптимизацию по задержке (важно, когда рассматривается время выполнения зависимых операций) и оптимизацию по пропускной способности (актуально при рассмотрении времени выполнения независимых операций). В течение последних двадцати лет распространилось множество аппаратных методов распараллеливания, таких как параллелизм на уровне команд (например, суперскалярные, неупорядоченные процессоры), параллелизм на уровне данных (например, инструкции SIMD) и параллелизм на уровне потоков (например, гиперпоточность, многоядерные процессоры) привел к акценту на вычислительную пропускную способность как более актуальный показатель.
Объединенная операция множественного добавления ( FMA ), впервые представленная IBM 25 лет назад и доступная на всех основных процессорных архитектурах, является важнейшим строительным блоком современных реализаций математической библиотеки. Это обеспечивает уменьшение ошибок округления, дает ограниченную защиту от вычитания и значительно упрощает арифметику двойных двойных чисел .
C99log()C99fma()источник