椭圆曲线ECM因子分解法——C语言实现

原创已于 2025-08-14 17:19:54 修改 · 855 阅读

20 ·

CC 4.0 BY-SA版权

文章标签：

#c语言 #前端 #javascript

于 2024-08-31 16:58:20 首次发布

ECM因子分解法概述

椭圆曲线方法（Elliptic Curve Method，ECM）是一种用于整数分解的算法，由Hendrik Lenstra于1985年提出。该方法通过椭圆曲线上的点运算寻找非平凡因子，尤其适用于中等大小的合数（通常为30到60位）。其核心思想是将Pollard的p-1算法推广到椭圆曲线群上，利用群的阶随曲线参数变化的特性提高分解效率。

ECM算法步骤

选择椭圆曲线和初始点
随机选择一条椭圆曲线 ( E: y^2 = x^3 + ax + b \mod n ) 和一个曲线上的点 ( P )。参数 ( a, b ) 需满足 ( 4a^3 + 27b^2 \neq 0 \mod n )，否则曲线奇异。

计算倍点
选择一个平滑边界 ( B )，计算 ( k = \text{lcm}(1, 2, \dots, B) )，然后计算 ( k \cdot P )（标量乘法）。过程中涉及模逆运算，若无法计算逆元且结果不为 ( n )，则找到 ( n ) 的一个非平凡因子。

检查因子
若标量乘法过程中出现模逆失败，且与 ( n ) 的最大公约数 ( d ) 满足 ( 1 < d < n )，则 ( d ) 为 ( n ) 的一个因子。否则需更换曲线或调整参数重新尝试。

关键点说明

随机曲线生成：通过随机参数 ( a ) 和点 ( P = (x, y) ) 构造曲线，确保每次尝试的群阶不同。
平滑边界选择：( B ) 的选择影响计算效率和成功率，通常根据目标数大小调整。
标量乘法优化：实际实现需使用高效的倍点算法（如Montgomery阶梯），避免逐次加法。

注意事项

GMP库依赖：示例需链接GMP库（-lgmp）编译。
简化逻辑：示例未完整实现椭圆曲线点运算，仅展示核心流程。实际应用需补充点加、倍点等函数。
边界处理：需检查曲线非奇异性和点的有效性，避免无效运算。

通过调整参数和优化实现，ECM可高效分解中等规模的合数，是数论和密码学中的重要工具。

椭圆曲线分解关键点：

算法核心思想：
- 在模n的椭圆曲线上进行点运算
- 当计算斜率需要模逆元时，如果分母与n不互质，则找到n的因子
标量乘法优化：
- 使用倍加算法高效计算k*P
- k选择为不超过B的所有素数幂的乘积
参数选择：
- 光滑界限B决定计算量和成功率
- 70191551的因子是7793和9007，需要B≥10001
- 实际中B=1000足够，因为算法可能在较小k值就找到因子

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

typedef struct {
    long long x;
    long long y;
} Point;

typedef struct {
    int found_factor;
    long long factor;
    Point result;
} PointOpResult;

// 前置声明函数
long long gcd(long long a, long long b);
long long extended_gcd(long long a, long long b, long long *x, long long *y);
long long mod_inverse(long long a, long long n, long long *g);
PointOpResult point_double(Point P, long long a, long long n);  // 添加前置声明
PointOpResult point_add(Point P, Point Q, long long a, long long n);
PointOpResult scalar_mult(Point P, long long k, long long a, long long n);
int generate_primes(int B, int *primes);
void compute_prime_powers(int B, int *primes, int num_primes, long long *prime_powers);
long long ecm_factor(long long n, int B, int max_tries);

long long gcd(long long a, long long b) {
    while (b) {
        long long t = b;
        b = a % b;
        a = t;
    }
    return a;
}

long long extended_gcd(long long a, long long b, long long *x, long long *y) {
    if (b == 0) {
        *x = 1;
        *y = 0;
        return a;
    }
    long long g = extended_gcd(b, a % b, y, x);
    *y -= a / b * *x;
    return g;
}

long long mod_inverse(long long a, long long n, long long *g) {
    long long x, y;
    *g = extended_gcd(a, n, &x, &y);
    if (*g != 1) {
        return 0;
    }
    return (x % n + n) % n;
}

PointOpResult point_double(Point P, long long a, long long n) {
    PointOpResult res;
    res.found_factor = 0;

    if (P.x == 0 && P.y == 0) {
        res.result = P;
        return res;
    }
    if (P.y == 0) {
        res.result.x = 0;
        res.result.y = 0;
        return res;
    }

    long long denom = (2 * P.y) % n;
    long long g;
    long long inv = mod_inverse(denom, n, &g);
    if (g != 1) {
        res.found_factor = 1;
        res.factor = g;
        return res;
    }

    long long lambda = (3 * P.x * P.x + a) % n;
    if (lambda < 0) lambda += n;
    lambda = (lambda * inv) % n;
    long long x3 = (lambda * lambda - 2 * P.x) % n;
    if (x3 < 0) x3 += n;
    long long y3 = (lambda * (P.x - x3) - P.y) % n;
    if (y3 < 0) y3 += n;

    res.result.x = x3;
    res.result.y = y3;
    return res;
}

PointOpResult point_add(Point P, Point Q, long long a, long long n) {
    PointOpResult res;
    res.found_factor = 0;

    // 处理无穷远点（单位元）
    if (P.x == 0 && P.y == 0) {
        res.result = Q;
        return res;
    }
    if (Q.x == 0 && Q.y == 0) {
        res.result = P;
        return res;
    }

    // 处理相同点或互为逆元的情况
    if (P.x == Q.x) {
        if (P.y == Q.y) {
            // P = Q，执行倍点运算
            return point_double(P, a, n);
        } else {
            // P和Q是y坐标相反的互逆点，返回无穷远点
            res.result.x = 0;
            res.result.y = 0;
            return res;
        }
    }

    // 计算斜率
    long long dx = (Q.x - P.x) % n;
    if (dx < 0) dx += n;
    long long g;
    long long inv = mod_inverse(dx, n, &g);
    if (g != 1) {
        res.found_factor = 1;
        res.factor = g;
        return res;
    }

    long long dy = (Q.y - P.y) % n;
    if (dy < 0) dy += n;
    long long lambda = (dy * inv) % n;
    long long x3 = (lambda * lambda - P.x - Q.x) % n;
    if (x3 < 0) x3 += n;
    long long y3 = (lambda * (P.x - x3) - P.y) % n;
    if (y3 < 0) y3 += n;

    res.result.x = x3;
    res.result.y = y3;
    return res;
}

PointOpResult scalar_mult(Point P, long long k, long long a, long long n) {
    PointOpResult res;
    res.found_factor = 0;
    Point R = {0, 0};  // 初始化为无穷远点
    Point T = P;

    // 使用倍加算法计算k*P
    while (k) {
        if (k & 1) {
            PointOpResult add_res = point_add(R, T, a, n);
            if (add_res.found_factor) {
                return add_res;
            }
            R = add_res.result;
        }
        k >>= 1;
        if (k == 0) break;

        PointOpResult double_res = point_double(T, a, n);
        if (double_res.found_factor) {
            return double_res;
        }
        T = double_res.result;
    }

    res.result = R;
    return res;
}

int generate_primes(int B, int *primes) {
    if (B < 2) return 0;
    int *is_prime = (int *)malloc((B + 1) * sizeof(int));
    for (int i = 2; i <= B; i++) is_prime[i] = 1;
    
    // 埃拉托斯特尼筛法
    for (int i = 2; i * i <= B; i++) {
        if (is_prime[i]) {
            for (int j = i * i; j <= B; j += i) {
                is_prime[j] = 0;
            }
        }
    }
    
    // 收集素数
    int count = 0;
    for (int i = 2; i <= B; i++) {
        if (is_prime[i]) {
            primes[count++] = i;
        }
    }
    free(is_prime);
    return count;
}

void compute_prime_powers(int B, int *primes, int num_primes, long long *prime_powers) {
    for (int i = 0; i < num_primes; i++) {
        int p = primes[i];
        long long power = p;
        long long next = power * p;
        // 计算不超过B的最大素数幂
        while (next <= B) {
            power = next;
            next = power * p;
        }
        prime_powers[i] = power;
    }
}

long long ecm_factor(long long n, int B, int max_tries) {
    srand(time(NULL));
    for (int try_count = 0; try_count < max_tries; try_count++) {
        // 随机生成曲线参数和起点
        long long a = rand() % n;
        long long x = rand() % n;
        long long y = rand() % n;
        long long b = (y * y - x * x * x - a * x) % n;
        if (b < 0) b += n;

        // 检查曲线是否非奇异
        long long disc = (4 * a * a * a + 27 * b * b) % n;
        if (disc < 0) disc += n;
        long long g = gcd(disc, n);
        if (g > 1) {
            if (g < n) return g;  // 找到因子
            continue;
        }

        // 生成不超过B的素数表
        int *primes = (int *)malloc(B * sizeof(int));
        int num_primes = generate_primes(B, primes);
        long long *prime_powers = (long long *)malloc(num_primes * sizeof(long long));
        compute_prime_powers(B, primes, num_primes, prime_powers);

        // 计算复合阶k = ∏ p^e(p,B)
        Point P = {x, y};
        Point R = P;

        // 逐个应用素数幂
        for (int i = 0; i < num_primes; i++) {
            PointOpResult mult_res = scalar_mult(R, prime_powers[i], a, n);
            if (mult_res.found_factor) {
                free(primes);
                free(prime_powers);
                return mult_res.factor;  // 找到因子
            }
            R = mult_res.result;
        }

        free(primes);
        free(prime_powers);
    }
    return 0;  // 未找到因子
}

int main() {
    long long n = 70191551;
    int B = 1000;  // 增大光滑界限以提高成功率
    int max_tries = 100;

    printf("Attempting to factor %lld using ECM...\n", n);
    clock_t start = clock();
    long long factor = ecm_factor(n, B, max_tries);
    double elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;

    if (factor > 0 && factor < n) {
        printf("Factorization found: %lld = %lld * %lld\n", n, factor, n / factor);
        printf("Time taken: %.3f seconds\n", elapsed);
    } else {
        printf("Failed to find a factor after %d tries.\n", max_tries);
    }

    return 0;
}