пусть у нас есть числа a и b и мы хотим быстро посчитать a/b с большой точностью (а складывать-умножать числа мы уже умеем)
можно сдвинуть числитель и знаменатель на степень двойки, так что достаточно научиться находить 1/b для b между, скажем, 1/2 и 1
в этом посте нет экспериментальной математики, но так понравилась история, что не могу не поделиться — спасибо рассказавшему про это А.Гасникову (все ошибки, естественно, на моей совести)
1.
когда надо вычислить значение функции, у меня первый рефлекс — разложить ее в ряд Тейлора, т.е. в данном случае просто воспользоваться бесконечной геометрической прогрессией:
если b=1-q, то 1/b=1+q+q²+q³+… — сиди и вычисляй столько членов, сколько тебе нужно (так как |q|<1/2, рано или поздно всё получится)… но сколько нужно? чтобы найти N (двоичных) знаков после запятой нужно взять ~N членов, т.е. сделать ~N умножений, и это не очень вдохновляет
в конце концов, уравнение f(x)=0 для монотонной функции b всегда можно решать методом деления пополам (выбираем ту половину, на концах которой у f разные знаки, повторяем процесс) — уже это позволяет искать N знаков после запятой за ~N действий (для деления такой способ так же известен как деление в столбик)
но оказывается, что делить можно и намного быстрее, чем в столбик!
2.
если функция достаточно хорошая, то уравнение f(x)=0 можно быстро приближенно решать при помощи метода Ньютона
(
напомню идею: если x — приближенное значение корня, то рядом ним график функции недалеко ушел от касательной, поэтому в качестве следующего приближения можно взять пересечения касательной с нулем, т.е. x→x-f(x)/f'(x)
и в некоторой окрестности корня метод Ньютона, если производная в этом корне не равна 0, сходится очень быстро: за итерацию погрешность ~возводится в квадрат
)
на первый взгляд, нам метод Ньютона не поможет, так как в него входит деление — но тут происходит чудо: если сформулировать нашу задачу как задачу поиска нуля функции f(x)=1/x-b, то в методе Ньютона все деления сокращаются: f/f'=(1/x-b)⋅(-x²)=-x⋅(1-bx)
и получается рецепт x→x⋅(2-bx), который позволяет получить N знаков числа 1/b всего за ~log(N) операций (за каждую операцию количество верных знаков ~удваивается)
можно проверить, как это работает:
from mpmath import *
mp.dps = 300
b = mpf(57)/100
x = mpf(1)
print("1/"+nstr(b,30))
print(0,":",nstr(x,80))
for k in range(10):
x = x*(2-b*x)
print(k+1,":",nstr(x,80),
"diff:",nstr(abs(1/b-x),2,min_fixed=1))
print("T :",nstr(1/b,80))
3.
что это всё-таки за странная формула x→x(2-bx), можно ли это связать с чем-то более знакомым?
оказывается, это просто способ быстро вычислять частичные суммы всё той же геометрической прогрессии! действительно, если b = 1-q, и x = (1-q^N)/(1-q), то x’ = (1-q^N)/(1-q) ⋅ (1+q^N) = (1-q^{2N})/(1-q) — т.е. за шаг мы переходим от суммы первых N членов к сумме первых 2N членов геометрической прогрессии — немножко похоже на быстрое возведение в степень, только еще формулы сокращенного умножения знать надо
===
что можно ускорять какие-то базовые операции, меня впечатляет; вот небольшой текст про это в Мат. составляющей (в т.ч. про быстрое умножение): https://book.etudes.ru/articles/fast-arithmetic/
метод Ньютона здесь уже появлялся раньше, и будет, думаю, обсуждаться еще