二次筛法Quadratic Sieve因子分解法----C语言实现

二次筛法概述

二次筛法是一种用于分解大整数的算法，属于现代通用整数分解方法之一。其核心思想是通过寻找模n的平方同余关系，构造出形如x² ≡ y² (mod n)的等式，从而利用n的因数分解性质（若x ≢ ±y mod n，则gcd(x±y, n)可得到非平凡因子）。

数学原理优化

核心公式围绕二次多项式筛选平滑数（smooth numbers）： [ Q(x) = (x + \lfloor \sqrt{n} \rfloor)^2 - n ] 其中x为小整数，目标为找到使Q(x)在选定因子基（factor base）上完全分解的x值。

优化点1：多项式选择

通过调整多项式形式减少计算量： [ Q(x) = ax^2 + 2bx + c ] 满足 ( b^2 - ac = n )，此时： [ aQ(x) = (ax + b)^2 - n ] 通过选择适当的a，使得Q(x)值更小，提高平滑概率。

优化点2：筛分过程

利用对数近似加速筛分： [ \log Q(x) \approx \log(a) + 2\log|x| ] 对每个素数p在因子基中，计算其对数贡献： [ \log p \cdot \left\lfloor \frac{\log Q(x)}{\log p} \right\rfloor ] 通过累加和阈值判断平滑性。

矩阵步骤优化

筛选后得到关系矩阵，优化方法包括：

• 稀疏矩阵技术：利用关系矩阵的稀疏性（如Block Lanczos算法）加速线性代数求解。
• 过滤冗余：提前剔除重复或线性相关的关系式。

参数调优

关键参数选择直接影响效率：

• 因子基大小：B ≈ exp(√(0.5 log n log log n))。
• 筛分区间：M ≈ B²，保证足够平滑数。

注：实际实现需处理多项技术细节（如大整数运算、并行筛分等）。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <string.h>
#include <gmp.h>#define MAX_PRIMES 10000
#define SIEVE_SIZE 10000// 全局变量
uint32_t smooth_relations = 0;
uint32_t required_relations = 0;
uint32_t *factor_base = NULL;
uint32_t factor_base_size = 0;
uint32_t *exponent_vectors = NULL;
int64_t *smooth_x = NULL;
mpz_t N;  // 要分解的数// 辅助函数：检查n是否为素数
int is_prime(uint32_t n) {if (n <= 1) return 0;if (n <= 3) return 1;if (n % 2 == 0 || n % 3 == 0) return 0;for (uint32_t i = 5; i * i <= n; i += 6) {if (n % i == 0 || n % (i + 2) == 0) return 0;}return 1;
}// 计算勒让德符号 (a/p)
int legendre_symbol(mpz_t a, uint32_t p) {if (p == 2) {return 1; // 假设N是奇数}mpz_t tmp, p_mpz;mpz_init(tmp);mpz_init_set_ui(p_mpz, p);mpz_mod(tmp, a, p_mpz);if (mpz_cmp_ui(tmp, 0) == 0) {mpz_clears(tmp, p_mpz, NULL);return 0;}mpz_sub_ui(tmp, p_mpz, 1);mpz_fdiv_q_2exp(tmp, tmp, 1);mpz_powm(tmp, a, tmp, p_mpz);int result = (mpz_cmp_ui(tmp, 1) == 0) ? 1 : -1;mpz_clears(tmp, p_mpz, NULL);return result;
}// 构建因子基
void build_factor_base(uint32_t B) {factor_base = (uint32_t *)malloc(MAX_PRIMES * sizeof(uint32_t));if (!factor_base) {fprintf(stderr, "内存分配失败\n");exit(1);}factor_base[0] = 2; // 总是包含2factor_base_size = 1;// 找到所有小于B的素数p，使得N是模p的二次剩余for (uint32_t p = 3; p <= B && factor_base_size < MAX_PRIMES; p += 2) {if (is_prime(p)) {if (legendre_symbol(N, p) == 1) {factor_base[factor_base_size++] = p;}}}required_relations = factor_base_size + 10; // 需要比因子基大小稍多的关系
}// 多项式Q(x) = (x + floor(sqrt(N)))^2 - N
void poly_Q(mpz_t result, int64_t x, mpz_t sqrtN) {mpz_t tmp;mpz_init(tmp);mpz_set_si(tmp, x);mpz_add(tmp, tmp, sqrtN); // x + floor(sqrt(N))mpz_mul(tmp, tmp, tmp);   // (x + floor(sqrt(N)))^2mpz_sub(result, tmp, N);  // - Nmpz_clear(tmp);
}// Tonelli-Shanks算法求平方根模素数
int tonelli_shanks(mpz_t n, uint32_t p, mpz_t result) {mpz_t Q, S, z, c, t, b, temp, p_mpz;mpz_inits(Q, S, z, c, t, b, temp, p_mpz, NULL);mpz_set_ui(p_mpz, p);if (p == 2) {mpz_set_ui(result, 1);mpz_clears(Q, S, z, c, t, b, temp, p_mpz, NULL);return 1;}// Q * 2^S = p - 1mpz_sub_ui(Q, p_mpz, 1);mpz_set_ui(S, 0);while (mpz_even_p(Q)) {mpz_fdiv_q_2exp(Q, Q, 1);mpz_add_ui(S, S, 1);}// 寻找二次非剩余zmpz_set_ui(z, 2);while (1) {int legendre = legendre_symbol(z, p);if (legendre == -1) break;mpz_add_ui(z, z, 1);}mpz_powm(c, z, Q, p_mpz);mpz_powm(t, n, Q, p_mpz);mpz_add_ui(temp, Q, 1);mpz_fdiv_q_2exp(temp, temp, 1);mpz_powm(result, n, temp, p_mpz);while (mpz_cmp_ui(t, 0) != 0 && mpz_cmp_ui(t, 1) != 0) {uint32_t i = 0;mpz_set(temp, t);while (mpz_cmp_ui(temp, 1) != 0 && i < mpz_get_ui(S)) {mpz_powm_ui(temp, temp, 2, p_mpz);i++;}if (i == 0) {mpz_clears(Q, S, z, c, t, b, temp, p_mpz, NULL);return 0;}mpz_set_ui(temp, 0);mpz_setbit(temp, mpz_get_ui(S) - i - 1);mpz_powm(b, c, temp, p_mpz);mpz_powm_ui(c, b, 2, p_mpz);mpz_mul(t, t, c);mpz_mod(t, t, p_mpz);mpz_mul(result, result, b);mpz_mod(result, result, p_mpz);mpz_set_ui(S, i);}mpz_clears(Q, S, z, c, t, b, temp, p_mpz, NULL);return (mpz_cmp_ui(t, 0) == 0) ? 0 : 1;
}// 筛选过程
void sieve(int64_t sieve_interval, uint32_t *sieve_array) {mpz_t sqrtN, Qx, p_mpz;mpz_init(sqrtN);mpz_init(Qx);mpz_init(p_mpz);mpz_sqrt(sqrtN, N); // floor(sqrt(N))// 初始化筛数组memset(sieve_array, 0, (2 * sieve_interval + 1) * sizeof(uint32_t));// 对因子基中的每个素数p进行筛选for (uint32_t i = 0; i < factor_base_size; i++) {uint32_t p = factor_base[i];mpz_set_ui(p_mpz, p);// 解方程 Q(x) ≡ 0 mod pmpz_t tmp, x1, x2;mpz_init(tmp);mpz_init(x1);mpz_init(x2);mpz_mod(tmp, N, p_mpz); // N mod pif (tonelli_shanks(tmp, p, x1)) { // 如果解存在// 第二个解 x2 = p - x1mpz_sub(x2, p_mpz, x1);// 转换为筛区间内的x值mpz_sub(x1, x1, sqrtN);mpz_mod(x1, x1, p_mpz);mpz_sub(x2, x2, sqrtN);mpz_mod(x2, x2, p_mpz);int64_t r1 = mpz_get_si(x1);int64_t r2 = mpz_get_si(x2);// 在筛区间内标记p的倍数for (int64_t x = r1; x <= sieve_interval; x += p) {if (x >= -sieve_interval) {sieve_array[x + sieve_interval] += log(p) + 0.5;}}for (int64_t x = r2; x <= sieve_interval; x += p) {if (x >= -sieve_interval) {sieve_array[x + sieve_interval] += log(p) + 0.5;}}}mpz_clear(tmp);mpz_clear(x1);mpz_clear(x2);}mpz_clear(sqrtN);mpz_clear(Qx);mpz_clear(p_mpz);
}// 检查Q(x)是否平滑
int is_smooth(mpz_t Qx, uint32_t *exponents) {mpz_t tmp;mpz_init(tmp);mpz_set(tmp, Qx);memset(exponents, 0, factor_base_size * sizeof(uint32_t));for (uint32_t i = 0; i < factor_base_size; i++) {uint32_t p = factor_base[i];while (mpz_divisible_ui_p(tmp, p)) {mpz_divexact_ui(tmp, tmp, p);exponents[i]++;}}int smooth = (mpz_cmp_ui(tmp, 1) == 0);mpz_clear(tmp);return smooth;
}// 高斯消元法模2
void gaussian_elimination() {for (uint32_t col = 0; col < factor_base_size; col++) {// 找到当前列的第一个非零行uint32_t pivot = (uint32_t)-1;for (uint32_t row = col; row < smooth_relations; row++) {if (exponent_vectors[row * factor_base_size + col] % 2 == 1) {pivot = row;break;}}if (pivot == (uint32_t)-1) continue; // 这一列全为零// 交换行if (pivot != col) {for (uint32_t i = 0; i < factor_base_size; i++) {uint32_t temp = exponent_vectors[pivot * factor_base_size + i];exponent_vectors[pivot * factor_base_size + i] = exponent_vectors[col * factor_base_size + i];exponent_vectors[col * factor_base_size + i] = temp;}// 交换x值int64_t temp_x = smooth_x[pivot];smooth_x[pivot] = smooth_x[col];smooth_x[col] = temp_x;}// 消去其他行for (uint32_t row = col + 1; row < smooth_relations; row++) {if (exponent_vectors[row * factor_base_size + col] % 2 == 1) {for (uint32_t i = 0; i < factor_base_size; i++) {exponent_vectors[row * factor_base_size + i] += exponent_vectors[col * factor_base_size + i];exponent_vectors[row * factor_base_size + i] %= 2;}}}}
}// 尝试分解
void try_factorization() {mpz_t sqrtN, a, b, gcd;mpz_init(sqrtN);mpz_init(a);mpz_init(b);mpz_init(gcd);mpz_sqrt(sqrtN, N);for (uint32_t i = 0; i < smooth_relations; i++) {// 构造a和bmpz_set_si(a, smooth_x[i]);mpz_add(a, a, sqrtN);mpz_mod(a, a, N);mpz_set_ui(b, 1);for (uint32_t j = 0; j < factor_base_size; j++) {if (exponent_vectors[i * factor_base_size + j] % 2 == 1) {mpz_mul_ui(b, b, factor_base[j]);}}mpz_mod(b, b, N);// 计算gcd(a-b, N)mpz_sub(gcd, a, b);mpz_gcd(gcd, gcd, N);if (mpz_cmp(gcd, N) != 0 && mpz_cmp_ui(gcd, 1) != 0) {gmp_printf("找到因子: %Zd\n", gcd);mpz_divexact(b, N, gcd);gmp_printf("另一个因子: %Zd\n", b);break;}}mpz_clear(sqrtN);mpz_clear(a);mpz_clear(b);mpz_clear(gcd);
}int main() {mpz_init_set_str(N, "70191551", 10);printf("开始分解N = ");mpz_out_str(stdout, 10, N);printf("\n");// 步骤1: 选择平滑界限B并构建因子基uint32_t B = 100000; // 平滑界限build_factor_base(B);printf("因子基大小: %u\n", factor_base_size);printf("需要的关系数: %u\n", required_relations);// 分配内存smooth_x = (int64_t *)malloc(required_relations * sizeof(int64_t));exponent_vectors = (uint32_t *)malloc(required_relations * factor_base_size * sizeof(uint32_t));uint32_t *sieve_array = (uint32_t *)malloc((2 * SIEVE_SIZE + 1) * sizeof(uint32_t));if (!smooth_x || !exponent_vectors || !sieve_array) {fprintf(stderr, "内存分配失败\n");exit(1);}// 步骤2: 筛选过程sieve(SIEVE_SIZE, sieve_array);// 步骤3: 寻找平滑数mpz_t sqrtN, Qx;mpz_init(sqrtN);mpz_init(Qx);mpz_sqrt(sqrtN, N);uint32_t *exponents = (uint32_t *)malloc(factor_base_size * sizeof(uint32_t));for (int64_t x = -SIEVE_SIZE; x <= SIEVE_SIZE && smooth_relations < required_relations; x++) {// 检查筛数组中的高值if (sieve_array[x + SIEVE_SIZE] > 0.9 * log(mpz_get_d(N))) {poly_Q(Qx, x, sqrtN);if (is_smooth(Qx, exponents)) {// 存储平滑关系smooth_x[smooth_relations] = x;memcpy(&exponent_vectors[smooth_relations * factor_base_size],exponents, factor_base_size * sizeof(uint32_t));smooth_relations++;if (smooth_relations % 10 == 0) {printf("找到 %u 个平滑关系\n", smooth_relations);}}}}printf("总共找到 %u 个平滑关系\n", smooth_relations);// 步骤4: 高斯消元gaussian_elimination();// 步骤5: 尝试分解try_factorization();// 清理free(factor_base);free(smooth_x);free(exponent_vectors);free(sieve_array);free(exponents);mpz_clear(N);mpz_clear(sqrtN);mpz_clear(Qx);return 0;
}