平方剩余因子分解法SQUFOF—C语言实现

算法概述：
  SQUFOF（Shanks' method of square forms factorization）即香克斯平方形式分解算法，是由数学家Daniel Shanks提出的用于分解大整数的一个因子分解算法。SQUFOF算法基于平方形式的连分数展开理论，通过寻找某些特殊形式的二次方程的解来找出整数的因子。这种算法尤其适用于分解较小的因子。

关键概念：
1. 平方形式（Square Form）：一个整数的平方形式是指能够写成a^2 - nb^2形式的整数，其中n是一个给定的正整数。在SQUFOF算法中，寻找这样的平方形式是分解因子的关键步骤。
2. 连分数（Continued Fraction）：连分数是一种特殊的分数表达形式，用来逼近实数。它由一系列的整数分子和分母构成，每个分子是前一个分数的分母，分母是两个整数的差。连分数在数论和代数中有着广泛的应用。
3. 二次互反律（Quadratic Reciprocity Law）：这是数论中的一个重要定理，用于判定一个数在模另一个数的乘法群中是否可逆，即是否是模另一个数的二次剩余。

SQUFOF算法原理：
香克斯平方形式分解算法利用二次数的特性，通过以下步骤实现因子分解：
1. 选取合适n值  ：从待分解的整数N中选取一个合适的n值，构成平方形式N = a^2 - nb^2。
2. 构造二次序列 ：利用n值构造一系列二次数的序列，并进行连分数展开。
3. 寻找周期性   ：在连分数展开过程中寻找周期性，利用周期性寻找二次方程的解。
4. 分解因子     ：找到满足条件的x、y值，使得N可以分解为(N = x^2 - ny^2)，从而得到N的因子。

#include <math.h>
#include <stdio.h>
 
#define nelems(x) (sizeof(x) / sizeof((x)[0]))
 
const unsigned long 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};
 
unsigned long long gcd(unsigned long long a, unsigned long long b)
{
    while (b != 0)
    {
        a %= b;
        a ^= b;
        b ^= a;
        a ^= b;
    }
 
    return a;
}
 
unsigned long long SQUFOF( unsigned long long N )
{
    unsigned long long D, Po, P, Pprev, Q, Qprev, q, b, r, s;
    unsigned long L, B, i;
    s = (unsigned long long)(sqrtl(N)+0.5);
    if (s*s == N) return s;
    for (int k = 0; k < nelems(multiplier) && N <= 0xffffffffffffffff/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 = (unsigned long long)((Po + P)/Q);
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            r = (unsigned long long)(sqrtl(Q)+0.5);
            if (!(i & 1) && r*r == Q) break;
            Qprev = q;
            Pprev = P;
        };
        if (i >= B) continue;
        b = (unsigned long long)((Po - P)/r);
        Pprev = P = b*r + P;
        Qprev = r;
        Q = (D - Pprev*Pprev)/Qprev;
        i = 0;
        do {
            b = (unsigned long long)((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;
}
 
int main(int argc, char *argv[]) {
    /*
    int i;
    const unsigned long long data[] = {
        2501,
        12851,
        13289,
        75301,
        120787,
        967009,
        997417,
        7091569,
        70191551,
        13290059,
        42854447,
        223553581,
        2027651281,
        11111111111,
        100895598169,
        1002742628021,
        60012462237239,
        287129523414791,
        9007199254740931,
        11111111111111111,
        314159265358979323,
        384307168202281507,
        419244183493398773,
        658812288346769681,
        922337203685477563,
        1000000000000000127,
        1152921505680588799,
        1537228672809128917,
        4611686018427387877};
    for(int i = 0; i < nelems(data); i++) {
        unsigned long long example, factor, quotient;
        example = data[i];
        factor = SQUFOF(example);
        if(factor == 0) {
            printf("%llu was not factored.\n", example);
        }
        else {
            quotient = example / factor;
            printf("Integer %llu has factors %llu and %llu\n",
               example, factor, quotient); 
        }
    }*/
    
    ;
    printf("q =  %ld\n",SQUFOF(70191551));
}