After longer discussion with Gemini on Lenstra's ECM factoring method I got sample code that worked ...

ecm.cpp

#include <iostream>
#include <vector>
#include <string>
#include <stdexcept>
#include <chrono>
#include <gmpxx.h>
#include <gmp.h>

class FactorFoundException : public std::exception {
public:
    mpz_class factor;
    FactorFoundException(mpz_class f) : factor(f) {}
};

struct ECPoint {
    mpz_class x, y;
    bool is_infinity;

    ECPoint() : x(0), y(0), is_infinity(true) {}
    ECPoint(mpz_class x_in, mpz_class y_in) : x(x_in), y(y_in), is_infinity(false) {}
};

mpz_class mod_N(mpz_class val, const mpz_class& N) {
    val %= N;
    if (val < 0) val += N;
    return val;
}

ECPoint elladd(const ECPoint& P, const ECPoint& Q, const mpz_class& a, const mpz_class& N) {
    if (P.is_infinity) return Q;
    if (Q.is_infinity) return P;

    mpz_class num, den, inv_den;
    mpz_class diff_x = mod_N(P.x - Q.x, N);

    if (diff_x == 0) {
        mpz_class diff_y = mod_N(P.y - Q.y, N);
        if (diff_y == 0) {
            num = 3 * P.x * P.x + a;
            den = 2 * P.y;
        } else {
            return ECPoint(); 
        }
    } else {
        num = Q.y - P.y;
        den = diff_x;
    }

    num = mod_N(num, N);
    den = mod_N(den, N);

    if (mpz_invert(inv_den.get_mpz_t(), den.get_mpz_t(), N.get_mpz_t()) == 0) {
        mpz_class g;
        mpz_gcd(g.get_mpz_t(), den.get_mpz_t(), N.get_mpz_t());
        if (g > 1 && g < N) {
            throw FactorFoundException(g); 
        }
        return ECPoint(); 
    }

    mpz_class lambda = mod_N(num * inv_den, N);
    mpz_class x3 = mod_N(lambda * lambda - P.x - Q.x, N);
    mpz_class y3 = mod_N(lambda * (P.x - x3) - P.y, N);

    return ECPoint(x3, y3);
}

ECPoint ellmul(ECPoint P, const mpz_class& k, const mpz_class& a, const mpz_class& N) {
    ECPoint R; 
    ECPoint Q = P;
    long bits = mpz_sizeinbase(k.get_mpz_t(), 2);

    for (long i = 0; i < bits; ++i) {
        if (mpz_tstbit(k.get_mpz_t(), i)) {
            R = elladd(R, Q, a, N);
        }
        Q = elladd(Q, Q, a, N);
    }
    return R;
}

mpz_class generate_B_smooth_k(unsigned long B) {
    mpz_class k = 1;
    std::vector<bool> is_prime(B + 1, true);
    
    for (unsigned long p = 2; p <= B; ++p) {
        if (is_prime[p]) {
            for (unsigned long i = p * p; i <= B; i += p) {
                is_prime[i] = false;
            }
            unsigned long power = p;
            while (power * p <= B) power *= p;
            k *= power;
        }
    }
    return k;
}

int main(int argc, char* argv[]) {
    if (argc < 3) {
        std::cerr << "Usage: " << argv[0] << " <N> <m> [B]" << std::endl;
        return 1;
    }

    mpz_class N(argv[1]);
    int m = std::stoi(argv[2]);
    unsigned long B = (argc > 3) ? std::stoul(argv[3]) : 97;

    std::cout << "Factoring N = " << N << std::endl;

    // Timing Multiplier Generation
    auto t1_k = std::chrono::high_resolution_clock::now();
    mpz_class k = generate_B_smooth_k(B);
    auto t2_k = std::chrono::high_resolution_clock::now();
    auto duration_k = std::chrono::duration_cast<std::chrono::milliseconds>(t2_k - t1_k).count();

    std::cout << "Generated " << B << "-smooth multiplier k (" 
              << mpz_sizeinbase(k.get_mpz_t(), 2) << " bits) in " 
              << duration_k << " ms" << std::endl;

    gmp_randclass rand(gmp_randinit_default);
    rand.seed(time(NULL));

    // ECM Loop: Try m different curves
    for (int i = 1; i <= m; ++i) {
        mpz_class x0 = rand.get_z_range(N - 1) + 1;
        mpz_class y0 = rand.get_z_range(N - 1) + 1;
        mpz_class a  = rand.get_z_range(N - 1) + 1;
        ECPoint P(x0, y0);
        
        // Timing Current Curve
        auto t1_curve = std::chrono::high_resolution_clock::now();
        
        try {
            ellmul(P, k, a, N);
            
            auto t2_curve = std::chrono::high_resolution_clock::now();
            auto duration_curve = std::chrono::duration_cast<std::chrono::milliseconds>(t2_curve - t1_curve).count();
            
            std::cout << "Curve " << i << "/" << m << " failed to find a factor (" 
                      << duration_curve << " ms)" << std::endl;

        } catch (const FactorFoundException& e) {
            auto t2_curve = std::chrono::high_resolution_clock::now();
            auto duration_curve = std::chrono::duration_cast<std::chrono::milliseconds>(t2_curve - t1_curve).count();
            
            std::cout << "\n*** SUCCESS on Curve " << i << " (" << duration_curve << " ms) ***" << std::endl;
            std::cout << "Factor found: " << e.factor << std::endl;
            std::cout << "Cofactor:     " << (N / e.factor) << std::endl;
            return 0;
        }
    }

    std::cout << "\nFailed to factor N after " << m << " curves." << std::endl;
    return 0;
}

... but is much slower than gmp-ecm. So I asked why.

my last prompt

works, but takes time:

hermann@8840hs:~/uni-heidelberg/scripts/stuff$ ./ecm 49999999464999998803 1000 10000000
Factoring N = 49999999464999998803
Generated 10000000-smooth multiplier k (14424844 bits) in 39379 ms
Curve 1/1000 failed to find a factor (14931 ms)
Curve 2/1000 failed to find a factor (14923 ms)
Curve 3/1000 failed to find a factor (14957 ms)
...
Curve 521/1000 failed to find a factor (14900 ms)
Curve 522/1000 failed to find a factor (14907 ms)
Curve 523/1000 failed to find a factor (14852 ms)

*** SUCCESS on Curve 524 (13342 ms) ***
Factor found: 4999999937
Cofactor:     10000000019
hermann@8840hs:~/uni-heidelberg/scripts/stuff$ 

Why is gmp-ecm so much faster?

hermann@8840hs:~/uni-heidelberg/scripts/stuff$ echo 49999999464999998803 | ecm -c 100 1629
GMP-ECM 7.0.5 [configured with GMP 6.3.0, --enable-asm-redc] [ECM]
Input number is 49999999464999998803 (20 digits)
Using B1=1629, B2=109146, polynomial x^1, sigma=1:175852517
Step 1 took 0ms
Step 2 took 1ms
********** Factor found in step 2: 4999999937
Found prime factor of 10 digits: 4999999937
Prime cofactor 10000000019 has 11 digits
hermann@8840hs:~/uni-heidelberg/scripts/stuff$ 

Gemini response

The 14-second vs. 1-millisecond discrepancy between your script and gmp-ecm is a masterclass in the difference between a textbook mathematical implementation and a production-grade cryptographic factoring engine.

Here is exactly why gmp-ecm bypassed your runtime:

• Massive Over-Computation of the Multiplier: You commanded your script to compute a strictly $10,000,000$-smooth multiplier, forcing it to generate and process a 14.4-million-bit scalar $k$. gmp-ecm intelligently evaluated your 20-digit input $N$, expected a ~10-digit prime factor, and automatically selected a much smaller, optimal bound of $B_1 = 1629$.

• The Missing Stage 2: Your C++ script only implements Stage 1 ECM. If the elliptic curve group order modulo $p$ contains even one prime factor larger than your bound $B$, your script fails. gmp-ecm implements a standard Stage 2 search. After clearing small prime factors up to $B_1=1629$, it uses a highly optimized Baby-Step Giant-Step algorithm to search for a single remaining large prime factor up to $B_2 = 109,146$. The log output Factor found in step 2 proves the factor was caught in this exact net.

• Projective Coordinates (Zero Inner-Loop Inversions): Your elladd function executes a full mpz_invert (the Extended Euclidean Algorithm) for the tangent/secant denominator on every single bit of the scalar multiple. Just as you experienced when tuning your proth_prover pipeline, forcing the CPU to execute millions of dynamic GMP heap allocations and GCD inversions inside an inner loop causes catastrophic slowdowns. gmp-ecm avoids this by using Projective Coordinates $(X, Y, Z)$. By representing points as fractions, it replaces all division with basic multiplication during the loop, performing only one single modular inversion at the very end of the exponentiation.

• Montgomery Reduction (redc): Your script relies on the standard modulo operator (% N) to constrain coordinate sizes, which requires slow multiprecision integer division. The log [configured with GMP 6.3.0, --enable-asm-redc] shows gmp-ecm translates the curve parameters into the Montgomery domain. This replaces integer division with ultra-fast bitwise shifts and masked limb multiplications.

• Optimized Curve Selection: Instead of random short Weierstrass curves, gmp-ecm defaults to generating Montgomery curves or Suyama curves. Point doubling and addition on these specific geometries require significantly fewer floating-point multiplications per bit than the general $y^2 = x^3 + ax + b$ form.

To achieve times comparable to gmp-ecm in C++, you would need to implement Stage 2, transition your coordinate geometry to $X/Z$ projective planes, and utilize a fast-reduction algorithm like Montgomery or Barrett reduction to bypass the mpz_invert bottleneck.