factor128.c

#include "ptypes.h"
#include "constants.h"
#include "factor.h"
#include "factor128.h"
#include "tinysiqs128.h"
#include "cache.h"
#include "sieve.h"
#include "util.h"       /* for verbose and next_prime */
#include "util_math.h"  /* for isqrt */
#include "moebius.h"

#if HAVE_FACTOR128

int u128_to_str(char str[40], uint128_t n) {
  int slen = 0;

  do {
    uint128_t d = n / 10;
    str[slen++] = '0' + (char)(n - d*10);
    n = d;
  } while (n);

  str[slen] = '\0';
  if (slen > 1) {
    char *L = str, *R = str+slen;
    while (--R > L) { char t = *R; *R = *L; *L++ = t; } /* Reverse digits. */
  }
  return slen;
}

bool str_to_u128(uint128_t *out, const char *s, size_t len)
{
  static const char uint128_max_str[] = "340282366920938463463374607431768211455";
  size_t i;
  uint128_t n = 0;

  if (len == 0) return 0;
  if (*s == '+') { s++; len--; if (len == 0) return 0; }
  if (*s == '-') return 0;

  while (len > 1 && *s == '0') { s++; len--; }

  if (len > 39 || (len == 39 && memcmp(s, uint128_max_str, 39) > 0))
    return 0;

  for (i = 0; i < len; i++) {
    if (s[i] < '0' || s[i] > '9')
      return 0;
    n = n * 10 + (uint8_t)(s[i] - '0');
  }

  *out = n;
  return 1;
}

/* Decimal string for a uint128_t, for debugging.
 * printf("  found factor %s of %s\n", u128_str(f), u128_str(t)); */
static const char *u128_str(uint128_t n) {
  static char bufs[4][42];
  static int  idx = 0;
  char *buf = bufs[idx++ & 3];
  u128_to_str(buf, n);
  return buf;
}

/*****************************************************************************
 * 128-bit modular arithmetic
 * All functions assume 0 <= a, b < n and n > 1.
 *****************************************************************************/

static INLINE uint128_t addmod128(uint128_t a, uint128_t b, uint128_t n) {
  /* Avoid relying on unsigned 128-bit overflow. */
  uint128_t nb = n - b;
  return (a >= nb) ? a - nb : a + b;
}

static INLINE uint128_t submod128(uint128_t a, uint128_t b, uint128_t n) {
  return (a >= b) ? a - b : n - (b - a);
}

/* (a + n) / 2 mod n, for odd n */
static INLINE uint128_t half_mod128(uint128_t a, uint128_t n) {
  return (a & 1)  ?  (a >> 1) + (n >> 1) + 1  :  a >> 1;
}

/* Binary double-and-add mulmod.  O(128) addmod128 calls.
 * Used only for primality testing (variable modulus). */
static uint128_t mulmod128(uint128_t a, uint128_t b, uint128_t n) {
  uint128_t r = 0;

  if (a >= n) a %= n;
  if (b >= n) b %= n;

  if (a <= (uint128_t)UINT64_MAX && b <= (uint128_t)UINT64_MAX)
    return ((uint128_t)(uint64_t)a * (uint64_t)b) % n;

  uint128_t x = a >= b ? a : b;
  uint128_t y = a >= b ? b : a;

  while (y > 0) {
    if (y & 1) r = addmod128(r, x, n);
    x = addmod128(x, x, n);
    y >>= 1;
  }
  return r;
}

#define sqrmod128(a, n)  mulmod128(a, a, n)

#if 0
static uint128_t powmod128(uint128_t a, uint128_t k, uint128_t n) {
  uint128_t r = 1;
  a %= n;
  while (k > 0) {
    if (k & 1) r = mulmod128(r, a, n);
    k >>= 1;
    if (k) a = sqrmod128(a, n);
  }
  return r;
}
#endif

static uint128_t pow2mod128(unsigned int e, uint128_t n) {
  uint128_t r = 1;
  while (e-- > 0)
    r = addmod128(r, r, n);
  return r;
}

/* Inlining this gives problems with gcc 14 and _BitInt(128) */
static NOINLINE uint128_t gcd128(uint128_t a, uint128_t b) {
  while (b) { uint128_t t = b; b = a % b; a = t; }
  return a;
}

/* Returns a^{-1} mod n, or 0 if gcd(a,n) > 1. */
static uint128_t modinv128(uint128_t a, uint128_t n) {
  uint128_t r0 = n, r1 = a % n;
  uint128_t t0 = 0, t1 = 1;    /* coefficients modulo n */

  if (r1 == 0) return 0;

  while (r1) {
    uint128_t q = r0 / r1;
    uint128_t r = r0 - q * r1;
    uint128_t qt = q % n;
    if (t1 != 1) qt = mulmod128(qt, t1, n);
    uint128_t t = submod128(t0, qt, n);

    r0 = r1;  r1 = r;
    t0 = t1;  t1 = t;
  }

  return (r0 == 1) ? t0 : 0;
}


/*****************************************************************************
 * Montgomery multiplication for fixed modulus.
 *
 * R = 2^128.  Given odd n, we keep values in Montgomery form: â = a*R mod n.
 * mont_mulmod(â, b̂) = â*b̂*R^{-1} mod n  (REDC algorithm).
 * Correct for all odd n < 2^128, including n > 2^127.
 *
 * Setup: ninv = -n^{-1} mod 2^128 (via Hensel lifting).
 *        r2   = R^2 mod n          (for converting to Montgomery form).
 *****************************************************************************/

typedef struct { uint128_t n, ninv, r2; } mont128_t;

/* Low 128 bits of the full 256-bit product a*b. */
static INLINE uint128_t lo128(uint128_t a, uint128_t b) {
  uint64_t a0 = (uint64_t)a,  a1 = (uint64_t)(a >> 64);
  uint64_t b0 = (uint64_t)b,  b1 = (uint64_t)(b >> 64);
  uint128_t lo = (uint128_t)a0 * b0;
  lo += (uint128_t)a0 * b1 << 64;
  lo += (uint128_t)a1 * b0 << 64;
  return lo;
}

/* Full 256-bit product a*b, split into lo/hi 128-bit halves. */
static INLINE void mul256(uint128_t a, uint128_t b,
                          uint128_t *lo_out, uint128_t *hi_out) {
  uint64_t a0 = (uint64_t)a,  a1 = (uint64_t)(a >> 64);
  uint64_t b0 = (uint64_t)b,  b1 = (uint64_t)(b >> 64);

  uint128_t p00 = (uint128_t)a0 * b0;
  uint128_t p01 = (uint128_t)a0 * b1;
  uint128_t p10 = (uint128_t)a1 * b0;
  uint128_t p11 = (uint128_t)a1 * b1;

  /* lo = p00 + ((p01 + p10) << 64), tracked with carry into hi */
  uint128_t mid  = p01 + p10;
  int       mc   = (mid < p01);          /* carry from mid addition */
  uint128_t lo   = p00 + (mid << 64);
  int       lc   = (lo < p00);           /* carry from lo addition */

  *lo_out = lo;
  *hi_out = p11 + (mid >> 64) + ((uint128_t)mc << 64) + lc;
}

/* High 128 bits of the full 256-bit product a*b (lo half discarded). */
static INLINE uint128_t mulhi128(uint128_t a, uint128_t b) {
  uint64_t a0 = (uint64_t)a,  a1 = (uint64_t)(a >> 64);
  uint64_t b0 = (uint64_t)b,  b1 = (uint64_t)(b >> 64);

  uint128_t p00 = (uint128_t)a0 * b0;
  uint128_t p01 = (uint128_t)a0 * b1;
  uint128_t p10 = (uint128_t)a1 * b0;
  uint128_t p11 = (uint128_t)a1 * b1;

  uint128_t mid = p01 + p10;
  int       mc  = (mid < p01);
  uint128_t lo  = p00 + (mid << 64);     /* only needed for carry lc */
  int       lc  = (lo < p00);

  return p11 + (mid >> 64) + ((uint128_t)mc << 64) + lc;
}

/* Montgomery REDC: returns a*b*R^{-1} mod n.
 * Inputs must already be in [0, n). */
static INLINE uint128_t mont_mulmod128(uint128_t a, uint128_t b,
                                       const mont128_t *ctx) {
  uint128_t lo, hi;
  mul256(a, b, &lo, &hi);

  /* m = lo * ninv  (only low 128 bits needed) */
  uint128_t m = lo128(lo, ctx->ninv);

  /* m·n ≡ -lo (mod R), so lo + mn_lo = R when lo≠0, 0 otherwise.
   * carry into the high half is simply (lo != 0); no need to compute mn_lo. */
  int carry = (lo != 0);
  uint128_t mn_hi = mulhi128(m, ctx->n);

  /* hi + mn_hi + carry may exceed 2^128 when n > 2^127 (true result still
   * in [0, 2n), so one conditional subtract suffices).  Track overflow bit. */
  uint128_t t = hi + mn_hi;
  int ov = (t < hi);
  t += (uint128_t)carry;
  ov += (t < (uint128_t)carry);

  /* When ov=1 the true value is t + 2^128 >= n; uint128_t wraparound in
   * t - ctx->n correctly gives (t + 2^128 - n) mod 2^128. */
  if (ov || t >= ctx->n) return t - ctx->n;
  return t;
}

#define mont_sqrmod128(a, ctx)  mont_mulmod128(a, a, ctx)

/* Convert a (normal) value into Montgomery form: a*R mod n. */
static INLINE uint128_t mont_enter128(uint128_t a, const mont128_t *ctx) {
  return mont_mulmod128(a % ctx->n, ctx->r2, ctx);
}

/* Convert from Montgomery form back to normal: â*R^{-1} mod n. */
static INLINE uint128_t mont_exit128(uint128_t a, const mont128_t *ctx) {
  return mont_mulmod128(a, 1, ctx);
}

/* Initialise a Montgomery context for modulus n (must be odd). */
static void mont_setup128(mont128_t *ctx, uint128_t n) {
  ctx->n = n;

  /* ninv = n^{-1} mod 2^128 via Hensel lifting.
   * Each step doubles the number of correct bits: x = x*(2 - n*x). */
  uint128_t x = (3*n)^2;                 /*   4 bits (n odd) */
  x *= 2 - n * x;                        /*   8 bits */
  x *= 2 - n * x;                        /*  16 bits */
  x *= 2 - n * x;                        /*  32 bits */
  x *= 2 - n * x;                        /*  64 bits */
  x *= 2 - n * x;                        /* 128 bits */
  /* REDC add-form needs m = T*(-n^{-1}) mod R; store the negation */
  ctx->ninv = -x;

  /* r2 = R^2 mod n = (2^128)^2 mod n = 2^256 mod n. */
  ctx->r2 = pow2mod128(256,n);
}

uint128_t muladdmod128_s(uint128_t a, uint128_t b, uint128_t c,
                         uint128_t n, int sub)
{
  uint128_t r;

  if (a == 0 || b == 0) r = 0;
  else if (a == 1)      r = b;
  else if (b == 1)      r = a;
  else                  r = mulmod128(a, b, n);

  return sub ? submod128(r, c, n) : addmod128(r, c, n);
}

/*****************************************************************************
 * Montgomery powmod — exponent k is a plain UV, a/return in Montgomery form.
 *****************************************************************************/

static uint128_t mont_powmod128(uint128_t a, UV k, const mont128_t *ctx) {
  uint128_t r = mont_enter128(1, ctx);
  while (k > 0) {
    if (k & 1) r = mont_mulmod128(r, a, ctx);
    k >>= 1;
    if (k) a = mont_sqrmod128(a, ctx);
  }
  return r;
}

static uint128_t mont_powmod128_u128(uint128_t a, uint128_t k,
                                     const mont128_t *ctx) {
  uint128_t r = mont_enter128(1, ctx);
  while (k > 0) {
    if (k & 1) r = mont_mulmod128(r, a, ctx);
    k >>= 1;
    if (k) a = mont_sqrmod128(a, ctx);
  }
  return r;
}


/*****************************************************************************
 * Primality — BPSW (Miller-Rabin base 2 + extra-strong Lucas)
 *****************************************************************************/

/* We could add a base-2 M-R function, that avoided computing R2, but it is
 * only a small gain for the added code.  Can revisit later if desired. */

/* Returns 1 if n is a strong pseudoprime base `base`, 0 otherwise.
 * Assumes n > 2, n odd. */
static bool miller_rabin128(uint128_t n, uint128_t base) {
  uint128_t d = n - 1;
  mont128_t ctx;
  int s = 0;

  while (!(d & 1)) { d >>= 1; s++; }

  mont_setup128(&ctx, n);
  uint128_t mont1 = mont_enter128(1, &ctx);
  uint128_t montbase = (base == 2) ? addmod128(mont1, mont1, n)
                                   : mont_enter128(base % n, &ctx);
  uint128_t x = mont_powmod128_u128(montbase, d, &ctx);
  uint128_t montm1 = n - mont1;
  if (x == mont1 || x == montm1) return 1;
  while (--s > 0) {
    x = mont_sqrmod128(x, &ctx);
    if (x == montm1) return 1;
  }
  return 0;
}

/* Kronecker symbol (a/b), with each argument given as sign and magnitude.
 * Keeping the signs separate permits the full unsigned 128-bit range. */
int kronecker128(uint128_t a, int asign, uint128_t b, int bsign)
{
  int s = 1;

  if (b == 0) return a == 1;

  /* (a/-b) differs from (a/b) only when a is negative. */
  if (asign < 0 && bsign < 0) s = -s;

  if (!(b & 1)) {
    int odd_power = 0;
    if (!(a & 1)) return 0;
    do {
      b >>= 1;
      odd_power = !odd_power;
    } while (!(b & 1));
    if (odd_power && ((a & 7) == 3 || (a & 7) == 5)) s = -s;
  }

  /* The denominator is now positive and odd. */
  if (asign < 0 && (b & 3) == 3) s = -s;

  while (a != 0) {
    int odd_power = 0;
    while (!(a & 1)) {
      a >>= 1;
      odd_power = !odd_power;
    }
    if (odd_power && ((b & 7) == 3 || (b & 7) == 5)) s = -s;
    if ((a & 3) == 3 && (b & 3) == 3) s = -s;
    {
      uint128_t rem = b % a;
      b = a;
      a = rem;
    }
  }
  return (b == 1) ? s : 0;
}

/* Jacobi symbol (a/n), n odd > 0.  a may be any signed 128-bit value. */
static int jacobi128(int128_t a_in, uint128_t n)
{
  uint128_t a;
  int asign = 1;

  if (a_in < 0) {
    /* This form is safe even if a_in is the minimum signed value. */
    a = (uint128_t)(-(a_in + 1)) + 1;
    asign = -1;
  } else {
    a = (uint128_t)a_in;
  }
  return kronecker128(a, asign, n, 1);
}

/* Lucas sequence mod n in Montgomery form.  Given (D, P=1, Q), index k,
 * computes (U_k, V_k, Q^k), all returned in Montgomery form. */
static void lucas_seq128_mont(uint128_t *U, uint128_t *V, uint128_t *Qk,
                              int128_t D_s, int64_t Q_s,
                              uint128_t k, const mont128_t *ctx) {
  uint128_t n = ctx->n;
  uint128_t D = (D_s >= 0) ? (uint128_t)D_s % n
                           : n - (uint128_t)(-D_s) % n;
  uint128_t Q = (Q_s >= 0) ? (uint128_t)Q_s % n
                           : n - (uint128_t)(-Q_s) % n;
  uint128_t mont1 = mont_enter128(1, ctx);
  uint128_t mont2 = addmod128(mont1, mont1, n);

  D = mont_enter128(D, ctx);
  Q = mont_enter128(Q, ctx);

  uint128_t u = mont1, v = mont1, q = Q;
  uint128_t u2, v2, q2;
  uint128_t bit;

  if (k == 0) { *U = 0; *V = mont2; *Qk = mont1; return; }

  bit = (uint128_t)1 << 126;
  while (bit > k) bit >>= 1;
  bit >>= 1;

  while (bit > 0) {
    /* Double: U_{2m}, V_{2m}, Q^{2m} */
    u2 = mont_mulmod128(u, v, ctx);
    v2 = submod128(mont_sqrmod128(v, ctx), addmod128(q, q, n), n);
    q2 = mont_sqrmod128(q, ctx);

    if (k & bit) {
      /* Advance by 1, with P=1.  half_mod128 preserves Montgomery form. */
      uint128_t pu2_v2 = addmod128(u2, v2, n);
      uint128_t du2_v2 = addmod128(mont_mulmod128(D, u2, ctx), v2, n);
      u = half_mod128(pu2_v2, n);
      v = half_mod128(du2_v2, n);
      q = mont_mulmod128(q2, Q, ctx);
    } else {
      u = u2; v = v2; q = q2;
    }
    bit >>= 1;
  }

  *U = u; *V = v; *Qk = q;
}

/* Strong Lucas probable prime test with Selfridge parameters.
 * n must be odd, > 2, not a perfect square.
 * D is chosen from sequence 5,-7,9,-11,13,... until jacobi(D,n)=-1.
 * P=1, Q=(1-D)/4. */
MAYBE_UNUSED static bool is_strong_lucas128(uint128_t n) {
  int64_t D;
  int64_t sign = 1;
  int r, abs_D;

  for (abs_D = 5; ; abs_D += 2, sign = -sign) {
    D = sign * abs_D;
    int j = jacobi128((int128_t)D, n);
    if (j == -1) break;
    if (j == 0)  return 0; /* gcd(D,n) > 1 => composite */
    /* Selfridge parameter search cannot terminate for a coprime square. */
    if (abs_D == 21 && is_perfect_square128(n)) return 0;
  }

  int64_t Q_s = (1 - D) / 4;
  uint128_t d = n + 1;
  int s = 0;
  while (!(d & 1)) { d >>= 1; s++; }

  uint128_t U, V, Qk;
  mont128_t ctx;
  mont_setup128(&ctx, n);
  lucas_seq128_mont(&U, &V, &Qk, D, Q_s, d, &ctx);

  if (U == 0) return 1;
  for (r = 0; r < s; r++) {
    if (V == 0) return 1;
    /* V_{2m} = V_m^2 - 2*Q^m  (mod n) */
    V = submod128(mont_sqrmod128(V, &ctx), addmod128(Qk, Qk, n), n);
    Qk = mont_sqrmod128(Qk, &ctx);
  }
  return 0;
}

/* Select Baillie's extra-strong Lucas parameters.  Q=1 and P is the first
 * value in 3,4,5,... for which jacobi(P^2-4,n) == -1. */
static bool select_extra_strong_parameters128(uint32_t *Pout, uint128_t n) {
  uint128_t g;
  uint64_t D;
  uint32_t P = 3;
  int j;

  while (1) {
    D = (uint64_t)P * P - 4;
    j = jacobi128((int128_t)D, n);
    if (j == 0) {
      g = gcd128((uint128_t)D, n);
      if (g != 1 && g != n) return 0;
    }
    if (j == -1) break;
    if (P == 23 && is_perfect_square128(n)) return 0;
    if (++P > 65535)
      croak("lucas_extrastrong_params: P exceeded 65535");
  }
  *Pout = P;
  return 1;
}

/* Full extra-strong Lucas probable-prime test with Baillie's parameters.
 * n must be odd, greater than 2, and less than the maximum uint128_t. */
static bool is_extra_strong_lucas128(uint128_t n) {
  mont128_t ctx;
  uint128_t U, V, d, bit, t2, unext, vnext;
  uint128_t mont1, mont2, montP, montD;
  uint64_t D;
  uint32_t P;
  int s = 0;

  if (!select_extra_strong_parameters128(&P, n)) return 0;
  D = (uint64_t)P * P - 4;

  d = n + 1;
  while (!(d & 1)) { d >>= 1; s++; }

  mont_setup128(&ctx, n);
  mont1 = mont_enter128(1, &ctx);
  mont2 = addmod128(mont1, mont1, n);
  montP = mont_enter128(P, &ctx);
  montD = mont_enter128(D, &ctx);
  U = mont1;
  V = montP;

  bit = (uint128_t)1 << 126;
  while (bit > d) bit >>= 1;
  bit >>= 1;
  while (bit > 0) {
    U = mont_mulmod128(U, V, &ctx);
    V = submod128(mont_sqrmod128(V, &ctx), mont2, n);
    if (d & bit) {
      t2 = mont_mulmod128(U, montD, &ctx);
      unext = addmod128(mont_mulmod128(U, montP, &ctx), V, n);
      vnext = addmod128(mont_mulmod128(V, montP, &ctx), t2, n);
      U = half_mod128(unext, n);
      V = half_mod128(vnext, n);
    }
    bit >>= 1;
  }

  if (U == 0 && (V == mont2 || V == n-mont2)) return 1;
  s--;
  while (s--) {
    if (V == 0) return 1;
    if (s) V = submod128(mont_sqrmod128(V, &ctx), mont2, n);
  }
  return 0;
}

/* BPSW primality test, no trial division.  Assumes n > 2011 and odd.
 * Called from the factoring loop where small factors are already removed. */
bool is_bpsw128(uint128_t n) {
  if (n < 7) return (n==2 || n==3 || n==5);
  if (!(n & 1) || n == (uint128_t)-1) return 0;
  if (!miller_rabin128(n, 2)) return 0;
  if (!is_extra_strong_lucas128(n)) return 0;
  return 1;
}

/* Returns 1 if n is (probably) prime, 0 if composite.
 * Uses trial division up to 2011, then BPSW. */
bool is_prime128(uint128_t n) {
  int i;

  if (n < 7) return (n==2 || n==3 || n==5);
  if (!(n & 1)) return 0;

  /* Trial division using primes_small[] */
  for (i = 2; i < NPRIMES_SMALL; i++) {
    uint64_t p = primes_small[i];
    if ((uint128_t)p * p > n) return 1;
    if (n % p == 0) return 0;
  }

  return is_bpsw128(n);
}

/* Integer square root of a 128-bit value: returns floor(sqrt(n)). */
static uint64_t isqrt128(uint128_t n) {
  if (n == 0) return 0;
  /* Initial estimate: r = 2^ceil(bits/2), an overestimate of sqrt(n) */
  uint128_t r = 1, tmp = n;
  while (tmp > 3) { tmp >>= 2; r <<= 1; }
  r <<= 1;   /* ensure we start from above so Newton descends to floor */
  /* Newton's method descends from above: r = (r + n/r) / 2 */
  while (1) {
    uint128_t rn = (r + n / r) >> 1;
    if (rn >= r) break;
    r = rn;
  }
  /* r may still be one too high; adjust down */
  while (r * r > n) r--;
  return (uint64_t)r;
}

bool is_perfect_square128_ret(uint128_t n, uint64_t *root) {
  uint64_t r;
  uint32_t m;

  /* Fast filters reject 95.0% of non-squares. */
  if ((UINT64_C(1) << ((uint64_t)n & 63)) &
      UINT64_C(0xfdfdfdedfdfcfdec))
    return 0;
  m = (uint32_t)(n % 45);
  if ((UINT64_C(1) << m) & UINT64_C(0xfffffeeb7df6f9ec))
    return 0;

  r = isqrt128(n);
  if (root != 0) *root = r;
  return (uint128_t)r * r == n;
}

/*****************************************************************************
 * Pollard P-1 factoring — Stage 1 + Stage 2 using Montgomery arithmetic.
 *
 * Stage 1: a = 2^E mod n  where E = ∏ q^e, q prime, q^e ≤ B1.
 *          Check gcd(a−1, n).  GCD batched every 32 primes; backtrack on
 *          gcd == n to isolate the exact factor.
 *
 * Stage 2: standard prime continuation.  bm = a after stage 1.
 *          For each prime p ∈ (B1, B2]: a = a · bm^(p−p_prev) mod n.
 *          Accumulate b = ∏(a−1) mod n; GCD every 64 primes.
 *          Small prime-gap powers bm^{2k} (k=1..111) are cached lazily.
 *****************************************************************************/

static uint128_t pminus1_128(uint128_t n, uint64_t B1_in, uint64_t B2_in) {
  mont128_t ctx;
  UV B1, B2;

  if (B1_in > (uint64_t)UV_MAX || B2_in > (uint64_t)UV_MAX)
    croak("internal: pminus1_128 bounds exceed UV_MAX");

  B1 = (UV)B1_in;
  B2 = (UV)B2_in;

  if (B1 < 7) return 0;
  mont_setup128(&ctx, n);

  /*--- Stage 1 -----------------------------------------------------------*/
  uint128_t a     = mont_enter128(2, &ctx);
  uint128_t savea = a;
  UV q = 2, saveq = 1;
  UV sqrtB1 = isqrt(B1);
  UV j = 15;   /* checkpoint counter, start offset like GMP */

  START_DO_FOR_EACH_PRIME(2,B1) {
    UV k = p;
    if (p <= sqrtB1) { UV km = B1/p; while (k <= km) k *= p; }
    a = mont_powmod128(a, k, &ctx);
    if ((j++ % 32) == 0) {
      uint128_t an = mont_exit128(a, &ctx);
      uint128_t g  = gcd128(an > 0 ? an - 1 : n - 1, n);
      if (g == n)  RETURN_FROM_EACH_PRIME(goto stage1_backtrack);
      if (g > 1)   RETURN_FROM_EACH_PRIME(return g);
      saveq = p;
      savea = a;
    }
  } END_DO_FOR_EACH_PRIME

  /* Final stage-1 GCD */
  {
    uint128_t an = mont_exit128(a, &ctx);
    uint128_t g  = gcd128(an > 0 ? an - 1 : n - 1, n);
    if (g == n)  goto stage1_backtrack;
    if (g > 1)   return g;
  }
  goto stage2;

stage1_backtrack:
  /* savea is the state after saveq; resume at the following prime. */
  a = savea;
  for (q = next_prime(saveq); q <= B1; q = next_prime(q)) {
    UV k = q;
    if (q <= sqrtB1) { UV km = B1/q; while (k <= km) k *= q; }
    a = mont_powmod128(a, k, &ctx);
    uint128_t an = mont_exit128(a, &ctx);
    uint128_t g  = gcd128(an > 0 ? an - 1 : n - 1, n);
    if (g == n)  return 0;
    if (g > 1)   return g;
  }
  return 0;

stage2:
  if (B2 <= B1) return 0;

  /*--- Stage 2 -----------------------------------------------------------*/
  /* bm = a (end-of-stage-1 value).  We precompute bm^2, bm^4, …, bm^44
   * eagerly, and lazily cache up to bm^222.  The index is:
   *   qdiff = (prime_gap)/2 − 1,  bm_pow[qdiff] = bm^(prime_gap).       */
  uint128_t bm      = a;
  uint128_t bm_sq2  = mont_sqrmod128(bm, &ctx);      /* bm^2 */
  uint128_t bm_pow[111];
  int       bm_init[111];
  int       ii;
  for (ii = 0; ii < 111; ii++) bm_init[ii] = 0;

  bm_pow[0]  = bm_sq2;  bm_init[0]  = 1;         /* bm^2  */
  for (ii = 1; ii < 22; ii++) {                  /* bm^4 … bm^44 */
    bm_pow[ii] = mont_mulmod128(bm_pow[ii-1], bm_sq2, &ctx);
    bm_init[ii] = 1;
  }

  q = next_prime(B1);  /* q is now the first prime > B1 */
  /* Advance a to the first stage-2 prime (= q, first prime > B1) */
  a = mont_powmod128(a, q, &ctx);

  uint128_t R_modn = mont_enter128(1, &ctx);   /* Mont. form of 1 */
  uint128_t b = R_modn;                        /* accumulator */
  j = 31;

  START_DO_FOR_EACH_PRIME(q,B2) {
    uint128_t step;
    UV gap   = p - q;  /* q is the previous prime */
    UV pdiff = gap / 2 - 1;
    if (pdiff < 111) {
      if (!bm_init[pdiff]) {
        /* lazily compute bm^gap = bm^(2*(qdiff+1)) */
        bm_pow[pdiff] = mont_powmod128(bm, gap, &ctx);
        bm_init[pdiff] = 1;
      }
      step = bm_pow[pdiff];
    } else {
      /* rare large gap — compute on the fly */
      step = mont_powmod128(bm, gap, &ctx);
    }
    a = mont_mulmod128(a, step, &ctx);
    /* (a − 1) in Montgomery form: a_mont − R_modn */
    uint128_t am1 = submod128(a, R_modn, n);
    b = mont_mulmod128(b, am1, &ctx);

    if ((j++ % 64) == 0) {
      uint128_t g = gcd128(mont_exit128(b, &ctx), n);
      if (g > 1 && g < n)  RETURN_FROM_EACH_PRIME(return g);
      if (g == n)          break;
    }
    q = p;
  } END_DO_FOR_EACH_PRIME
  /* Final stage-2 GCD */
  {
    uint128_t g = gcd128(mont_exit128(b, &ctx), n);
    if (g > 1 && g < n)  return g;
  }
  return 0;
}

/*****************************************************************************
 * Tiny ECM — elliptic curve factoring.
 *
 * Targets factors in the 40–55 bit range when earlier methods fail.
 * Uses Suyama's parameterization, the Montgomery ladder, and the existing
 * mont_mulmod128 infrastructure.
 *****************************************************************************/

#if 0  /* Set to 1 for tuning ECM ladder, assumes local variables. */
#define ECM128_REPORT(stage, g) \
  do { \
    if ((g) > 1 && (g) < n) \
      printf("ECM128 %luk/%luk " stage " sigma %u found factor %s of %s\n", \
             (unsigned long)B1_in/1000, (unsigned long)B2_in/1000, \
             sigma, u128_str(g), u128_str(n)); \
  } while (0)
#else
#define ECM128_REPORT(stage, g) do { } while (0)
#endif


/* Projective (X:Z) point on a Montgomery curve.  Values in Montgomery form. */
typedef struct { uint128_t X, Z; } ecpt128_t;

/* Point doubling: R = 2P on By²=x³+Ax²+x, A24 = (A+2)/4. */
static INLINE void ecm_double128(ecpt128_t *R, const ecpt128_t *P,
                                 uint128_t A24, const mont128_t *ctx) {
  uint128_t u = mont_sqrmod128(submod128(P->X, P->Z, ctx->n), ctx); /* (X-Z)² */
  uint128_t v = mont_sqrmod128(addmod128(P->X, P->Z, ctx->n), ctx); /* (X+Z)² */
  uint128_t w = submod128(v, u, ctx->n);                             /* 4XZ   */
  R->X = mont_mulmod128(u, v, ctx);
  R->Z = mont_mulmod128(w, addmod128(u, mont_mulmod128(A24, w, ctx), ctx->n), ctx);
}

/* Differential addition: R = P+Q given P-Q = Pm.  R may alias P or Q. */
static INLINE void ecm_dadd128(ecpt128_t *R,
                               const ecpt128_t *P, const ecpt128_t *Q,
                               const ecpt128_t *Pm, const mont128_t *ctx) {
  uint128_t u = mont_mulmod128(submod128(P->X, P->Z, ctx->n),
                               addmod128(Q->X, Q->Z, ctx->n), ctx);
  uint128_t v = mont_mulmod128(addmod128(P->X, P->Z, ctx->n),
                               submod128(Q->X, Q->Z, ctx->n), ctx);
  uint128_t s = addmod128(u, v, ctx->n);
  uint128_t d = submod128(u, v, ctx->n);
  R->X = mont_mulmod128(mont_sqrmod128(s, ctx), Pm->Z, ctx);
  R->Z = mont_mulmod128(mont_sqrmod128(d, ctx), Pm->X, ctx);
}

/* Montgomery-ladder scalar multiply: R = k*P.
 * R1-R0 = P is maintained, so the original P is always the Pm arg.
 * Safe to call with R aliasing P: R0 is a local copy made before any writes. */
static void ecm_mul128(ecpt128_t *R, const ecpt128_t *P, UV k,
                       uint128_t A24, const mont128_t *ctx) {
  if (k == 1) { *R = *P; return; }
  ecpt128_t R0 = *P, R1;
  ecm_double128(&R1, &R0, A24, ctx);
  if (k == 2) { *R = R1; return; }
  UV bit = (UV)1 << (8*sizeof(UV) - 1);
  while (!(bit & k)) bit >>= 1;
  for (bit >>= 1; bit; bit >>= 1) {
    if (k & bit) {
      ecm_dadd128(&R0, &R0, &R1, P, ctx);
      ecm_double128(&R1, &R1, A24, ctx);
    } else {
      ecm_dadd128(&R1, &R0, &R1, P, ctx);
      ecm_double128(&R0, &R0, A24, ctx);
    }
  }
  *R = R0;
}

/* Fixed σ values for Suyama's parameterization.
 * All σ ≥ 11 so u = σ²-5 > v = 4σ (non-degenerate).
 * For this table, all prime factors of u are <= 2011, hence removed by
 * trial division.  The denominator 16u³v also contains sigma itself;
 * if that divides n, curve setup can expose it via gcd(den,n). */
static const uint16_t ecm_sigmas[] = {
   11,   13,   17,   19,   23,   29,   31,   37,   41,   43,
   47,   53,   59,   61,   67,   71,   73,   79,   83,   89,
  103,  127,  139,  149,  151,  157,  163,  181,  191,  197,
  199,  211,  223,  227,  233,  239,  257,  271,  293,  313, /*  40 */
  331,  337,  347,  359,  379,  389,  397,  401,  409,  421,
  443,  449,  457,  479,  487,  509,  521,  523,  547,  557,
  587,  641,  653,  659,  673,  677,  683,  691,  719,  727,
  739,  751,  769,  797,  809,  853,  919,  929,  941,  997,
 1049, 1051, 1063, 1091, 1093, 1109, 1117, 1129, 1153, 1201,
 1217, 1229, 1283, 1327, 1361, 1381, 1427, 1447, 1459, 1471, /* 100 */
 1481, 1489, 1543, 1549, 1571, 1621, 1709, 1723, 1753, 1759,
 1801, 1811, 1867, 1987, 2039, 2099, 2113, 2131, 2251, 2309,
 2347, 2381, 2399, 2447, 2473, 2551, 2557, 2663, 2677, 2689,
 2713, 2719, 2749, 2857, 2879, 2887, 2939, 3001, 3061, 3067,
 3121, 3137, 3187, 3251, 3259, 3271, 3307, 3359, 3371, 3373,
 3467, 3593, 3607, 3623, 3643, 3709, 3733, 3793, 3851, 3923, /* 160 */
 3989, 4019, 4049, 4129, 4231, 4253, 4283, 4339, 4349, 4441,
 4523, 4649, 4787, 4987, 4999, 5171, 5237, 5273, 5297, 5333,
 5387, 5471, 5479, 5647, 5749, 5791, 6101, 6163, 6257, 6299,
 6337, 6451, 6491, 6659, 6793, 6823, 6967, 7013, 7229, 7253, /* 200 */
 7333, 7369, 7477, 7621, 7793, 7817, 8059, 8167, 8209, 8263,
 8311, 8377, 8573, 8641, 8741, 8837, 8863, 8963, 9001, 9151,
 9203, 9433, 9697, 9743, 9781, 9883,10007,10069,10099,10139,
10163,10193,10267,10429,10457,10487,10691,10837,10949,11087,
11243,11321,11411,11681,11813,11903,12011,12263,12277,12401,
12409,12437,12479,12569,12619,12739,12911,13331,13367,13537,
13721,13789,13841,13873,14051,14149,14221,14419,14431,14827,
14887,15077,15289,15467,15511,15649,15773,15797,15859,15901,
16057,16141,16217,16529,16547,16553,16619,17299,17393,17419,
17449,17737,17921,18049,18223,19073,19183,19477,20021,20323,
20347,20759,20929,21023,21157,21587,21611,21613,21673,21751,
21799,21821,22109,22469,22651,22943,23327,23459,23567,23767,
23911,23957,24001,24197,24281,24407,24799,24851,25147,25183,
25469,25679,25703,26561,26683,26701,26821,27073,27191,27271,
27277,27427,27487,27539,27617,27647,27673,27749,27983,28319,
28789,28843,29017,29123,29209,29669,29803,29921,30323,30809,
30851,30911,30983,31397,31541,31963,32369,32561,32771,32969,
33029,33083,33487,33637,33757,34057,34381,34513,34613,34807,
35083,35171,35311,35381,36013,36251,36493,36529,36551,36913,
36919,37363,37517,37699,37907,38047,38177,38273,38749,38903  /* 400 */
};
#define NECM128_SIGMAS ((int)(sizeof(ecm_sigmas)/sizeof(ecm_sigmas[0])))

/* Batch convert projective Montgomery points to affine x-coordinates, still
 * in Montgomery form.  Returns 1 on success.  On failure, *fout is either a
 * non-trivial factor or 0 when this curve should be abandoned. */
static int ecm_batch_normalize_x128(uint128_t *xout, uint128_t *fout,
                                    const ecpt128_t *P, UV npoints,
                                    const mont128_t *ctx) {
  uint128_t n = ctx->n, R = mont_enter128(1, ctx);
  uint128_t acc = R, inv, g;
  uint128_t *prefix;
  UV i;

  if (npoints == 0) { *fout = 0; return 1; }

  prefix = (uint128_t*)mpu_aligned_alloc(npoints, sizeof(uint128_t),
                                         sizeof(uint128_t));

  for (i = 0; i < npoints; i++) {
    acc = mont_mulmod128(acc, P[i].Z, ctx);
    prefix[i] = acc;
  }

  g = gcd128(acc, n);
  if (g > 1) {
    *fout = (g < n) ? g : 0;
    /* Different coordinates can expose different factors whose product has
     * gcd n.  Check them individually only on this uncommon failure path. */
    if (g == n) {
      for (i = 0; i < npoints; i++) {
        g = gcd128(P[i].Z, n);
        if (g > 1 && g < n) { *fout = g; break; }
      }
    }
    mpu_aligned_free(prefix);
    return 0;
  }

  inv = modinv128(mont_exit128(acc, ctx), n);
  if (inv == 0) {
    *fout = 0;
    mpu_aligned_free(prefix);
    return 0;
  }
  inv = mont_enter128(inv, ctx);

  for (i = npoints; i > 0; i--) {
    UV j = i - 1;
    uint128_t prev = (j == 0) ? R : prefix[j-1];
    uint128_t zinv = mont_mulmod128(inv, prev, ctx);
    inv = mont_mulmod128(inv, P[j].Z, ctx);
    xout[j] = mont_mulmod128(P[j].X, zinv, ctx);
  }

  *fout = 0;
  mpu_aligned_free(prefix);
  return 1;
}

/* Brent-Suyama style ECM stage 2.  Q is the stage-1 output point. */
static uint128_t tinyecm128_stage2(const ecpt128_t *Q, uint128_t A24,
                                   const mont128_t *ctx,
                                   uint64_t B1_in, uint64_t B2_in) {
  uint128_t n = ctx->n, f, *nqx, *Sx;
  ecpt128_t *nq, *S;
  UV B1, B2, D, twoD, m, mend, nwindows, w, i;
  UV const B2max = UV_MAX - (BITS_PER_WORD == 64 ? 3037000500U : 46340U);
  uint128_t gprod;

  if (B1_in > (uint64_t)UV_MAX || B2_in > (uint64_t)B2max)
    croak("internal: tinyecm128 stage 2 bounds exceed UV_MAX");

  B1 = (UV)B1_in;
  B2 = (UV)B2_in;

  if (B2 <= B1 || B2 < 2) return 0;
  D = isqrt(B2 >> 1);
  if (D & 1) D++;
  twoD = 2*D;

  mend = B2 + D;
  nwindows = 1 + (mend - 2) / twoD;

  nq = (ecpt128_t*)mpu_aligned_alloc(twoD+1, sizeof(ecpt128_t),
                                     sizeof(uint128_t));
  nqx = (uint128_t*)mpu_aligned_alloc(D+1, sizeof(uint128_t),
                                      sizeof(uint128_t));

  nq[1] = *Q;
  for (i = 2; i <= twoD; i++) {
    if (i & 1) {
      ecm_dadd128(&nq[i], &nq[(i-1)/2], &nq[(i+1)/2], Q, ctx);
    } else {
      ecm_double128(&nq[i], &nq[i/2], A24, ctx);
    }
  }

  nqx[0] = 0;
  if (!ecm_batch_normalize_x128(nqx+1, &f, nq+1, D, ctx)) {
    mpu_aligned_free(nqx);
    mpu_aligned_free(nq);
    return f;
  }

  S = (ecpt128_t*)mpu_aligned_alloc(nwindows, sizeof(ecpt128_t),
                                    sizeof(uint128_t));
  Sx = (uint128_t*)mpu_aligned_alloc(nwindows, sizeof(uint128_t),
                                     sizeof(uint128_t));

  S[0] = *Q;
  {
    ecpt128_t Xm = nq[twoD-1];
    for (w = 1; w < nwindows; w++) {
      ecpt128_t oldS = S[w-1];
      ecm_dadd128(&S[w], &nq[twoD], &S[w-1], &Xm, ctx);
      Xm = oldS;
    }
  }

  if (!ecm_batch_normalize_x128(Sx, &f, S, nwindows, ctx)) {
    mpu_aligned_free(Sx);
    mpu_aligned_free(S);
    mpu_aligned_free(nqx);
    mpu_aligned_free(nq);
    return f;
  }

  gprod = mont_enter128(1, ctx);
  m = 1;
  for (w = 0; w < nwindows; w++) {
    UV hi = (m > B2 || B2 - m < D) ? B2 : m + D;
    UV lo = (m > D) ? m - D : 0;

    if (hi > B1) {
      if (lo <= B1) lo = B1 + 1;
      if (lo < 2) lo = 2;
      if (lo <= hi) {
        START_DO_FOR_EACH_PRIME(lo,hi) {
          UV idx;
          if (p < m) {
            idx = m - p;
          } else if (p > m) {
            if (m <= UV_MAX/2) {
              UV mm = 2*m;
              if (p <= mm) {
                UV mirror = mm - p;
                if (mirror > B1 && mirror >= lo && is_prime(mirror))
                  continue;
              }
            }
            idx = p - m;
          } else {
            continue;
          }
          if (idx <= D) {
            uint128_t diff = submod128(Sx[w], nqx[idx], n);
            gprod = mont_mulmod128(gprod, diff, ctx);
          }
        } END_DO_FOR_EACH_PRIME
        f = gcd128(mont_exit128(gprod, ctx), n);
        if (f > 1) {
          mpu_aligned_free(Sx);
          mpu_aligned_free(S);
          mpu_aligned_free(nqx);
          mpu_aligned_free(nq);
          return (f < n) ? f : 0;
        }
      }
    }
    if (w + 1 < nwindows) m += twoD;
  }

  f = gcd128(mont_exit128(gprod, ctx), n);
  mpu_aligned_free(Sx);
  mpu_aligned_free(S);
  mpu_aligned_free(nqx);
  mpu_aligned_free(nq);
  return (f > 1 && f < n) ? f : 0;
}

/* Returns a non-trivial factor of n, or 0.
 * ncurves curves, Suyama parameterization, stage-1 bound B1, optional B2. */
static uint128_t tinyecm128(uint128_t n, uint64_t B1_in, uint64_t B2_in,
                            int ncurves, int sigma_offset) {
  mont128_t ctx;
  UV B1, sqrtB1;
  int ci;

  if (B1_in > (uint64_t)UV_MAX)
    croak("internal: tinyecm128 B1 exceeds UV_MAX");

  B1 = (UV)B1_in;
  sqrtB1 = isqrt(B1);

  mont_setup128(&ctx, n);

  for (ci = 0; ci < ncurves && sigma_offset+ci < NECM128_SIGMAS; ci++) {
    uint32_t sigma = ecm_sigmas[sigma_offset+ci];

    /* Suyama parameterization. */
    uint64_t ui    = (uint64_t)sigma * sigma - 5; /* u = σ²-5  */
    uint64_t vi    = 4 * (uint64_t)sigma;         /* v = 4σ    */
    uint64_t umv   = ui - vi;                     /* u-v > 0 for σ ≥ 6 */
    uint64_t t3uv  = 3 * ui + vi;
    uint128_t u3   = (uint128_t)ui * ui * ui;     /* u³ */
    uint128_t v3   = (uint128_t)vi * vi * vi;     /* v³ */
    uint128_t umv3 = (uint128_t)umv * umv * umv;

    /* A24 = (v-u)³(3u+v) / (16u³v) = -(u-v)³(3u+v) / (16u³v)  mod n */
    uint128_t abs_num = umv3 * t3uv;
    uint128_t abs_num_r = abs_num % n;
    uint128_t num = abs_num_r ? n - abs_num_r : 0;

    uint128_t den     = 16 * u3 * vi;
    uint128_t den_r   = den % n;
    uint128_t den_inv = modinv128(den_r, n);
    if (den_inv == 0) {
      uint128_t g = gcd128(den_r, n);
      if (g > 1 && g < n) return g;
      continue;
    }

    uint128_t mnum  = mont_enter128(num,     &ctx);   /* num * R mod n */
    uint128_t mdinv = mont_enter128(den_inv, &ctx);   /* den_inv * R mod n */
    uint128_t mA24  = mont_mulmod128(mnum, mdinv, &ctx);

    ecpt128_t P    = { mont_enter128(u3, &ctx),
                       mont_enter128(v3, &ctx) };

    /* Stage 1: apply prime powers ≤ B1; check gcd(Z, n) every 64 primes. */
    uint32_t j = 0;
    START_DO_FOR_EACH_PRIME(2,B1) {
      UV k = p;
      if (p <= sqrtB1) { UV pm = B1 / p; while (k <= pm) k *= p; }
      /* A GMP-style PRAC multiply could replace the generic ladder here. */
      ecm_mul128(&P, &P, k, mA24, &ctx);
      if ((j++ % 64) == 0) {
        uint128_t g = gcd128(P.Z, n);
        ECM128_REPORT("S1", g);
        if (g > 1 && g < n) RETURN_FROM_EACH_PRIME(return g);
        if (g == n)         RETURN_FROM_EACH_PRIME(goto next_curve);
      }
    } END_DO_FOR_EACH_PRIME
    { uint128_t g = gcd128(P.Z, n);
      ECM128_REPORT("S1", g);
      if (g > 1 && g < n) return g; }

    if (B2_in > B1_in) {
      uint128_t g = tinyecm128_stage2(&P, mA24, &ctx, B1_in, B2_in);
      ECM128_REPORT("S2", g);
      if (g > 1 && g < n) return g;
    }

next_curve:;
  }
  return 0;
}

#if 0
/*****************************************************************************
 * Brent's rho factoring
 *****************************************************************************/

/* Returns a non-trivial factor of n, or 0 on failure.
 * c is the polynomial constant x^2+c; try several values of c on failure.
 * Uses Montgomery multiplication for speed. */
static uint128_t brent128(uint128_t n, uint128_t c, UV maxr) {
  mont128_t ctx;
  mont_setup128(&ctx, n);

  UV m = 128;
  uint128_t r = 1, p;

  /* Convert to Montgomery form */
  uint128_t cy = mont_enter128(2, &ctx);   /* y in Montgomery form */
  uint128_t cc = mont_enter128(c, &ctx);   /* c in Montgomery form */
  uint128_t cq = mont_enter128(1, &ctx);   /* accumulator, starts at 1 */
  uint128_t cx, cys;

#define MBRENT_F(v)  addmod128(mont_sqrmod128(v, &ctx), cc, n)

  cx = cy;
  do {
    UV i, k = 0;
    cx = cy;
    for (i = 0; i < r; i++)
      cy = MBRENT_F(cy);
    do {
      cys = cy;
      for (i = 0; i < (m < r - k ? m : r - k); i++) {
        cy = MBRENT_F(cy);
        uint128_t diff = (cx > cy) ? cx - cy : cy - cx;
        cq = mont_mulmod128(cq, diff, &ctx);
      }
      p = gcd128(cq, n);   /* gcd(q_mont,n)=gcd(q_value,n) since gcd(R,n)=1 */
      k += m;
    } while (k < r && p == 1);
    if (r > maxr) return 0;
    r <<= 1;
  } while (p == 1);

  if (p == n) {
    do {
      cys = MBRENT_F(cys);
      uint128_t diff = (cx > cys) ? cx - cys : cys - cx;
      p = gcd128(diff, n);
    } while (p == 1);
  }
#undef MBRENT_F

  return (p == n) ? 0 : p;
}
#endif

/*****************************************************************************
 * Result management
 *****************************************************************************/

/* Insert a 64-bit prime factor p into nf (sorted).  Increments exponent if
 * already present.  Caller must ensure nfactors < MPU_MAX_DFACTORS128. */
static void insert_factor128(factored128_t *nf, uint64_t p) {
  uint16_t i, j, nf_ = nf->nfactors;
  for (i = 0; i < nf_; i++) {
    if (nf->f[i] == p) { nf->e[i]++; return; }
    if (nf->f[i] >  p) break;
  }
  /* Shift right and insert */
  for (j = nf_; j > i; j--) {
    nf->f[j] = nf->f[j-1];
    nf->e[j] = nf->e[j-1];
  }
  nf->f[i] = p;
  nf->e[i] = 1;
  nf->nfactors++;
}

/*****************************************************************************
 * Public helper functions
 *****************************************************************************/

uint32_t factored128p_total_factors(const factored128_t *nf) {
  uint32_t i, nfacs = nf->flarge ? 1 : 0;
  for (i = 0; i < nf->nfactors; i++)
    nfacs += nf->e[i];
  return nfacs;
}
uint32_t factored128p_distinct_factors(const factored128_t *nf) {
  return nf->nfactors + (nf->flarge ? 1 : 0);
}
bool factored128p_is_square_free(const factored128_t *nf) {
  uint32_t i;
  for (i = 0; i < nf->nfactors; i++)
    if (nf->e[i] > 1)
      return 0;
  return 1;
}
signed char factored128p_moebius(const factored128_t *nf) {
  uint32_t i, nfactors = nf->flarge ? 1 : 0;
  for (i = 0; i < nf->nfactors; i++) {
    if (nf->e[i] > 1)
      return 0;
    nfactors++;
  }
  return nfactors % 2 ? -1 : 1;
}

signed char moebius128(uint128_t n) {
  factored128_t nf;
  uint32_t i;

  if (n <= (uint128_t)UV_MAX)
    return moebius((UV)n);

  for (i = 1; i < NPRIMES_SMALL && primes_small[i] <= 37; i++) {
    uint128_t p = primes_small[i];
    if ((n % (p*p)) == 0)
      return 0;
  }

  factorintp128(&nf, n);
  return factored128p_moebius(&nf);
}

bool is_semiprime128(uint128_t n) {
  factored128_t nf;
  uint64_t r;
  uint32_t i;

  if (n <= (uint128_t)UV_MAX)
    return is_semiprime((UV)n);

  if (!(n & 1)) return is_prime128(n >> 1);
  if (!(n % 3)) return is_prime128(n / 3);
  if (!(n % 5)) return is_prime128(n / 5);

  for (i = 4; i < NPRIMES_SMALL; i++) {
    uint128_t p = primes_small[i];
    if (p*p*p > n)
      break;
    if ((n % p) == 0)
      return is_prime128(n / p);
  }

  if (is_prime128(n))
    return 0;

  if (is_perfect_square128_ret(n, &r))
    return is_prime128(r);

  factorintp128(&nf, n);
  return factored128p_total_factors(&nf) == 2;
}

/*****************************************************************************
 * Public entry point
 *****************************************************************************/

void factorintp128(factored128_t *nf, uint128_t n) {
  int const verbose = _XS_get_verbose();
  int i;

  nf->n        = n;
  nf->nfactors = 0;
  nf->flarge   = 0;

  if (n < 4) {
    nf->f[0] = (uint64_t)n;
    nf->e[0] = 1;
    nf->nfactors = 1 - (n == 1);
    return;
  }

  /* Trial division using primes_small[] (primes 2..2011) */
  for (i = 1; i < NPRIMES_SMALL; i++) {
    uint64_t p = primes_small[i];
    if ((uint128_t)p * p > n) break;
    while (n % p == 0) {
      insert_factor128(nf, p);
      n /= p;
    }
  }
  if (n == 1) return;

  /* n is now > 2011; remaining cofactor stack */
  uint128_t stack[128];
  int top = 0;
  stack[top++] = n;

  while (top > 0) {
    uint128_t t = stack[--top];
    int show = verbose;
    /* printf("  factor128 working on %s\n", u128_str(t)); */

    /* Shortcut: fits in UV */
    if (t <= (uint128_t)UV_MAX) {
      UV uv = (UV)t;
      UV factors[MPU_MAX_FACTORS];
      int nf2 = factor(uv, factors);
      for (i = 0; i < nf2; i++)
        insert_factor128(nf, factors[i]);
      continue;
    }

    if (is_bpsw128(t)) {
      /* Store primes that fit in the 64-bit factor slots.  A larger prime
       * can only appear once, because the original input is at most 128 bits. */
      if (t <= (uint128_t)UINT64_MAX) insert_factor128(nf, (uint64_t)t);
      else                            nf->flarge = t;
      continue;
    }

    /* Composite: find a factor */
    uint128_t f = 0;


    /* P-1 stage 1+2: inexpensive first try */
    if (!f) f = pminus1_128(t, 2000, 10000);
    if (f && show) {show=0;printf("p-1/128 found factor %s\n",u128_str(f));}


    /* tinysiqs128 gets the factor quickly. */
    if (!f) f = tinysiqs128(t);
    if (f && show) {show=0;printf("tinysiqs128 found factor %s\n",u128_str(f));}


    /* tinyecm128 used to be our main factoring method for 65-128 bit inputs,
     * but now it's only used as a fallback. */

#define TINYECM_STEP(eb1,b2mult,curves) \
    if (!f) { \
      uint64_t B1_ = (uint64_t)(eb1) * 1000; \
      uint64_t B2_ = B1_ * (uint64_t)(b2mult); \
      f = tinyecm128(t, B1_, B2_, curves, NC); \
      if (f && show) {show=0;printf("tinyecm128 %uk/%uk found factor %s\n",eb1,eb1*b2mult,u128_str(f));} \
      NC += curves; \
    }
    int NC = 0;
    /*              B1  *=B2  # curves */
    TINYECM_STEP(   2U,   20,  40)
    TINYECM_STEP(  10U,   20,  40)
    TINYECM_STEP(  40U,   20,  40)
    TINYECM_STEP(  80U,   20,  40)
    TINYECM_STEP( 100U,   40, 120)
    TINYECM_STEP( 800U,   40,  80) /* Failsafe.  Extremely unlikely to hit. */
#undef TINYECM_STEP

    if (f <= 1 || f >= t)
      croak("internal: factorintp128 failed to factor %s\n", u128_str(t));
    if ((t % f) != 0)
      croak("internal: factorintp128 found non-factor %s for %s\n", u128_str(f), u128_str(t));

    /* Push both halves (f and t/f) for further factoring */
    stack[top++] = f;
    stack[top++] = t / f;
  }
  if (verbose) fflush(stdout);
}

#endif /* HAVE_FACTOR128 */