Использование серии Тейлора - не самый простой и не самый быстрый способ сделать это. Большинство профессиональных реализаций используют аппроксимирующие полиномы. Я покажу вам, как сгенерировать один в Maple (это программа компьютерной алгебры), используя алгоритм Remez .
Для 3 цифр точности выполните следующие команды в Maple:
with(numapprox):
Digits := 8
minimax(ln(x), x = 1 .. 2, 4, 1, 'maxerror')
maxerror
Его ответ следующий полином:
-1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x
С максимальной погрешностью: 0,000061011436
Мы сгенерировали полином, который приближает ln (x), но только внутри интервала [1..2]. Увеличение интервала не является разумным, потому что это увеличит максимальную ошибку еще больше. Вместо этого сделайте следующее разложение:
Итак, сначала найдите наибольшую степень 2, которая по-прежнему меньше, чем число (см .: Какой самый быстрый / самый эффективный способ найти старший установленный бит (мсб) в целом числе в C? ). Это число на самом деле является логарифмом основания-2. Разделите это значение, и результат попадет в интервал 1..2. В конце нам нужно будет добавить n * ln (2), чтобы получить окончательный результат.
Пример реализации для чисел> = 1:
float ln(float y) {
int log2;
float divisor, x, result;
log2 = msb((int)y); // See: https://stackoverflow.com/a/4970859/6630230
divisor = (float)(1 << log2);
x = y / divisor; // normalized value between [1.0, 2.0]
result = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
result += ((float)log2) * 0.69314718; // ln(2) = 0.69314718
return result;
}
Хотя, если вы планируете использовать его только в интервале [1.0, 2.0], функция выглядит так:
float ln(float x) {
return -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
}