Compiler Explorer

Source code

#include <cassert>
#include <cstdint>
#include <vector>

#ifdef _MSC_VER
#include <bit>
#include <__msvc_int128.hpp>
using uint128 = std::_Unsigned128;
#define __builtin_ctzll(x) std::countr_zero(x)
#else
using uint128 = unsigned __int128;
#endif

static_assert(sizeof(uintmax_t) == 8);
static_assert(uintmax_t(-1) >> 63 == 1);

/* 2*3*5*7*11...*101 is 128 bits, and has 26 prime factors */
#define MAX_NFACTS 26

#define FIRST_OMITTED_PRIME 503

struct factors {
  uintmax_t     plarge[2]; /* Can have a single large factor */
  uintmax_t     p[MAX_NFACTS];
  unsigned char e[MAX_NFACTS];
  unsigned char nfactors;
};

constexpr void factor(uintmax_t t0, struct factors *factors);

constexpr uintmax_t gcd_odd(uintmax_t a, uintmax_t b) {
  if ((b & 1) == 0) {
    uintmax_t t = b;
    b = a;
    a = t;
  }
  if (a == 0)
    return b;

/* Take out least significant one bit, to make room for sign */
  b >>= 1;

while (true) {
    a >>= __builtin_ctzll(a) + 1;
    if (a == b) {
      return (a << 1) + 1;
    } else if (a > b) {
      a = (a - b);
    } else {
      uintmax_t olda = a;
      a = (b - a);
      b = olda;
    }
  }
}

static constexpr void
factor_insert_multiplicity(struct factors *factors, uintmax_t prime, int m) {
  int nfactors = factors->nfactors;
  uintmax_t *p = factors->p;
  unsigned char *e = factors->e;

/* Locate position for insert new or increment e.  */
  int i;
  for (i = nfactors - 1; i >= 0; --i) {
    if (p[i] <= prime)
      break;
  }

if (i < 0 || p[i] != prime) {
    for (int j = nfactors - 1; j > i; --j) {
      p[j + 1] = p[j];
      e[j + 1] = e[j];
    }
    p[i + 1] = prime;
    e[i + 1] = m;
    factors->nfactors = nfactors + 1;
  } else {
    e[i] += m;
  }
}

static constexpr unsigned char primes_diff[] = {
#define P(a,b,c,d) a,
  PRIMES_H_CONTENTS
#undef P
  0,0,0,0,0,0,0                           /* 7 sentinels for 8-way loop */
};

#define PRIMES_PTAB_ENTRIES \
  (sizeof(primes_diff) / sizeof(primes_diff[0]) - 8 + 1)

static constexpr unsigned char primes_diff8[] = {
#define P(a,b,c,d) b,
  PRIMES_H_CONTENTS
#undef P
  0,0,0,0,0,0,0                           /* 7 sentinels for 8-way loop */
};

struct primes_dtab {
  uintmax_t binv;
  uintmax_t lim;
};

static constexpr struct primes_dtab primes_dtab[] = {
#define P(a,b,c,d) {c,d},
  PRIMES_H_CONTENTS
#undef P
  {1,0},{1,0},{1,0},{1,0},{1,0},{1,0},{1,0} /* 7 sentinels for 8-way loop */
};

static constexpr void
factor_insert_refind(struct factors *factors, uintmax_t p, int i, int off) {
  for (int j = 0; j < off; ++j)
    p += primes_diff[i + j];
  factor_insert_multiplicity(factors, p, 1);
}

/* Trial division with odd primes uses the following trick.

Let p be an odd prime, and B = 2^{64}.  For simplicity,
   consider the case t < B (this is the second loop below).

From our tables we get

binv = p^{-1} (mod B)
     lim = floor ((B-1) / p).

First assume that t is a multiple of p, t = q * p.  Then 0 <= q <= lim
   (and all quotients in this range occur for some t).

Then t = q * p is true also (mod B), and p is invertible we get

q = t * binv (mod B).

Next, assume that t is *not* divisible by p.  Since multiplication
   by binv (mod B) is a one-to-one mapping,

t * binv (mod B) > lim,

because all the smaller values are already taken.

This can be summed up by saying that the function

q(t) = binv * t (mod B)

is a permutation of the range 0 <= t < B, with the curious property
   that it maps the multiples of p onto the range 0 <= q <= lim, in
   order, and the non-multiples of p onto the range lim < q < B.
 */

static constexpr uintmax_t
factor_using_division(uintmax_t t0, struct factors *factors) {
  if (t0 % 2 == 0) {
    int cnt = __builtin_ctzll(t0);
    t0 >>= cnt;
    factor_insert_multiplicity(factors, 2, cnt);
  }

#define DIVBLOCK(I) do { \
    while (true) { \
      uintmax_t q = t0 * pd[I].binv;                                \
      if ((q > pd[I].lim)) [[likely]]                               \
        break;                                                      \
      t0 = q;                                                       \
      factor_insert_refind(factors, p, i + 1, I);                   \
    }                                                               \
  } while (0)

uintmax_t p = 3;
  for (size_t i = 0; i < PRIMES_PTAB_ENTRIES; i += 8) {
    const struct primes_dtab *pd = &primes_dtab[i];
    DIVBLOCK(0);
    DIVBLOCK(1);
    DIVBLOCK(2);
    DIVBLOCK(3);
    DIVBLOCK(4);
    DIVBLOCK(5);
    DIVBLOCK(6);
    DIVBLOCK(7);
    p += primes_diff8[i];
    if (p * p > t0)
      break;
  }
  return t0;
}

/* Entry i contains (2i+1)^(-1) mod 2^8.  */
static constexpr unsigned char  binvert_table[128] = {
  0x01, 0xAB, 0xCD, 0xB7, 0x39, 0xA3, 0xC5, 0xEF,
  0xF1, 0x1B, 0x3D, 0xA7, 0x29, 0x13, 0x35, 0xDF,
  0xE1, 0x8B, 0xAD, 0x97, 0x19, 0x83, 0xA5, 0xCF,
  0xD1, 0xFB, 0x1D, 0x87, 0x09, 0xF3, 0x15, 0xBF,
  0xC1, 0x6B, 0x8D, 0x77, 0xF9, 0x63, 0x85, 0xAF,
  0xB1, 0xDB, 0xFD, 0x67, 0xE9, 0xD3, 0xF5, 0x9F,
  0xA1, 0x4B, 0x6D, 0x57, 0xD9, 0x43, 0x65, 0x8F,
  0x91, 0xBB, 0xDD, 0x47, 0xC9, 0xB3, 0xD5, 0x7F,
  0x81, 0x2B, 0x4D, 0x37, 0xB9, 0x23, 0x45, 0x6F,
  0x71, 0x9B, 0xBD, 0x27, 0xA9, 0x93, 0xB5, 0x5F,
  0x61, 0x0B, 0x2D, 0x17, 0x99, 0x03, 0x25, 0x4F,
  0x51, 0x7B, 0x9D, 0x07, 0x89, 0x73, 0x95, 0x3F,
  0x41, 0xEB, 0x0D, 0xF7, 0x79, 0xE3, 0x05, 0x2F,
  0x31, 0x5B, 0x7D, 0xE7, 0x69, 0x53, 0x75, 0x1F,
  0x21, 0xCB, 0xED, 0xD7, 0x59, 0xC3, 0xE5, 0x0F,
  0x11, 0x3B, 0x5D, 0xC7, 0x49, 0x33, 0x55, 0xFF
};

/* Compute n^(-1) mod B, using a Newton iteration.  */
#define binv(inv,n) do { \
    uintmax_t  __n = (n); \
    uintmax_t __inv = binvert_table[(__n / 2) & 0x7F]; \
    __inv = 2 * __inv - __inv * __inv * __n;     \
    __inv = 2 * __inv - __inv * __inv * __n;     \
    __inv = 2 * __inv - __inv * __inv * __n;     \
    (inv) = __inv;                               \
  } while (0)

/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B.
   Both a and b must be in redc form, the result will be in redc form too.  */
[[gnu::const]] constexpr uintmax_t
mulredc(uintmax_t a, uintmax_t b, uintmax_t m, uintmax_t mi)
{
  uint128 product = uint128(a) * uint128(b);
  uintmax_t rh = uintmax_t(product >> 64); 
  uintmax_t rl = uintmax_t(product); 
  uintmax_t q = rl * mi;
  product = uint128(q) * uint128(m);
  uintmax_t th = uintmax_t(product >> 64);
  if (rh < th) {
    return (rh - th) + m;
  } else {
    return (rh - th);
  }
}

[[gnu::const]] constexpr uintmax_t
powm(uintmax_t b, uintmax_t e, uintmax_t n, uintmax_t ni, uintmax_t one) {
  uintmax_t y = one;
  if (e & 1)
    y = b;
  while (e != 0) {
    b = mulredc(b, b, n, ni);
    e >>= 1;
    if (e & 1)
      y = mulredc(y, b, n, ni);
  }
  return y;
}

[[gnu::const]] constexpr bool
millerrabin(uintmax_t n, uintmax_t ni, uintmax_t b, uintmax_t q, int k, uintmax_t one) {
  uintmax_t y = powm(b, q, n, ni, one);
  uintmax_t nm1 = n - one;      /* -1, but in redc representation.  */
  if (y == one || y == nm1)
    return true;
  for (int i = 1; i < k; ++i) {
    y = mulredc(y, y, n, ni);
    if (y == nm1)
      return true;
    if (y == one)
      return false;
  }
  return false;
}

/* Lucas' prime test.  The number of iterations vary greatly, up to a few dozen
   have been observed.  The average seem to be about 2.  */
[[gnu::const]] constexpr bool prime_p(uintmax_t n) {
  uintmax_t ni;
  struct factors factors;

if (n <= 1)
    return false;

/* We have already cast out small primes.  */
  if (n < FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME)
    return true;

/* Precomputation for Miller-Rabin.  */
  uintmax_t q = n - 1;
  int k = __builtin_ctzll(q);
  q >>= k;

uintmax_t a = 2;
  binv(ni, n);                 /* ni <- 1/n mod B */
  uintmax_t one = uintmax_t((uint128(1) << 64) % n);
  uintmax_t a_prim = uintmax_t((uint128(one) + one) % n);

/* Perform a Miller-Rabin test, finds most composites quickly.  */
  if (!millerrabin(n, ni, a_prim, q, k, one))
    return false;

/* Factor n-1 for Lucas.  */
  factor(n-1, &factors);

/* Loop until Lucas proves our number prime, or Miller-Rabin proves our
     number composite.  */
  for (size_t r = 0; r < PRIMES_PTAB_ENTRIES; ++r) {
    bool is_prime = true;

for (int i = 0; i < factors.nfactors && is_prime; ++i) {
      is_prime = powm(a_prim, (n - 1) / factors.p[i], n, ni, one) != one;
    }
    if (is_prime)
      return true;

a += primes_diff[r];      /* Establish new base.  */

a_prim = uintmax_t((uint128(one) * a) % n);
    if (!millerrabin(n, ni, a_prim, q, k, one))
      return false;

// Miller-Rabin for a=2,3,5,7,11,13,17,19,23,29,31,37 is enough for n < 2^64
    if (r == 10) return true;
  }

assert(!"Lucas prime test failure.  This should not happen");
}

constexpr void
factor_using_pollard_rho(uintmax_t n, unsigned long int a, struct factors *factors) {
  uintmax_t t, ni, g;

unsigned long int k = 1;
  unsigned long int l = 1;

uintmax_t P = uintmax_t((uint128(1) << 64) % n);
  uintmax_t x = uintmax_t((uint128(P) + P) % n);
  uintmax_t y = x;
  uintmax_t z = x;

while (n != 1) {
      assert(a < n);
      binv(ni, n);             /* FIXME: when could we use old 'ni' value?  */
      [&]() {
        while (true) {
          do {
            x = mulredc(x, x, n, ni);
            x = uintmax_t((uint128(x) + a) % n);
            t = uintmax_t((uint128(z) - x) % n);
            P = mulredc(P, t, n, ni);
            if (k % 32 == 1) {
              if (gcd_odd(P, n) != 1)
                return; // goto factor_found;
              y = x;
            }
          } while (--k != 0);

z = x;
          k = l;
          l = 2 * l;
          for (unsigned long int i = 0; i < k; i++) {
            x = mulredc(x, x, n, ni);
            x = uintmax_t((uint128(x) + a) % n);
          }
          y = x;
        }
      }();

// factor_found:
      do {
        y = mulredc(y, y, n, ni);
        y = uintmax_t((uint128(y) + a) % n);
        t = uintmax_t((uint128(z) - y) % n);
        g = gcd_odd(t, n);
      } while (g == 1);

if (n == g) {
        /* Found n itself as factor.  Restart with different params.  */
        factor_using_pollard_rho(n, a + 1, factors);
        return;
      }

n /= g;

if (prime_p(g)) {
        factor_insert_multiplicity(factors, g, 1);
      } else {
        factor_using_pollard_rho(g, a + 1, factors);
      }

if (prime_p(n)) {
        factor_insert_multiplicity(factors, n, 1);
        break;
      }

x %= n;
      z %= n;
      y %= n;
    }
}

/* Compute the prime factors of the 128-bit number (T1,T0), and put the
   results in FACTORS.  */
constexpr void factor(uintmax_t t0, struct factors *factors) {
  factors->nfactors = 0;
  factors->plarge[1] = 0;
  if (t0 < 2)
    return;
  t0 = factor_using_division(t0, factors);
  if (t0 < 2)
    return;
  if (prime_p(t0))
    factor_insert_multiplicity(factors, t0, 1);
  else
    factor_using_pollard_rho(t0, 1, factors);
}

constexpr
std::vector<uintmax_t> factors_of(uintmax_t t0) {
    struct factors factors;
    factor(t0, &factors);
    std::vector<uintmax_t> result;
    for (int j = 0; j < factors.nfactors; ++j) {
        for (int k = 0; k < factors.e[j]; ++k) {
            result.push_back(factors.p[j]);
        }
    }
    return result;
}

using V = std::vector<uintmax_t>;

// https://oeis.org/A014233
static_assert((factors_of(2047) == V{23, 89}));
static_assert((factors_of(1373653) == V{829, 1657}));
static_assert((factors_of(25326001) == V{2251, 11251}));
static_assert((factors_of(3215031751) == V{151, 751, 28351}));
static_assert((factors_of(2152302898747) == V{6763, 10627, 29947}));
static_assert((factors_of(3474749660383) == V{1303, 16927, 157543}));
static_assert((factors_of(341550071728321) == V{10670053, 32010157}));
static_assert((factors_of(3825123056546413051) == V{149491, 747451, 34233211}));

// Miscellaneous test cases
static_assert((factors_of(100) == V{2, 2, 5, 5}));
static_assert((factors_of(2147462143) == V{2147462143}));
static_assert((factors_of(2147462142) == V{2, 3, 7, 7, 757, 9649}));
static_assert((factors_of(15032235001) == V{7, 2147462143}));
static_assert((factors_of(214746214200) == V{2, 2, 2, 3, 5, 5, 7, 7, 757, 9649}));
static_assert((factors_of(18361375334787046697uLL) == V{18361375334787046697uLL}));
//static_assert((factors_of(4611593655618152449) == V{2147462143, 2147462143}));
//static_assert((factors_of(14862346268936648689uLL) == V{3855171367, 3855171367}));