连分数因子分解法——C语言实现

原创已于 2025-08-17 14:09:58 修改 · 773 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#c语言 #算法 #开发语言

于 2024-09-06 21:55:55 首次发布

参考网址：连分数分解法寻找整数的因子(Python)-优快云博客

连分数的基本概念

连分数是一种表示实数的方式，通过嵌套的分数形式展开。一般形式为：
[ a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + \cdots}}} ]
其中 ( a_0 ) 为整数部分，( a_1, a_2, \ldots ) 为正整数，称为部分分母。若序列有限，则为有限连分数；无限则为无限连分数。

连分数的构造方法

对于实数 ( x )，其连分数展开可通过以下迭代过程生成：

取整数部分 ( a_n = \lfloor x_n \rfloor )，其中初始 ( x_0 = x )。
计算剩余部分 ( x_{n+1} = \frac{1}{x_n - a_n} )。
重复步骤直至 ( x_n ) 为整数或达到所需精度。

示例：对 ( \sqrt{2} ) 的连分数展开：

( a_0 = \lfloor \sqrt{2} \rfloor = 1 )
( x_1 = \frac{1}{\sqrt{2}-1} \approx 2.414 )
( a_1 = \lfloor x_1 \rfloor = 2 )
后续部分重复 ( a_n = 2 )，形成无限循环 ( [1; 2, 2, 2, \ldots] )。

收敛性与近似分数

连分数的截断称为渐近分数（Convergents），记为 ( \frac{p_n}{q_n} )，通过递推关系计算：
[ \begin{cases} p_n = a_n p_{n-1} + p_{n-2} \ q_n = a_n q_{n-1} + q_{n-2} \end{cases} ]
初始条件为 ( p_{-2} = 0, p_{-1} = 1 ) 和 ( q_{-2} = 1, q_{-1} = 0 )。

性质：

渐近分数是最佳有理逼近，即对任意更简单的分数 ( \frac{a}{b} )，有 ( |x - \frac{p_n}{q_n}| \leq |x - \frac{a}{b}| )。
奇数项渐近分数偏小，偶数项偏大，交替逼近真实值。

应用场景

数论与丢番图逼近：用于寻找无理数的有理近似解，如 Pell 方程 ( x^2 - Dy^2 = 1 ) 的求解。
算法优化：在密码学中，RSA 攻击利用连分数分解大整数。
工程计算：硬件设计中使用连分数近似实现复杂函数（如三角函数）的高效计算。

示例：黄金比例的连分数

黄金比例 ( \phi = \frac{1+\sqrt{5}}{2} ) 的连分数展开为：
[ \phi = 1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{1 + \cdots}}} ]
其渐近分数为斐波那契数列的相邻项之比 ( \frac{F_{n+1}}{F_n} )。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <openssl/bn.h>
#include <openssl/err.h>
#include <openssl/evp.h>
#include <openssl/crypto.h>

#define MAX_ITERATIONS 5000000
#define SQRT_CACHE_SIZE 512
#define PROGRESS_INTERVAL 100000
#define EARLY_ABORT_BITS 300

#ifdef _OPENMP
#include <omp.h>
#define MAX_THREADS 4
#else
#define omp_get_thread_num() 0
#endif

typedef struct {
    const char *number;
    const char *description;
    const char *expected_factor;
} TestCase;

// 前向声明所有函数
BN_CTX *global_ctx = NULL;
void init_global_ctx();
void cleanup_global_ctx();
BIGNUM *cached_sqrt(BIGNUM *n);
int montgomery_is_square(BIGNUM *n, BIGNUM **sqrt_result);
int is_square_optimized(BIGNUM *num, BIGNUM **sqrt_result);
void factor_number(const TestCase *test);

// 初始化全局上下文
void init_global_ctx() {
    if (!global_ctx) {
        global_ctx = BN_CTX_new();
        if (!global_ctx) {
            fprintf(stderr, "Error creating BN_CTX\n");
            exit(EXIT_FAILURE);
        }
        BN_CTX_start(global_ctx);
    }
}

// 清理全局上下文
void cleanup_global_ctx() {
    if (global_ctx) {
        BN_CTX_end(global_ctx);
        BN_CTX_free(global_ctx);
        global_ctx = NULL;
    }
}

// 使用Montgomery模乘加速平方验证
int montgomery_is_square(BIGNUM *n, BIGNUM **sqrt_result) {
    BN_CTX *ctx = BN_CTX_new();
    BN_CTX_start(ctx);

    BN_MONT_CTX *mont = BN_MONT_CTX_new();
    BN_MONT_CTX_set(mont, n, ctx);

    BIGNUM *sqrt = cached_sqrt(n);
    BIGNUM *tmp = BN_CTX_get(ctx);

    BN_mod_mul_montgomery(tmp, sqrt, sqrt, mont, ctx);
    BN_from_montgomery(tmp, tmp, mont, ctx);

    int ret = (BN_cmp(tmp, n) == 0);
    if (ret && sqrt_result) {
        *sqrt_result = BN_dup(sqrt);
    }

    BN_MONT_CTX_free(mont);
    BN_CTX_end(ctx);
    BN_CTX_free(ctx);
    return ret;
}

// 优化的平方根计算（带缓存）
BIGNUM *cached_sqrt(BIGNUM *n) {
    static struct {
        BIGNUM *key;
        BIGNUM *value;
    } cache[SQRT_CACHE_SIZE] = {0};
    static int cache_index = 0;

    // 查找缓存
    for (int i = 0; i < SQRT_CACHE_SIZE; i++) {
        if (cache[i].key && BN_cmp(cache[i].key, n) == 0) {
            return BN_dup(cache[i].value);
        }
    }

    // 计算平方根
    BIGNUM *result = BN_new();
    if (BN_is_zero(n) || BN_is_one(n)) {
        BN_copy(result, n);
    } else {
        BN_CTX *ctx = BN_CTX_new();
        BN_CTX_start(ctx);

        BIGNUM *x = BN_CTX_get(ctx);
        BIGNUM *y = BN_CTX_get(ctx);
        BIGNUM *tmp = BN_CTX_get(ctx);
        BIGNUM *two = BN_CTX_get(ctx);
        BN_set_word(two, 2);

        // 优化的初始估计
        int bits = BN_num_bits(n);
        BN_set_bit(x, bits / 2);

        // 牛顿迭代
        for (int i = 0; i < 20; i++) {
            BN_div(y, NULL, n, x, ctx);
            BN_add(y, y, x);
            BN_div(y, NULL, y, two, ctx);

            if (BN_cmp(x, y) == 0) break;
            BN_copy(x, y);
        }
        BN_copy(result, x);

        // 验证并调整
        BN_sqr(tmp, result, ctx);
        if (BN_cmp(tmp, n) > 0) {
            BN_sub(result, result, BN_value_one());
        }

        BN_CTX_end(ctx);
        BN_CTX_free(ctx);
    }

    // 更新缓存
    if (cache_index >= SQRT_CACHE_SIZE) cache_index = 0;
    if (cache[cache_index].key) {
        BN_free(cache[cache_index].key);
        BN_free(cache[cache_index].value);
    }
    cache[cache_index].key = BN_dup(n);
    cache[cache_index].value = BN_dup(result);
    cache_index++;

    return result;
}

// 增强的平方检测
int is_square_optimized(BIGNUM *num, BIGNUM **sqrt_result) {
    // 快速过滤
    static const char last_digit_filter[] = {0,1,4,5,6,9};
    char last_char = BN_bn2dec(num)[BN_num_bytes(num)*2-1];
    int last_digit = last_char - '0';
    int valid = 0;
    for (int i = 0; i < sizeof(last_digit_filter); i++) {
        if (last_digit == last_digit_filter[i]) {
            valid = 1;
            break;
        }
    }
    if (!valid) return 0;

    // 小整数快速路径
    if (BN_num_bits(num) <= 64) {
        unsigned long val = BN_get_word(num);
        unsigned long root = (unsigned long)sqrt(val);
        if (root * root == val) {
            if (sqrt_result) {
                *sqrt_result = BN_new();
                BN_set_word(*sqrt_result, root);
            }
            return 1;
        }
        return 0;
    }

    // 完整计算
    return montgomery_is_square(num, sqrt_result);
}

// 主分解函数
void factor_number(const TestCase *test) {
    printf("\n===== Factoring %s =====\n", test->description);
    printf("Number: %s\n", test->number);

    clock_t start = clock();
    init_global_ctx();

    BIGNUM *n = BN_new();
    BN_dec2bn(&n, test->number);

    BIGNUM *sqrt_n = cached_sqrt(n);
    printf("Square root: %s\n", BN_bn2dec(sqrt_n));

    // 初始化序列变量
    BIGNUM *P = BN_new();
    BIGNUM *Q = BN_new();
    BIGNUM *a = BN_new();
    BIGNUM *p0 = BN_new();
    BIGNUM *p1 = BN_new();
    BIGNUM *p2 = BN_new();
    BIGNUM *tmp1 = BN_new();
    BIGNUM *tmp2 = BN_new();

    BN_zero(P);
    BN_one(Q);
    BN_copy(a, sqrt_n);
    BN_one(p0);
    BN_copy(p1, a);

    int found_factor = 0;
    for (int i = 1; i <= MAX_ITERATIONS && !found_factor; i++) {
        // 提前终止检查
        if (BN_num_bits(Q) > EARLY_ABORT_BITS) {
            printf("Q too large (%d bits), aborting...\n", BN_num_bits(Q));
            break;
        }

        // 保存前一次的Q值
        BIGNUM *Q_prev = BN_dup(Q);

        // 计算P_i
        BN_mul(tmp1, a, Q_prev, global_ctx);
        BN_sub(P, tmp1, P);
        if (BN_is_negative(P)) BN_set_negative(P, 0);

        // 计算Q_i
        BN_sqr(tmp1, P, global_ctx);
        BN_sub(tmp2, n, tmp1);
        BN_div(Q, NULL, tmp2, Q_prev, global_ctx);
        if (BN_is_negative(Q)) BN_set_negative(Q, 0);

        // 计算a_i
        BN_add(tmp1, P, sqrt_n);
        BN_div(a, NULL, tmp1, Q, global_ctx);

        // 更新收敛子序列
        if (i == 1) {
            BN_mul(tmp1, a, p1, global_ctx);
            BN_add(p2, tmp1, p0);
        } else {
            BN_copy(p0, p1);
            BN_copy(p1, p2);
            BN_mul(tmp1, a, p1, global_ctx);
            BN_add(p2, tmp1, p0);
        }

        // 动态调整检查频率
        int check_interval = 1;
        int q_bits = BN_num_bits(Q);
        if (q_bits > 100) check_interval = 100;
        if (q_bits > 150) check_interval = 500;
        if (q_bits > 200) check_interval = 1000;

        if (i % check_interval == 0) {
            BIGNUM *s = NULL;
            if (is_square_optimized(Q, &s)) {
                printf("\n[Step %d] Found square Q (%d bits)\n", i, q_bits);

                // 并行计算GCD
                BIGNUM *f1 = BN_new();
                BIGNUM *f2 = BN_new();

                #pragma omp parallel sections
                {
                    #pragma omp section
                    {
                        BN_add(tmp1, p1, s);
                        BN_gcd(f1, tmp1, n, global_ctx);
                    }

                    #pragma omp section
                    {
                        BN_sub(tmp1, p1, s);
                        BN_gcd(f2, tmp1, n, global_ctx);
                    }
                }

                // 检查因子
                if (!BN_is_one(f1) && BN_cmp(f1, n) != 0) {
                    printf("Factor found: %s\n", BN_bn2dec(f1));
                    if (test->expected_factor &&
                        strcmp(BN_bn2dec(f1), test->expected_factor) == 0) {
                        printf("*** Correct factor verified! ***\n");
                        }
                        found_factor = 1;
                } else if (!BN_is_one(f2) && BN_cmp(f2, n) != 0) {
                    printf("Factor found: %s\n", BN_bn2dec(f2));
                    if (test->expected_factor &&
                        strcmp(BN_bn2dec(f2), test->expected_factor) == 0) {
                        printf("*** Correct factor verified! ***\n");
                        }
                        found_factor = 1;
                }

                BN_free(s);
                BN_free(f1);
                BN_free(f2);
            }
        }

        // 进度报告
        if (i % PROGRESS_INTERVAL == 0) {
            printf("Step %d: Q = %d bits\n", i, q_bits);
        }

        BN_free(Q_prev);
    }

    if (!found_factor) {
        printf("\nNo factor found after %d iterations\n", MAX_ITERATIONS);
    }

    double elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;
    printf("Time elapsed: %.2f seconds\n", elapsed);

    // 清理
    BN_free(n);
    BN_free(sqrt_n);
    BN_free(P);
    BN_free(Q);
    BN_free(a);
    BN_free(p0);
    BN_free(p1);
    BN_free(p2);
    BN_free(tmp1);
    BN_free(tmp2);
}

int main() {
    // 初始化OpenSSL
    OpenSSL_add_all_algorithms();
    ERR_load_crypto_strings();

    // 测试用例
    TestCase tests[] = {
        {"70191551", "70191551 (7793 × 9007)", "7793"},
        {"147573952589676412927", "147573952589676412927 (193707721 × 761838257287)", "193707721"},
        {"74207624142945242263057035287110983967646020057307828709587969646701361764263",
            "256-bit semiprime", NULL}
    };

    // 运行测试
    for (size_t i = 0; i < sizeof(tests)/sizeof(tests[0]); i++) {
        factor_number(&tests[i]);
    }

    // 清理OpenSSL
    EVP_cleanup();
    ERR_free_strings();
    CRYPTO_cleanup_all_ex_data();

    return 0;
}