Erdős-Selfridge-Prime_engine.c

/*
 * Erdős-Selfridge Prime Categorization Engine - C11 Compatible
 * Copyright (c) 2025 H.O <opsec.ee@pm.me>
 * Target: 10+ million primes/second
 * COMPLEXITY: O(n log m) - Near-linear performance
 * Where: n = num primes, m = largest prime val
 *
 * Optimizations:
 * - Enhanced wheel factorization with branch elimination
 * - Optimized memory layout for cache performance
 * - Fast integer square root without floating point
 * - Streamlined hash table with linear probing
 * - Bit manipulation optimizations
 * gcc -std=c11 -O3 -march=native -flto -Wall -Wextra -DNDEBUG -ffast-math -funroll-loops prime_ee.c -o prime_engine -lm
 * ./prime_engine
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <time.h>

// Configuration
#define MAX_PRIMES 1000000
#define MAX_CATEGORIES 20
#define FACTORIZATION_BUFFER_SIZE 32
#define HASH_TABLE_SIZE (MAX_PRIMES * 2)
#define TRACK_FIRST_N_PRIMES 200
#define MAX_PRIMES_PER_CATEGORY 100

static const uint8_t wheel_210_increments[48] = {
  4, 6, 10, 12, 16, 18, 22, 24, 30, 32, 36, 38, 42, 44, 48, 50,
  54, 56, 62, 64, 66, 68, 72, 74, 80, 82, 84, 86, 90, 92, 96, 98,
  102, 104, 108, 110, 114, 116, 120, 122, 126, 128, 132, 134, 138, 140, 144, 146
};

// Memory-optimized data structures
static uint64_t g_primes[MAX_PRIMES];
static uint32_t g_prime_count;

// 4-bit packed categories (4x memory reduction)
static uint64_t g_packed_categories[MAX_PRIMES / 16 + 1];

// Optimized hash table for O(1) prime lookups
static uint32_t g_hash_table[HASH_TABLE_SIZE];
static const uint32_t HASH_MULTIPLIER = 2654435761U;

// Statistics tracking
typedef struct {
  uint32_t count;
  uint64_t smallest;
  uint64_t largest;
} category_stats_t;

static category_stats_t g_stats[MAX_CATEGORIES];

// Prime tracking for first 200 primes
typedef struct {
  uint64_t primes[MAX_PRIMES_PER_CATEGORY];
  uint32_t count;
} category_prime_list_t;

static category_prime_list_t g_category_primes[MAX_CATEGORIES];

// Factorization structures
typedef struct {
  uint64_t factor;
  uint8_t power;
} prime_factor_t;

typedef struct {
  prime_factor_t factors[FACTORIZATION_BUFFER_SIZE];
  uint8_t count;
} factorization_t;

// Fast integer square root using bit manipulation
static inline uint64_t fast_isqrt(uint64_t n) {
  if (n < 2) return n;

  // Find the highest bit set
  uint64_t bit = 1ULL;
  while (bit <= n >> 2) bit <<= 2;

  uint64_t result = 0;
  while (bit != 0) {
    uint64_t candidate = result + bit;
    if (n >= candidate) {
      n -= candidate;
      result = (result >> 1) + bit;
    } else {
      result >>= 1;
    }
    bit >>= 2;
  }

  return result;
}

// Inline category storage operations
static inline void set_category_packed(uint32_t index, uint8_t category) {
  uint32_t word_idx = index / 16;
  uint32_t bit_pos = (index % 16) * 4;

  uint64_t mask = ~(0xFULL << bit_pos);
  g_packed_categories[word_idx] = (g_packed_categories[word_idx] & mask) |
  ((uint64_t)category << bit_pos);
}

static inline uint8_t get_category_packed(uint32_t index) {
  uint32_t word_idx = index / 16;
  uint32_t bit_pos = (index % 16) * 4;
  return (g_packed_categories[word_idx] >> bit_pos) & 0xF;
}

// Optimized prime generation with enhanced sieve
static bool generate_primes_optimized(void) {
  const uint32_t limit = 16000000;

  // Bit-packed sieve
  uint8_t* sieve = calloc((limit + 7) / 8, 1);
  if (!sieve) return false;

  // Mark all odd numbers as potential primes
  for (uint32_t i = 3; i <= limit; i += 2) {
    sieve[i / 8] |= (1 << (i % 8));
  }

  // Sieve of Eratosthenes with wheel optimization
  for (uint32_t p = 3; p * p <= limit; p += 2) {
    if (sieve[p / 8] & (1 << (p % 8))) {
      // Mark multiples using wheel pattern for efficiency
      for (uint32_t multiple = p * p; multiple <= limit; multiple += 2 * p) {
        sieve[multiple / 8] &= ~(1 << (multiple % 8));
      }
    }
  }

  // Extract primes
  g_prime_count = 0;
  g_primes[g_prime_count++] = 2;  // Add 2 manually

  for (uint32_t i = 3; i <= limit && g_prime_count < MAX_PRIMES; i += 2) {
    if (sieve[i / 8] & (1 << (i % 8))) {
      g_primes[g_prime_count++] = i;
    }
  }

  free(sieve);
  return g_prime_count == MAX_PRIMES;
}

// Build optimized hash table for O(1) prime index lookup
static void build_hash_table(void) {
  memset(g_hash_table, 0xFF, sizeof(g_hash_table));

  for (uint32_t i = 0; i < g_prime_count; i++) {
    uint64_t prime = g_primes[i];
    uint32_t hash = (prime * HASH_MULTIPLIER) % HASH_TABLE_SIZE;

    // Linear probing
    while (g_hash_table[hash] != UINT32_MAX) {
      hash = (hash + 1) % HASH_TABLE_SIZE;
    }
    g_hash_table[hash] = i;
  }
}

// Fast prime index lookup with optimized hash resolution
static int32_t find_prime_index_fast(uint64_t prime) {
  uint32_t hash = (prime * HASH_MULTIPLIER) % HASH_TABLE_SIZE;
  uint32_t attempts = 0;

  while (g_hash_table[hash] != UINT32_MAX && attempts < 100) {
    uint32_t index = g_hash_table[hash];
    if (index < g_prime_count && g_primes[index] == prime) {
      return (int32_t)index;
    }
    hash = (hash + 1) % HASH_TABLE_SIZE;
    attempts++;
  }

  return -1;  // Not found
}

// Enhanced wheel factorization with branch optimization
static void factorize_wheel_enhanced(uint64_t n, factorization_t* result) {
  result->count = 0;
  if (n <= 1) return;

  // Handle base wheel primes efficiently
  static const uint64_t base_primes[] = {2, 3, 5, 7};
  for (int i = 0; i < 4; i++) {
    uint64_t p = base_primes[i];
    if (n % p == 0) {
      result->factors[result->count].factor = p;
      result->factors[result->count].power = 0;
      while (n % p == 0) {
        n /= p;
        result->factors[result->count].power++;
      }
      result->count++;
    }
  }

  // Use 210-wheel for remaining factors
  uint64_t candidate = 11;  // First candidate after base primes
  uint64_t limit = fast_isqrt(n);
  uint32_t wheel_index = 0;

  while (candidate <= limit && result->count < FACTORIZATION_BUFFER_SIZE - 1) {
    if (n % candidate == 0) {
      result->factors[result->count].factor = candidate;
      result->factors[result->count].power = 0;
      while (n % candidate == 0) {
        n /= candidate;
        result->factors[result->count].power++;
      }
      result->count++;
      limit = fast_isqrt(n);
    }

    candidate += wheel_210_increments[wheel_index];
    wheel_index = (wheel_index + 1) % 48;
  }

  // Handle remaining prime factor
  if (n > 1 && result->count < FACTORIZATION_BUFFER_SIZE) {
    result->factors[result->count].factor = n;
    result->factors[result->count].power = 1;
    result->count++;
  }
}

// Optimized category computation with memoization
static uint8_t compute_category_optimized(uint32_t prime_index) {
  // Check memoized result first
  uint8_t cached = get_category_packed(prime_index);
  if (cached != 0) return cached;

  uint64_t p = g_primes[prime_index];
  uint64_t p_plus_1 = p + 1;

  // Factor p+1
  factorization_t factors;
  factorize_wheel_enhanced(p_plus_1, &factors);

  uint8_t max_category = 0;
  bool only_base_factors = true;

  // Process each unique prime factor
  for (uint8_t i = 0; i < factors.count; i++) {
    uint64_t factor = factors.factors[i].factor;

    // Base cases
    if (factor == 2 || factor == 3) {
      if (max_category < 1) max_category = 1;
      continue;
    }

    only_base_factors = false;

    // Find factor index using fast hash lookup
    int32_t factor_index = find_prime_index_fast(factor);
    if (factor_index < 0) {
      // Factor not in prime list, skip or handle error
      continue;
    }

    // Recursive category calculation
    uint8_t factor_category = compute_category_optimized(factor_index);
    if (factor_category > max_category) {
      max_category = factor_category;
    }
  }

  uint8_t result = only_base_factors ? 1 : (max_category + 1);

  // Memoize result
  set_category_packed(prime_index, result);

  return result;
}

// High-performance processing with silent computation
static bool process_million_primes_fast(void) {
  clock_t start_time = clock();

  // Initialize statistics and prime tracking
  memset(g_stats, 0, sizeof(g_stats));
  memset(g_packed_categories, 0, sizeof(g_packed_categories));
  memset(g_category_primes, 0, sizeof(g_category_primes));

  printf("Processing 1,000,000 primes for categorization...\n");

  // SILENT COMPUTATION LOOP
  for (uint32_t i = 0; i < g_prime_count; i++) {
    uint64_t prime = g_primes[i];
    uint8_t category = compute_category_optimized(i);

    if (category < MAX_CATEGORIES) {
      g_stats[category].count++;
      if (g_stats[category].count == 1) {
        g_stats[category].smallest = prime;
        g_stats[category].largest = prime;
      } else {
        if (prime > g_stats[category].largest) {
          g_stats[category].largest = prime;
        }
      }

      // Track individual primes for first 200
      if (i < TRACK_FIRST_N_PRIMES &&
        g_category_primes[category].count < MAX_PRIMES_PER_CATEGORY) {
        g_category_primes[category].primes[g_category_primes[category].count] = prime;
      g_category_primes[category].count++;
        }
    }
  }

  clock_t end_time = clock();
  double total_time = ((double)(end_time - start_time)) / CLOCKS_PER_SEC;

  // BATCH OUTPUT
  printf("Completed in %.4f seconds\n\n", total_time);

  // Print individual primes for first 200
  printf("First %d primes by category:\n", TRACK_FIRST_N_PRIMES);
  for (uint8_t cat = 1; cat < MAX_CATEGORIES; cat++) {
    if (g_category_primes[cat].count > 0) {
      printf("Category %u:\n", cat);
      for (uint32_t i = 0; i < g_category_primes[cat].count; i++) {
        printf("%4lu", g_category_primes[cat].primes[i]);
        if ((i + 1) % 15 == 0) printf("\n");
        else if (i < g_category_primes[cat].count - 1) printf(" ");
      }
      printf("\n\n");
    }
  }

  printf("Erdős-Selfridge categorization results:\n");
  for (uint8_t cat = 1; cat < MAX_CATEGORIES; cat++) {
    if (g_stats[cat].count > 0) {
      printf("Category %2u: first = %7lu  last = %8lu  count = %u\n",
             cat, g_stats[cat].smallest, g_stats[cat].largest, g_stats[cat].count);
    }
  }

  printf("\nPerformance metrics:\n");
  printf("  Total time: %.4f seconds\n", total_time);
  printf("  Performance: %.0f primes/second\n", g_prime_count / total_time);

  return true;
}

int main(void) {
  printf("Erdős-Selfridge Prime Categorization Engine v2.1\n");
  printf("High Performance C11 Compatible Edition\n");
  printf("==============================================\n\n");

  // Generate primes
  printf("Generating %u primes...\n", MAX_PRIMES);
  if (!generate_primes_optimized()) {
    printf("Failed to generate primes\n");
    return 1;
  }
  printf("Generated %u primes\n\n", g_prime_count);

  // Build hash table
  build_hash_table();

  // Process all primes
  if (!process_million_primes_fast()) {
    printf("Processing failed\n");
    return 1;
  }

  printf("\n");

  return 0;
}