c - Shanks 的平方因数分解实现

Question

最近我一直在研究Shanks的平方因式分解，就是从这个wiki page开始的。该页面上提供了 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;
}

这是该页面上给定实现的错误吗？该算法是否无法找到某些数字的因子？

Answer 1

维基百科页面中的算法不会检查任何整数除法的分母是否等于零。

我测试了slightly modified version of it (我也使用了C++，这可能是OP实际使用的，因为他们在评论中提到了“C++ STL的gcd函数”的使用)对于N的前32768个值，获得:

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

You'll notice that the first cases are prime numbers and are also present in the multipliers array.

Before the inner loop where the first division is performed, the variable Q is calculated essentially as

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

因此，当 D 是完全平方数时，Q 为零。您可以通过添加几行 code 来处理这些极端情况。 :

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