Реализация факторизации квадратной формы Шенкса - PullRequest
Новая
       34

Реализация факторизации квадратной формы Шенкса

0 голосов
/ 10 октября 2018

Недавно я изучал факторизацию квадратной формы Шенкса с этой вики-страницы На этой странице представлена ​​реализация на C.Я тестировал эту функцию и заметил, что функция не может найти коэффициент 27.

Это заданная функция C: 

#include <inttypes.h>
#define nelems(x) (sizeof(x) / sizeof((x)[0]))

const int multiplier[] = {1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11};

uint64_t SQUFOF( uint64_t N )
{
    uint64_t D, Po, P, Pprev, Q, Qprev, q, b, r, s;
    uint32_t L, B, i;
    s = (uint64_t)(sqrtl(N)+0.5);
    if (s*s == N) return s;
    for (int k = 0; k < nelems(multiplier) && N <= UINT64_MAX/multiplier[k]; k++) {
        D = multiplier[k]*N;
        Po = Pprev = P = sqrtl(D);
        Qprev = 1;
        Q = D - Po*Po;
        L = 2 * sqrtl( 2*s );
        B = 3 * L;
        for (i = 2 ; i < B ; i++) {
            b = (uint64_t)((Po + P)/Q);
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            r = (uint64_t)(sqrtl(Q)+0.5);
            if (!(i & 1) && r*r == Q) break;
            Qprev = q;
            Pprev = P;
        };
        if (i >= B) continue;
        b = (uint64_t)((Po - P)/r);
        Pprev = P = b*r + P;
        Qprev = r;
        Q = (D - Pprev*Pprev)/Qprev;
        i = 0;
        do {
            b = (uint64_t)((Po + P)/Q);
            Pprev = P;
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            Qprev = q;
            i++;
        } while (P != Pprev);
        r = gcd(N, Qprev);
        if (r != 1 && r != N) return r;
    }
    return 0;
}

Это ошибка данной реализации наэта страница?Может ли этот алгоритм не найти фактор для некоторых чисел?

1 Ответ

0 голосов
/ 11 октября 2018

Алгоритм на странице Википедии не проверяет, равны ли знаменатели какого-либо из целочисленных делений нулю.

Я протестировал слегка измененную версию этого (я такжеиспользовал C ++, что, вероятно, на самом деле и использовал OP, поскольку в комментариях упоминается об использовании «функции gcd C ++ STL») для первых 32768 значений N, получая: 

N:      3  Q == 0 in the first division
N:      5  Q == 0 in the first division
N:      7  Q == 0 in the first division
N:     11  Q == 0 in the first division
N:     27  Q == 0 in the first division
N:    363  Q == 0 in the first division
N:    867  Q == 0 in the first division
N:   1445  Q == 0 in the first division
N:   5043  Q == 0 in the first division
N:   6845  Q == 0 in the first division
N:   7803  Q == 0 in the first division
N:  10443  Q == 0 in the first division
N:  11163  Q == 0 in the first division
N:  13467  Q == 0 in the first division
N:  14283  Q == 0 in the first division
N:  18491  Q == 0 in the first division
N:  18723  Q == 0 in the first division
N:  23763  Q == 0 in the first division
N:  30603  Q == 0 in the first division
N:  31827  Q == 0 in the first division

You 'Заметим, что первые случаи являются простыми числами и также присутствуют в массиве multipliers.

Перед внутренним циклом, в котором выполняется первое деление, переменная Q вычисляется по существу как 

D = multiplier[k]*N;
Po = sqrtl(D);
Q = D - Po*Po;

Итак, Q - ноль, когда D - идеальный квадрат.Вы можете позаботиться об этих угловых случаях, добавив несколько строк код : 

const uint64_t multiplier[] = {
    1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7,
    3*5*11, 3*7*11, 5*7*11, 3*5*7*11
};
const uint64_t results[] = {
    1, 0, 0, 0, 0, 3, 3, 3, 5, 5, 7, 3, 3, 3, 5, 3
};
    // ... inside the k loop
    if ( multiplier[k] == N )   
        return results[k];
    // ... calculate Q as before
    if (Q == 0) 
        return std::gcd(multiplier[k], N);
    // ... rest of the loop
...