Факторизация рекуррентным методом

[Dmitriy40]

Учитывая что по таблице выше метод Ферма выигрывает всегда когда отношение сомножителей менее примерно 10-20, то вполне можно дойти именно им до примерно такого соотношения, а вот потом продолжать уже вашим.

[Побережный Александр]

Dmitriy40 в сообщении #1401415 писал(а):

Очень уж подозрительно выглядит почти линейное возрастание $b_i$, напрашивается какой-то метод более быстрого изменения $b_i$

Меня этот вопрос заинтересовал тоже. Предлагаю усовершенствование метода рекуррентной факторизации

Начальные условия прежние.

Побережный Александр в сообщении #1401136 писал(а):

Обозначения.
$n=[\sqrt{N}]$, где []- целая часть числа N.
$R=(n+1)^2-N$
$R'=N-n^2$, причем $R+R'=2n+1$.

$a_0=[\sqrt{R}]$, $b_0=[\frac{R'+a_0^2}{2(n-a_0)}]+1$

Начало цикла.
$N_i=(n-a_i)(n+a_i+2b_i)$, (причем должно выполняться $N_i>N$)
$s_i=n-2a_i-b_i$, ($s_i$ может быть и отрицательным)

$k_i=\lceil\frac{\sqrt{(s_i)^2+(N_i-N)}+s_i}3\rceil$, (где $\lceil \rceil$ округление в большую сторону)

$a_{i+1}=\lceil\sqrt{(n+b_i+k_i)^2-N}-(b_i+k_i)\rceil$

$b_{i+1}=\lceil\frac{R^\prime+a_{i+1}^2}{2(n-a_{i+1})}\rceil$

Цикл заканчивается, когда числа $a$ и $b$ станут натуральными

Dmitriy40 в сообщении #1401415 писал(а):

Например, даже для $N=101\cdot907$ надо более $120$ итераций.

Рассмотрим этот пример.

$N=101\cdot907=91607$

$n=302$

$R=202$

$R^\prime=403$

Следовательно, $a_0=14, b_0=2$

Проведем первый цикл.

$N_1=(302-14)(302+14+2\cdot2)=288\cdot320=92160$

$s_1=n-2\cdot14-2=272$
$ k_1=\lceil\frac{\sqrt{(272)^2+(92160-91607)}+272}3\rceil=182$

$a_1=\lceil\sqrt{(302+2+182)^2-91607}-(2+182)\rceil=197$

$b_1=\lceil\frac{403+197^2}{2(302-197)}\rceil=187$

Первый цикл закончен.
Алгоритм позволил перешагнуть сразу на $185$ единиц для числа $b$.

[Побережный Александр]

Уважаемые софорумники!
Прошу оценить эффективность модифицированного алгоритма факторизации рекуррентным методом (предыдущее сообщение).
К сожалению, у меня нет возможности работать с большими числами. Надеюсь усовершенствование ускоряет работу алгоритма.
Для меня Ваше мнение очень важно!

[mathpath]

В вашем алгоритме в теле цикла в двух местах извлекается квадратный корень
Какой формат должен быть у извлекаемого корня - целочисленный или с плавающей точкой ?

[Dmitriy40]

mathpath в сообщении #1410848 писал(а):

В вашем алгоритме в теле цикла в двух местах извлекается квадратный корень
Какой формат должен быть у извлекаемого корня - целочисленный или с плавающей точкой ?

Очевидно с плавающей точкой, иначе (целочисленный) как минимум не было бы смысла округлять (причём вверх!) значения $a_i$, они всегда автоматом были бы целыми.

[mathpath]

Насколько мне неизвестно, для сильно больших чисел не существует эффективных алгоритмов извлечения корня
с плавающей точкой. Эффективные - значит очень быстрые.

[Dmitriy40]

Ну здесь можно обойтись и целочисленным корнем, если проверять точно ли он извлёкся, а в $k_i$ вместо округления корня проверять подкоренное выражение в зависимости от остатка $s_i \bmod 3$. Или использовать целочисленный корень, а в конце проверить подстановкой нужно ли увеличивать результат на 1. Но это всё усложнение алгоритма.
И кстати тот же метод Ньютона при подходящем выборе начального приближения даёт квадратичную сходимость, т.е. не так уж и много придётся сделать длинных медленных делений.
В общем это решаемо.

[Побережный Александр]

Уважаемый Dmitriy40,
если вы работали с модифицированным алгоритмом, можете ли вы описать, как он работает с составными числами со 100 десятичными знаками?
Какое максимальное время затрачивается на "неудобные" числа?
И вообще, может возникают новые проблемы при работе алгоритма!

[Dmitriy40]

Что-то он (от 10 июля выше) не на всех числах работает (первые 5 из 85 итераций):

N=101*577:
i=0, a=16, b=2, n=241, N=58277, R=287, R'=196
i=1, a=156, b=145, Ni=58725, s=207, k=139
i=2, a=157, b=148, Ni=58395, s=18446744073709551400, k=1
i=3, a=158, b=152, Ni=58296, s=18446744073709551395, k=1
i=4, a=159, b=156, Ni=58349, s=18446744073709551389, k=1
i=5, a=160, b=160, Ni=58384, s=18446744073709551383, k=1

(Код PARI/GP)

n=sqrtint(N);
R=(n+1)^2-N; Rs=N-n^2;
a=sqrtint(R); b=floor((Rs+a^2)/2/(n-a))+1; i=0;
print("i=",i,", a=",a,", b=",b,", n=",n,", N=",N,", R=",R,", R'=",Rs);
while(N%(n-a)!=0,
   Ni=(n-a)*(n+a+2*b);
   s=n-2*a-b;
   k=ceil((sqrt(s^2+(Ni-N))+s)/3);
   a=ceil(sqrt((n+b+k)^2-N)-(b+k));
   b=ceil((Rs+a^2)/2/(n-a));
   i++;
   printf("i=%u, a=%u, b=%u, Ni=%u, s=%u, k=%u\n",i,a,b,Ni,s,k);
);
print(n-a,"*",N\(n-a));