Últimos N lucky numbers bajo M. Si encuentras un algoritmo mejor, please, cuéntamelo!

fable.last20a.cpp

// fable.last20a.cpp
// Últimos 20 lucky numbers <= N. Misma interfaz y salida que la referencia:
//   ./lucky_turbo N      -> imprime los últimos 20 lucky <= N, uno por línea, ascendente.
//
// Compilar:
//   g++ -O3 -march=native -mbmi2 -std=c++17 -fopenmp fable.last20a.cpp -o fable.last20a.cpp.exe
//
//
//
// $ head /proc/cpuinfo
// vendor_id	: AuthenticAMD
// model name	: AMD Ryzen 5 5500
// cpu MHz		: 3289.786
// cache size	: 512 KB
// $ time -f "%e, %M" ./fable.last20a.cpp.exe 100000000000
// 99999999561
// 99999999571
// 99999999585
// 99999999645
// 99999999687
// 99999999691
// 99999999721
// 99999999739
// 99999999765
// 99999999771
// 99999999781
// 99999999783
// 99999999793
// 99999999837
// 99999999897
// 99999999921
// 99999999951
// 99999999955
// 99999999967
// 99999999973
// 62.96, 3981404
//
// ESTRATEGIA (cambio algorítmico, no micro-optimización):
//
//   La referencia criba TODO [1..N]. Pero solo necesitamos los últimos 20.
//   Un n impar se puede testear individualmente sin cribar [1..N]:
//     p = (n+1)/2  (posición 1-indexada entre los impares)
//     para cada lucky L = 3, 7, 9, 13, ... en orden ascendente:
//        si L > p        -> n ES lucky (ningún paso futuro puede tocarlo: p < L y p%L != 0)
//        si p % L == 0   -> n NO es lucky (eliminado en la pasada de L)
//        si no           -> p -= p/L   (nueva posición tras quitar cada L-ésimo)
//
//   El test de un n cerca de N solo necesita los lucky hasta ~pi_lucky(N) ~ N/ln N,
//   es decir, una criba ~ln N veces MÁS PEQUEÑA (para N=1e9: 6e7 en vez de 1e9;
//   el bitmap pasa de 62.5 MB en DRAM a 3.8 MB residente en L3).
//
//   Coste: criba(1.25*N/ln N) + ~20/densidad ~ 200 walks (paralelos, división exacta
//   por recíprocos de Lemire: q = hi64(p*M), exacto si p*L < 2^64, garantizado porque
//   en la fase de recíprocos p < 2^32 y L <= p).
//
//   Corrección del recorte de la criba: los supervivientes de cribar [1..M] son
//   EXACTAMENTE los lucky <= M, porque una pasada con paso k solo elimina valores
//   >= 2k-1 (el k-ésimo superviviente vale al menos 2k-1), así que las pasadas con
//   k grande del proceso completo jamás tocan [1..M].
//
//   Red de seguridad: si algún test agota la lista (no decidible), se duplica M y
//   se repite. Con M = 1.25*N/ln N + 4096 no ocurre en la práctica (pi_lucky(x) es
//   ligeramente menor que pi_primos(x) ~ 1.05*x/ln x), pero garantiza corrección
//   incondicional.
//
//   La criba pequeña reutiliza el mismo motor de dos regímenes de la referencia
//   (streaming byte-table/pdep para k<=4096, Fenwick select+delete para k grande),
//   que ahora corre entero en caché L3.

#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#include <immintrin.h>

#ifdef _OPENMP
  #include <omp.h>
#endif

static inline void die(const char* msg) { std::perror(msg); std::exit(1); }

static inline void* xaligned_alloc(size_t align, size_t size) {
  size = (size + align - 1) & ~(align - 1);
  void* p = nullptr;
  if (posix_memalign(&p, align, size) != 0) die("posix_memalign");
  return p;
}

static inline uint32_t pop64(uint64_t x) { return (uint32_t)__builtin_popcountll(x); }
static inline uint32_t ctz64(uint64_t x) { return (uint32_t)__builtin_ctzll(x); }
static inline uint32_t clz64(uint64_t x) { return (uint32_t)__builtin_clzll(x); }

// --- jerarquía: words(64) -> groups(16 words) -> superbloques(128 groups) ---
static constexpr uint32_t GROUP_WORDS = 16;
static constexpr uint32_t SB_GROUPS   = 128;

static constexpr uint32_t TAB_MAX   = 256;
static constexpr uint32_t SMALL_MAX = 4096;

#ifndef DIRECT_MAX_OVERRIDE
static constexpr uint64_t DIRECT_MAX = 1ull << 20; // por debajo, criba directa (ya es trivial)
#else
static constexpr uint64_t DIRECT_MAX = DIRECT_MAX_OVERRIDE;
#endif

struct ByteEntry {
  uint8_t out;
  uint8_t outpop;
  uint16_t next; // 0..k-1
};

static void build_byte_table(uint32_t k, std::vector<ByteEntry>& tab) {
  tab.resize((size_t)k * 256);
  for (uint32_t s = 0; s < k; ++s) {
    for (uint32_t m = 0; m < 256; ++m) {
      uint32_t st = s;
      uint8_t out = (uint8_t)m;
      for (uint32_t b = 0; b < 8; ++b) {
        if (m & (1u << b)) {
          if (st == k - 1) { out = (uint8_t)(out & ~(1u << b)); st = 0; }
          else { ++st; }
        }
      }
      ByteEntry e;
      e.out = out;
      e.outpop = (uint8_t)__builtin_popcount((unsigned)out);
      e.next = (uint16_t)st;
      tab[s * 256 + m] = e;
    }
  }
}

static inline uint64_t word_delete_every_kth(uint64_t w, uint32_t pc_in, uint32_t k, uint32_t& st_inout) {
  const uint32_t s0 = st_inout;
  const uint32_t q  = s0 + pc_in;
  const uint32_t d  = q / k;
  st_inout = q - d * k;

  if (d == 0 || w == 0) return w;

  uint64_t del = 0;
  uint32_t rank = (s0 == 0) ? k : (k - s0);
  for (uint32_t i = 0; i < d; ++i) {
    del |= _pdep_u64(1ULL << (rank - 1), w);
    rank += k;
  }
  return w & ~del;
}

static inline uint64_t find_next_set_bit(const uint64_t* bits, uint64_t total_bits, uint64_t start_bit) {
  if (start_bit >= total_bits) return ~0ULL;
  uint64_t wi = start_bit >> 6;
  uint32_t bi = (uint32_t)(start_bit & 63);
  uint64_t w = bits[wi] & (~0ULL << bi);
  while (true) {
    if (w) return (wi << 6) + (uint64_t)ctz64(w);
    ++wi;
    uint64_t bitpos = wi << 6;
    if (bitpos >= total_bits) break;
    w = bits[wi];
  }
  return ~0ULL;
}

// ---- Fenwick sobre superbloques ----
struct Fenwick {
  uint32_t n;
  std::vector<uint64_t> bit;

  explicit Fenwick(uint32_t n_ = 0) : n(n_), bit((size_t)n_ + 1, 0) {}

  void reset(uint32_t n_) { n = n_; bit.assign((size_t)n_ + 1, 0); }

  void build_from(const uint32_t* sbcnt) {
    std::fill(bit.begin(), bit.end(), 0);
    for (uint32_t i = 1; i <= n; ++i) bit[i] = sbcnt[i - 1];
    for (uint32_t i = 1; i <= n; ++i) {
      uint32_t j = i + (i & -i);
      if (j <= n) bit[j] += bit[i];
    }
  }

  inline void add(uint32_t idx1, int32_t delta) {
    if (delta >= 0) {
      uint64_t d = (uint64_t)delta;
      for (uint32_t i = idx1; i <= n; i += (i & -i)) bit[i] += d;
    } else {
      uint64_t d = (uint64_t)(-(int64_t)delta);
      for (uint32_t i = idx1; i <= n; i += (i & -i)) bit[i] -= d;
    }
  }

  inline uint32_t find_by_order(uint64_t rank, uint64_t& rem) const {
    uint32_t idx = 0;
    uint32_t bitmask = 1u << (31 - __builtin_clz(n));
    uint64_t r = rank;
    while (bitmask) {
      uint32_t nxt = idx + bitmask;
      if (nxt <= n && bit[nxt] < r) {
        idx = nxt;
        r -= bit[nxt];
      }
      bitmask >>= 1;
    }
    rem = r;
    return idx + 1;
  }
};

// ---- criba completa hasta N (el motor de la referencia, con k en 64 bits) ----
struct Sieve {
  uint64_t total_bits = 0, n_words = 0, n_groups = 0;
  uint32_t n_sbs = 0;
  uint64_t* bits  = nullptr;
  uint8_t*  wcnt  = nullptr;
  uint16_t* gcnt  = nullptr;
  uint32_t* sbcnt = nullptr;

  void release() {
    std::free(bits);  bits  = nullptr;
    std::free(wcnt);  wcnt  = nullptr;
    std::free(gcnt);  gcnt  = nullptr;
    std::free(sbcnt); sbcnt = nullptr;
  }
};

static void run_lucky_sieve(uint64_t N, Sieve& S) {
  const uint64_t total_bits = (N + 1) >> 1;
  const uint64_t n_words  = (total_bits + 63) >> 6;
  const uint64_t n_groups = (n_words + GROUP_WORDS - 1) / GROUP_WORDS;
  const uint32_t n_sbs    = (uint32_t)((n_groups + SB_GROUPS - 1) / SB_GROUPS);

  S.total_bits = total_bits; S.n_words = n_words; S.n_groups = n_groups; S.n_sbs = n_sbs;
  S.bits  = (uint64_t*)xaligned_alloc(64, n_words * sizeof(uint64_t));
  S.wcnt  = (uint8_t*) xaligned_alloc(64, n_words * sizeof(uint8_t));
  S.gcnt  = (uint16_t*)xaligned_alloc(64, n_groups * sizeof(uint16_t));
  S.sbcnt = (uint32_t*)xaligned_alloc(64, (size_t)n_sbs * sizeof(uint32_t));

  uint64_t* bits  = S.bits;
  uint8_t*  wcnt  = S.wcnt;
  uint16_t* gcnt  = S.gcnt;
  uint32_t* sbcnt = S.sbcnt;

  std::memset(bits, 0xFF, n_words * sizeof(uint64_t));
  uint64_t remb = total_bits & 63;
  if (remb) bits[n_words - 1] = (1ULL << remb) - 1ULL;

#ifdef _OPENMP
  #pragma omp parallel for schedule(static)
#endif
  for (int64_t i = 0; i < (int64_t)n_words; ++i) wcnt[i] = (uint8_t)pop64(bits[i]);

#ifdef _OPENMP
  #pragma omp parallel for schedule(static)
#endif
  for (int64_t g = 0; g < (int64_t)n_groups; ++g) {
    uint32_t sum = 0;
    uint64_t w0 = (uint64_t)g * GROUP_WORDS;
    uint64_t w1 = std::min<uint64_t>(n_words, w0 + GROUP_WORDS);
    for (uint64_t w = w0; w < w1; ++w) sum += wcnt[w];
    gcnt[g] = (uint16_t)sum;
  }

#ifdef _OPENMP
  #pragma omp parallel for schedule(static)
#endif
  for (int64_t sb = 0; sb < (int64_t)n_sbs; ++sb) {
    uint32_t sum = 0;
    uint64_t g0 = (uint64_t)sb * SB_GROUPS;
    uint64_t g1 = std::min<uint64_t>(n_groups, g0 + SB_GROUPS);
    for (uint64_t g = g0; g < g1; ++g) sum += gcnt[g];
    sbcnt[sb] = sum;
  }

  auto sum_alive = [&]() -> uint64_t {
    uint64_t s = 0;
    for (uint32_t sb = 0; sb < n_sbs; ++sb) s += sbcnt[sb];
    return s;
  };

  uint64_t alive = sum_alive();

  Fenwick fw(n_sbs);
  fw.build_from(sbcnt);

  uint64_t last_k_idx = 0;

  std::vector<uint64_t> sb_prefix(n_sbs + 1);
  std::vector<uint32_t> sb_start_mod(n_sbs);
  std::vector<ByteEntry> byte_tab;

  while (true) {
    uint64_t idx = find_next_set_bit(bits, total_bits, last_k_idx + 1);
    if (idx == ~0ULL) break;
    const uint64_t k = idx * 2 + 1; // 64 bits: sin overflow para N grandes
    if (k > alive) break;
    last_k_idx = idx;

    if (k <= SMALL_MAX) {
      const uint32_t k32 = (uint32_t)k;
      sb_prefix[0] = 0;
      for (uint32_t sb = 0; sb < n_sbs; ++sb) sb_prefix[sb + 1] = sb_prefix[sb] + sbcnt[sb];
      for (uint32_t sb = 0; sb < n_sbs; ++sb) sb_start_mod[sb] = (uint32_t)(sb_prefix[sb] % k32);

      const bool use_tab = (k32 <= TAB_MAX);
      if (use_tab) build_byte_table(k32, byte_tab);

#ifdef _OPENMP
      #pragma omp parallel for schedule(static)
#endif
      for (int64_t sb = 0; sb < (int64_t)n_sbs; ++sb) {
        uint32_t st = sb_start_mod[(size_t)sb];

        const uint64_t g0 = (uint64_t)sb * SB_GROUPS;
        const uint64_t g1 = std::min<uint64_t>(n_groups, g0 + SB_GROUPS);

        uint32_t sb_new = 0;

        for (uint64_t g = g0; g < g1; ++g) {
          const uint64_t w0 = g * GROUP_WORDS;
          const uint64_t w1 = std::min<uint64_t>(n_words, w0 + GROUP_WORDS);

          uint32_t g_new = 0;

          for (uint64_t w = w0; w < w1; ++w) {
            uint64_t x = bits[w];
            uint32_t pc_in = (uint32_t)wcnt[w];
            if (!x || pc_in == 0) { wcnt[w] = 0; continue; }

            uint64_t y = 0;
            if (use_tab) {
              uint32_t state = st;
              uint32_t pc_out = 0;
              uint64_t outw = 0;
              for (uint32_t bi = 0; bi < 8; ++bi) {
                uint8_t inb = (uint8_t)(x >> (bi * 8));
                const ByteEntry e = byte_tab[state * 256 + inb];
                outw |= (uint64_t)e.out << (bi * 8);
                pc_out += e.outpop;
                state = e.next;
              }
              y = outw;
              st = state;
              wcnt[w] = (uint8_t)pc_out;
              g_new += pc_out;
            } else {
              uint32_t state = st;
              y = word_delete_every_kth(x, pc_in, k32, state);
              st = state;
              uint32_t pc_out = pop64(y);
              wcnt[w] = (uint8_t)pc_out;
              g_new += pc_out;
            }

            bits[w] = y;
          }

          gcnt[g] = (uint16_t)g_new;
          sb_new += g_new;
        }

        sbcnt[(size_t)sb] = sb_new;
      }

      alive = sum_alive();
      fw.build_from(sbcnt);

    } else {
      const uint64_t m = alive / k;
      if (m == 0) continue;

      for (uint64_t t = m; t >= 1; --t) {
        uint64_t rank = t * k;
        uint64_t rem_in_sb = 0;
        uint32_t sb1 = fw.find_by_order(rank, rem_in_sb);
        uint32_t sb = sb1 - 1;

        uint64_t g0 = (uint64_t)sb * SB_GROUPS;
        uint64_t g1 = std::min<uint64_t>(n_groups, g0 + SB_GROUPS);

        uint64_t r = rem_in_sb;
        uint64_t g_found = g0;
        for (uint64_t g = g0; g < g1; ++g) {
          uint32_t c = gcnt[g];
          if (r > c) r -= c;
          else { g_found = g; break; }
        }

        uint64_t w0 = g_found * GROUP_WORDS;
        uint64_t w1 = std::min<uint64_t>(n_words, w0 + GROUP_WORDS);

        uint64_t rr = r;
        uint64_t w_found = w0;
        for (uint64_t w = w0; w < w1; ++w) {
          uint32_t c = wcnt[w];
          if (rr > c) rr -= c;
          else { w_found = w; break; }
        }

        uint64_t word = bits[w_found];
        uint64_t mask = _pdep_u64(1ULL << (uint32_t)(rr - 1), word);

        bits[w_found] = word & ~mask;
        wcnt[w_found]--;
        gcnt[g_found]--;
        sbcnt[sb]--;
        fw.add(sb1, -1);

        if (t == 1) break;
      }

      alive -= m;
    }
  }
}

static void print_last20_from_bitmap(const Sieve& S) {
  std::vector<uint64_t> last;
  last.reserve(20);
  for (int64_t wi = (int64_t)S.n_words - 1; wi >= 0 && (int)last.size() < 20; --wi) {
    uint64_t w = S.bits[wi];
    while (w && (int)last.size() < 20) {
      uint32_t b = 63u - clz64(w);
      uint64_t idx = ((uint64_t)wi << 6) + b;
      if (idx < S.total_bits) last.push_back(idx * 2 + 1);
      w &= ~(1ULL << b);
    }
  }
  std::reverse(last.begin(), last.end());
  for (uint64_t x : last) std::printf("%llu\n", (unsigned long long)x);
}

// ---- test individual de lucky por "walk" de posiciones ----
// Devuelve 1 = lucky, 0 = no lucky, -1 = lista insuficiente (no decidible).
static int test_lucky(uint64_t n, const uint64_t* __restrict Ls,
                      const uint64_t* __restrict Ms, size_t nl) {
  uint64_t p = (n + 1) >> 1;
  size_t i = 1; // saltar el lucky 1 (nunca se usa como paso)

  // Fase 64 bits (p >= 2^32): división hardware; los recíprocos no son
  // exactos en general aquí. Solo se activa para N >~ 8.5e9.
  while (p > 0xFFFFFFFFull) {
    if (i >= nl) return -1;
    const uint64_t L = Ls[i];
    if (L > p) return 1;
    const uint64_t q = p / L;
    if (p - q * L == 0) return 0;
    p -= q;
    ++i;
  }

  // Fase recíprocos: p < 2^32 y L <= p  =>  p*L < 2^64  =>  q exacto
  // (Lemire: M = floor(2^64/L)+1, q = hi64(p*M)).
  for (; i < nl; ++i) {
    const uint64_t L = Ls[i];
    if (L > p) return 1;
    const uint64_t q = (uint64_t)(((unsigned __int128)p * Ms[i]) >> 64);
    if (p == q * L) return 0; // resto cero -> eliminado
    p -= q;
  }
  return -1; // lista agotada con L <= p: hay que ampliar M
}

int main(int argc, char** argv) {
  if (argc != 2) { std::fprintf(stderr, "Usage: %s N\n", argv[0]); return 1; }
  const uint64_t N = std::strtoull(argv[1], nullptr, 10);
  if (N == 0) return 0;

  if (N <= DIRECT_MAX) {
    Sieve S;
    run_lucky_sieve(N, S);
    print_last20_from_bitmap(S);
    return 0;
  }

  const double lnN = std::log((double)N);
  uint64_t Mb = (uint64_t)(1.25 * (double)N / lnN) + 4096;

  for (;;) {
    if (Mb >= N) { // red de seguridad extrema: equivale a la criba completa
      Sieve S;
      run_lucky_sieve(N, S);
      print_last20_from_bitmap(S);
      return 0;
    }

    Sieve S;
    run_lucky_sieve(Mb, S);

    // extraer la lista de lucky <= Mb y sus recíprocos
    std::vector<uint64_t> wpref(S.n_words + 1);
    wpref[0] = 0;
    for (uint64_t w = 0; w < S.n_words; ++w) wpref[w + 1] = wpref[w] + S.wcnt[w];
    const uint64_t nl = wpref[S.n_words];

    std::vector<uint64_t> Ls(nl), Ms(nl);
#ifdef _OPENMP
    #pragma omp parallel for schedule(static)
#endif
    for (int64_t w = 0; w < (int64_t)S.n_words; ++w) {
      uint64_t x = S.bits[w];
      uint64_t j = wpref[w];
      while (x) {
        uint32_t b = ctz64(x);
        uint64_t L = ((((uint64_t)w << 6) + b) << 1) + 1;
        Ls[j] = L;
        Ms[j] = (L >= 2) ? (UINT64_MAX / L + 1) : 0; // floor(2^64/L)+1 (L impar)
        ++j;
        x &= x - 1;
      }
    }
    S.release();

    // testear candidatos n = N, N-2, ... en paralelo hasta reunir 20 lucky
    std::vector<uint64_t> found;
    found.reserve(64);
    bool undecided = false;
    int64_t n = (int64_t)((N & 1) ? N : N - 1);
    constexpr int CH = 512;
    int8_t st[CH];

    while ((int)found.size() < 20 && n >= 1 && !undecided) {
      const int C = (int)std::min<int64_t>(CH, (n + 1) / 2);
#ifdef _OPENMP
      #pragma omp parallel for schedule(dynamic, 4)
#endif
      for (int i = 0; i < C; ++i)
        st[i] = (int8_t)test_lucky((uint64_t)(n - 2 * (int64_t)i), Ls.data(), Ms.data(), nl);

      for (int i = 0; i < C; ++i) {
        if (st[i] < 0) { undecided = true; break; }
        if (st[i] == 1) {
          found.push_back((uint64_t)(n - 2 * (int64_t)i));
          if ((int)found.size() == 20) break;
        }
      }
      n -= 2 * (int64_t)C;
    }

    if (!undecided) {
      if ((int)found.size() > 20) found.resize(20);
      std::reverse(found.begin(), found.end());
      for (uint64_t x : found) std::printf("%llu\n", (unsigned long long)x);
      return 0;
    }

    Mb *= 2; // ampliar lista de lucky y reintentar (en la práctica no ocurre)
  }
}