26 kwietnia 2013

Faktoryzacja przez dzielenie przez sume cyfr, implementacja

Faktoryzację polegającą na sprawdzaniu reszt z dzielenia przez kolejne liczby można zaimplementować mniej standardowo, sprawdzanych przypadków jest więcej, lecz złożoność arytmetyczna się polepsza.
Teoria jest taka, jak opisana w sąsiednim poście, kiedy szukałem ciągu konwersji wskazującego dzielnik.
Liczbę n przedstawiam w systemie o podstawie 2^i, oraz mając te m cyfr sprawdzam cechami podzielności przez liczby 2^i-k, k nieparzyste malejące od 2^i-1 do 1. Jeśli uzyskana suma cyfr dzieli / jest wielokrotnością podstawy, mamy dzielnik.

Algorytm wykorzystuje konwersję podwajającą podstawę. Oto jej implementacja.
Liczba dana jest tablicą unsigned int a[8], w której cyfra najmniej znacząca jest ostatnim elementem tablicy. Podstawą systemu jest unsigned int p.
[\code=cpp]
void divp( unsigned int * a, unsigned int p ) {
  unsigned int k=0x80; // 2^i
  unsigned int i=0; // iterator tablicy cyfr a
  unsigned int b=0; // zmienna pomocnicza
  while(i++<length) { // tu length = 8, tablica 8
    b = b*p + a[i];
    a[i] = b/k; // cyfra to czesc calkowita
    b %= k; // reszta przenoszona dalej
    k /= 2; // nastepny bit
  }
}
[\code]

Sam rozkład zaczyna się przedstawienia liczby w systemie szesnastkowej. Następuje sprawdzenie cech podzielności przez 2, 3, 5 (suma cyfr jest podzielna przez te wartości). Dalej podział na przedziały binarne, zaś w każdym z nich dla kolejnych wartości nieparzystych wykorzystywana cecha podzielności przez p.
[\code=cpp]
int rozklad() {
  unsigned int a[] = ; // inicjowanie
  unsigned int s = ; // suma dla inicjacji
  if(s) ; // sprawdzenie malych dzielnikow
  unsigned int p=9; // zaczynamy od systemu siodemkowego 16-9
  unsigned int k=7; // 16-9 = 7, do cechy podzielnosci
  while(1) {
    unsigned int i= a.length()-1; // na cyfrę najmniej znaczącą
    s = a[i]; // sumowanie, inicjacja
    unsigned int t=1; // zmienna wspomagająca znajdowanie wspolczynnikow
    while(--i) {
      t = (k*t)%p; // dla cyfry lub aktualnej lub wspolczynnika przy i-1
      s = s + (a[i]*t)%p; // aktualna cyfra
    }
    if( p==s ) return p; // dzielnikiem jest p
    if( 2>k ) {
      k=p; // przeliczenie postaci liczby a[]
      divp(a, k+1); // zmienia dlugosc a
    } else k -= 2; // do następnego dzielenia
    p += 2; // następna podstawa
    if() ;// wyjscie w przypadku liczby pierwszej
  }
}
[\code]

Złożoność modułu divp jest liniowa względem długości cyfr liczby n, czyli logarytmiczna. Pętla faktoryzacji rozdziela przypadki na przedziały [2^i, 2^(i+1)), których jest logarytm z n.
W każdym takim przedziale następuje sumowanie cyfr liczby modulo p, znowu mające krotność logarytmu z n. Pętla ma dokładnie 2^(i-1) przejść, zakończonych uruchomieniem funkcji divp() (liniowe, gdyż 2^(2i)<n<2^(2i+2)).
Najwięcej kosztują przypadki dla dużych k (ale wtedy a jest stosunkowo 'krótka'). Ogólnie złożoność pesymistyczna jest O(n ln n).
Z kolei iloczyny cząstkowe wykonywane modulo nie mają wartości większych niż potęga 2 przekraczająca n.
Jeden z czynników występujących iloczynów ma wykres piłokształtny (liniowo malejący o 2, gwałtowny wzrost do 2^i po osiągnięciu takiej podstawy).

23 kwietnia 2013

Heurystyka pierwszości zamiast rozkładu

Testowanie rozkładu liczby na czynniki ma bardzo dużo przypadków do sprawdzenia. I cóż z tego, że procesory testują miliony przypadków na sekundę. Obliczenia trwają godziny, dni, miesiące.
Zatem należy poszukać takich zależności, które zmniejszą tę krotność.

Jednym z pomysłów było zastosowanie cech podzielności. Znamy podzielność przez 9 (suma cyfr cyfry dziesiętnej dzieli się przez 9), przez 11 (naprzemienna suma cyfr liczby dziesiętnej dzieli się przez 11).
Dla systemów o podstawie p parzystej, ich analogi brzmią następująco:
liczba n przedstawiona w systemie p dzieli się przez p-1, gdy suma jej cyfr dzieli się przez p-1;
liczba n przedstawiona w systemie p dzieli się przez p+1, gdy naprzemienna suma jej cyfr dzieli się przez p+1;

Przykłady liczbowe.
7 ' 595 ' 485 _ p=1088
w systemie o podstawie p ma sumę 'cyfr' 1087, naprzemienną sumę -103. Zatem istnieje dzielnik i jest równy 1088-1 = 1087.
1 ' 3904 ' 3903 _p=10280
ma sumę 'cyfr' 7808, naprzemienną sumę 0, zatem dzielnikiem jest 10281.

Aby znajdować dzielniki, wystarczy rozpatrywać podstawy parzyste(!) 4*m, co już zmniejszy krotność przypadków o połowę.

Uda się wyciągnąć coś jeszcze. Przechodząc po małych podstawach: 4, 8, 16, 20, 24, 30, 36 oznaczmy przez a sumę, zaś przez b sumę naprzemienną. Dla tak małych podstaw dla liczb złożonych okazuje się, że największy wspólny dzielnik d = nwd(p,a) jest równy mniejszej z p, a.
Dla sumy naprzemiennej wystarczy, że ta suma będzie 0, 1 lub -1, aby istniały dzielniki. Dla testowanej liczby pierwszej sumy od razu miały większe wartości.
Iloraz n/p dla takich p jest bardzo bliski liczbom całkowitym. Nie zawsze jest to dzielnik, nawet przy zerowaniu się sumy naprzemiennej.
Kolejna cecha namierzajaca to zmiana znaku sumy naprzemiennej. Towarzyszy rozmienieniu wartości na drobne, co występuje także przy dzielniku.
Przy czym nie możemy wchodzić w minima globalne, lecz lokalne sum naprzemiennych, które dla kolejnych podstaw mogą mieć wartości (4, -49, 3 -> dzielnik, -18, 1, 25, 0)

Czy można zatem zaprząc sumy do poszukiwania dzielnika?
Udaje się dla początkowych wartości namierzyć: w tym kierunku jest dzielnik. Lecz przybliżając się do niego za pomocą konwersji:
a[i]*(2^ip^i) + ... + a[1]*(2p) + a[0] = 2^ia[i]*p^i + ... + 2a[1]*p+a[0]
dla coraz większych p (za każdym razem p się podwaja) mamy coraz  większe trudności z utrzymaniem dzielnika 'na celowniku'.
Pojawiające się zerowania lub cechy wspomagające namierzanie dzielnika nie trzymają się jednej wartości, lecz przemieszczają się.

W jednym tescie został namierzony dzielnik dla p=136. Przechodząc do p=272 dzielnik 'zniknął'. Został znowu namierzony dla p=256, przesunięcie 20. Zaś praktycznie dzielnik w tym przypadku był równy 136*8.

Zatem heurystyka niskim kosztem wskazuje: mamy dzielnik. Powinien być bliski takiej to a takiej wartości pomnożonej przez potęgę 2.
Lecz ze znalezieniem należy przeszukiwać brutalnie - chociaż można po wartościach podzielnych przez 4.

Metoda jest hipotezą.

12 kwietnia 2013

Iteracyjne przechodzenie między systemami

Mając faktoryzować liczbę postaci
n = (a*p+b)*p+c,
zastanawiałem się, czy można tak przekształcić tę liczbę, aby szybciej znaleźć jej ewentualne dzielniki. Bez użycia brute force dla poszczególnych liczb pierwszych lub liczb nieparzystych (jak nie znam wszystkich liczb pierwszych).

Pierwsze podejrzenie padło na postać:
n = e*p + (d*p+f)^2 = d*p^2 + (e+2*d*f)*p + f^2
Wzór sugeruje, że jeśli cyfra najmniej znacząca jest kwadratem, to wyłaczając odpowiedni kwadrat z n znajdę mniejszą wartość związaną z dzielnikiem. Pozostaje iloczyn 'e*p', który jednak podczas konwersji na odpowiednią podstawę modyfikuje wyraz wolny.

Kolejne pytanie, czy mogę tak przekształcać liczbę, aby mieć kontrolę nad zmianami wyrazu wolnego?
Wtedy odpowiedni ciąg konwersji zacznie przybliżać wyraz wolny do zera, dla którego istnieje dzielnik.

Podczas faktoryzacji przez zmianę systemów szybko uzyskuję liczby trzycyfrowe, zatem co robią konwersje z wyrazem wolnym liczb trzycyfrowych?
Okazało się, że przy konwersji z systemu o podstawie p na system o podstawie p+k na współczynniki liczb trzycyfrowych znalazłem wzory iteracyjne.
Załóżmy, że zwiększamy podstawę (k dodatnie):
Współczynnik a jest stały.
Współczynnik b zmienia się liniowo od a albo k: b = b-2*a*k. 
Współczynnik c zmienia się kwadratowo od k: c = a*k^2 - b*k + c.
Po wstawieniu nowych wartości
[a,b,c] = [a, b-2ak, akk-bk+c]
uzyskiwałem liczby, które nie były 'cyframi' przedziału [0,p), ale nie zawsze były miejsze, czasem także większe od p. Po 'naprawie' cyfr uzyskiwałem prawidłowe postaci cyfr w kolejnych systemach. 

Ale dobrze widać było tylko sąsiednie podstawy. Kiedy zwiększałem k, nawet tylko do 2, przy większych a bardzo często naprawa modyfikowała b, co powodowało szybko narastający błąd szacowań. Zachowania się wyrazu wolnego nie byłem w stanie przewidzieć.
Zatem ten sposób nie jest najlepszy. Może służyć jako zastąpienie wersji rekursywnej programu wersją iteracyjną (k=2), gdyż wtedy następuje do dwu pożyczek/przeniesień. Z kolei nie wiemy bez sprawdzenia, które z nich będzie.