Nowa metoda faktoryzacji 'trial division'

Cechą podzielności n przez liczbę a jest możliwość zapisania n = a*p+a.
Zastanawiało mnie, czy można łatwo przechodzić z liczby postaci a*p+b do (a+1)*q+c, gdzie q jest bliskie p zaś b<2a, c<2(a+1).
Można, co doprowadziło do uzyskania kolejnego siłowego algorytmu faktoryzacji typu 'trial division'.
W postaci n=a*p+b największą trudność wydaje się zachowywać b w przedziale (0,2a), lecz gorsze jest znajdowanie kolejnej podstawy systemu q. Znajdujemy ją za pomocą dzielenia n/(a+1), co nie zawsze trafia we właściwą wartość.

Oto jedna z wersji tego algorytmu rozkładu liczby n:
a = b = 1; p=n-1;
while( a<p ) {
  a++;
  b = b-p%a;
  p = p-p/a;
  if( b<a ) { p--; b = b+a; }
  if( a==b ) return dzielniki a, p+1
}
return liczba pierwsza;

Sposób ten niezwykle szybko pozwala znajdować kolejne reszty b, ale tylko dla sąsiedniej wartości a.
Oto początek rozkładu liczby 8934053:
a  b  p
1  1  8934052 
2  3  4467025  (4467025:3=1489008, 4467025-1489008-1 = 2978016)
3  5  2978016  (2978016:4=744504, 2978016-744504 = 2233512)
4  5  2233512  (2233512:5=446702, 2233512-446702-1 = 1786809)
5  8  1786809
6 11 1489007
7  9  1276292
...

Wartości p początkowo szybko maleją dążąc do pierwiastka kwadratowego z n, potem różnica między p na wejściu oraz wyjściu z pętli zmniejsza się do jedynki. Nie jest znana dokładnie, lecz błąd jest rzędu a, czyli stosunkowo mały. Dzielenia są prostsze niż w standardowym algorytmie dzielenia przez kolejne liczby.

Jeśli chcemy mieć mniejsze reszty b, należy zapamiętać wartość c=b z początku pętli oraz zmienić warunek pod if na tę kopię if(c<a)... Wtedy zamiast
5 8 1786809
mamy
5 3 1786810
Przebieg ten sam. Algorytm rekursywny (nie równoważny powyższemu, robiący to samo)
factd(a,p,b) {
  if( a==b ) return dzielniki a, p+1;
  if( a>p ) return liczba pierwsza;
  q=p-p/(a+1);
  c=p%b;
  if( b<a+1 ) {q--; c=c+a+1; }
  return factd(a+1, q, c);
}

Liczby złożone nieparzyste charakteryzują się nieparzystymi a i b=a oraz parzystymi p. Pojawia się chęć przeskoków tylko po nieparzystych a, jak w algorytmie standardowym. Jest to możliwe, lecz wtedy należy użyć konwersji. Tracimy łatwość wyznaczania reszt b. Pojawia się dodatkowe mnożenie, oraz konieczność korekty, gdyż b ma często zbyt dużą wartość. Wyliczane kolejne ilorazy nie są równe. Tworzą zaszumiony ciąg statystycznie nierosnący (fragmenty ciągu kolejnych ilorazów mogą być postaci: ..., 273, 270, 268, 265,..., 2, 1, 1, 1, 2, 1, 1, ).