Метод Ньютона та квадратні корені — квадратична збіжність і швидкий обернений корінь
Метод Ньютона — елегантний алгоритм чисельної математики: наближення дотичною уточнює оцінку кореня, подвоюючи точність з кожною ітерацією. Для квадратного кореня це стародавній вавилонський метод; для оберненого кореня з бітовою хитрістю — легендарний алгоритм освітлення Quake III Arena.
1. Виведення методу Ньютона-Рафсона
Щоб знайти корінь рівняння f(x) = 0: починаємо з припущення x₀, проводимо дотичну пряму до кривої в точці (x₀, f(x₀)) і знаходимо, де вона перетинає нуль.
2. Доведення квадратичної збіжності
3. Квадратний корінь методом Ньютона
4. Вавилонський метод
Ітерація x_{n+1} = (x + S/x)/2 для квадратних коренів була відома вавилонянам приблизно 1700 року до н. е. і записана на табличці YBC 7289. На табличці показано √2 ≈ 1.41421296..., обчислене з точністю до 6 шістдесяткових розрядів (за основою 60).
Герон Александрійський (бл. 50 р. н. е.) описав той самий алгоритм, і в західній традиції його іноді приписують йому. Метод також називають «методом Герона». Це найдавніший відомий приклад ітераційного чисельного алгоритму.
5. Швидкий обернений квадратний корінь із Quake III
1999 року Quake III Arena від id Software містила функцію
Q_rsqrt, яка обчислювала 1/√x із точністю ≈ 1.6% всього за кілька
операцій — без ділення й виклику sqrt. Вона була потрібна для обчислення
нормалей освітлення (нормалізація векторів вимагає ділення на їхню
довжину).
// Початковий код Quake III (C), приблизно 1999 рік
// Приписується Джону Кармаку; магічна стала, ймовірно, від Cleve Moler / Greg Walsh
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long *)&y; // бітова реінтерпретація float як цілого
i = 0x5f3759df - (i >> 1); // початкове припущення з магічної сталої
y = *(float *)&i; // бітова реінтерпретація цілого назад у float
y = y * (threehalfs - (x2 * y * y)); // ОДНА ітерація методу Ньютона-Рафсона
return y;
}
Ключова ідея: число з рухомою комою IEEE 754 подається в пам'яті як (1+m)×2^e. Бітове подання — це (e + 127) у старших бітах і m у молодших. Тож log₂(x) ≈ i/L − 127, де i — це цілочисловий бітовий шаблон, а L = 2²³. Зсув праворуч (i >> 1) приблизно ділить log₂(x) на 2 — що у звичайному просторі ділить порядок на 2, тобто бере квадратний корінь. Віднімання від магічної сталої перетворює результат на 1/√x. Одна ітерація методу Ньютона-Рафсона (для f(y) = 1/y² − x = 0) дає y = y(3/2 − x/2 · y²), уточнюючи до похибки ~1.6%.
0x5f3759df залишається одним із найвідоміших
чисел в історії розробки ігор.
6. Коли метод Ньютона не працює
- Кратні корені: поблизу подвійного кореня f(r) = f'(r) = 0 збіжність погіршується до лінійної (а не квадратичної). Використовуйте модифікований метод Ньютона: x_{n+1} = xₙ − m·f(xₙ)/f'(xₙ) із кратністю m, щоб відновити квадратичну збіжність.
- Осциляція: деякі функції породжують цикли x₀ → x₁ → x₀ → ... Періодичні цикли відображень Ньютона — це шлях до хаосу. Приклад: f(x) = x³ − 2x + 2 має 2-цикл для x₀ = 0.
- Розбіжність: якщо f'(xₙ) близьке до нуля, наступне наближення може опинитися дуже далеко. Завжди потрібне гарне початкове припущення; глобальна збіжність не гарантована.
- Фрактали Ньютона: для комплексних коренів многочленів басейни притягання для різних коренів утворюють гарні фрактальні межі — невеликі зміни в x₀ можуть призвести до збіжності до зовсім різних коренів.