Re: Forprime and posix multithreading

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <pthread.h>
#include <errno.h>
#include <getopt.h>

#ifdef _WIN32
#include <windows.h>
#define strtok_r strtok_s  // Windows equivalent of strtok_r
#else
#include <sys/sysinfo.h>
#endif

// PARI includes (with workaround for 'long' definition)
#ifdef long
#undef long
#endif
#include <pari/pari.h>

#define VERSION "1.0.8"
#define MAX_LINE_LENGTH 256
#define MAX_FILENAME_LENGTH 100
#define DEFAULT_STACK (4UL * 1024UL * 1024UL * 1024UL) // 4GB
#define THREAD_STACK (64UL * 1024UL * 1024UL) // 64MB
#define BUFFER_SIZE (4UL * 1024UL * 1024UL) // 4MB

// Global debug flag
int debug_mode = 1;

// Windows-compatible crash handler
void crash_handler(int sig) {
    fprintf(stderr, "\n=== CRASH DETECTED ===\n");
    fprintf(stderr, "Error: signal %d\n", sig);
    exit(1);
}

// Debug print macro
#define DEBUG_PRINT(fmt, ...) \
    do { if (debug_mode) fprintf(stderr, "DEBUG: %s:%d:%s(): " fmt, \
         __FILE__, __LINE__, __func__, ##__VA_ARGS__); } while (0)

typedef struct {
    GEN n1;
    GEN n2;
    char filename[MAX_FILENAME_LENGTH];
    pthread_mutex_t* file_mutex;
    size_t stack_size;
    int thread_id;
    FILE* output_file;
} thread_args_t;

static struct option options[] = {
    {"file", required_argument, NULL, 'f'},
    {"pari-size", optional_argument, NULL, 'p'},
    {"pari-max", optional_argument, NULL, 'P'},
    {"threads", optional_argument, NULL, 't'},
    {"thread-stack", optional_argument, NULL, 's'},
    {"debug", no_argument, NULL, 'd'},
    {"help", no_argument, NULL, 'h'},
    {"version", no_argument, NULL, 'v'},
    {NULL, 0, NULL, 0}
};

int get_cpu_cores() {
#ifdef _WIN32
    SYSTEM_INFO sysinfo;
    GetSystemInfo(&sysinfo);
    return (int)sysinfo.dwNumberOfProcessors;
#else
    int cores = (int)sysconf(_SC_NPROCESSORS_ONLN);
    return cores > 0 ? cores : 1;
#endif
}

void* prime_thread(void* arg) {
     thread_args_t* args = (thread_args_t*)arg;
    pthread_mutex_lock(args->file_mutex);
    FILE* output_file = fopen(args->filename, "a");
    pthread_mutex_unlock(args->file_mutex);
    
    if (!output_file) {
        perror("Failed to open output file");
        return NULL;
    }

    // Initialize PARI in this thread
    pari_sp ltop = avma;
    forprime_t iter;
    forprime_init(&iter, args->n1, args->n2);

    char prime_buffer[1024]; // Static buffer for each prime
    GEN p;
    int count = 0;

    while ((p = forprime_next(&iter))) {
        // Convert prime to string using a fresh stack frame
        pari_sp btop = avma;
        char* prime_str = GENtostr(p);
        if (!prime_str) {
            fprintf(stderr, "Prime conversion failed\n");
            break;
        }

        // Copy to our local buffer immediately
        strncpy(prime_buffer, prime_str, sizeof(prime_buffer)-1);
        prime_buffer[sizeof(prime_buffer)-1] = '\0';
        avma = btop; // Free PARI string immediately

        // Write to file
        pthread_mutex_lock(args->file_mutex);
        fprintf(output_file, "%s\n", prime_buffer);
        if (++count % 100 == 0) {
            fflush(output_file); // Periodic flush
        }
        pthread_mutex_unlock(args->file_mutex);
    }

    fclose(output_file);
    avma = ltop; // Reset PARI stack
    return NULL;
}

void parallel_primes(GEN n1, GEN n2, const char* filename, 
                   int threads, pthread_mutex_t* mutex, size_t stack_size) {


    DEBUG_PRINT("Starting parallel prime generation (%d threads)\n", threads);

    pthread_t* workers = malloc(threads * sizeof(pthread_t));
    thread_args_t* args = malloc(threads * sizeof(thread_args_t));
    
    if (!workers || !args) {
        fprintf(stderr, "Memory allocation failed for %d threads\n", threads);
        free(workers);
        free(args);
        return;
    }

    GEN range = gsub(n2, n1);
    GEN step = gdivent(range, stoi(threads));
    GEN current = n1;


    for (int i = 0; i < threads; i++) {
        args[i].n1 = gcopy(current);
        args[i].n2 = (i == threads-1) ? gcopy(n2) : gadd(gcopy(current), step);
        strncpy(args[i].filename, filename, MAX_FILENAME_LENGTH-1);
        args[i].filename[MAX_FILENAME_LENGTH-1] = '\0';
        args[i].file_mutex = mutex;
        args[i].stack_size = stack_size;
        args[i].thread_id = i;

        pthread_mutex_lock(args->file_mutex);
        args->output_file = fopen(args->filename, "a");
        pthread_mutex_unlock(args->file_mutex);

        DEBUG_PRINT("Creating thread %d for range ", i);
        if (debug_mode) {
            pari_fprintf(stderr, "%Ps to %Ps\n", args[i].n1, args[i].n2);
        }

        if (pthread_create(&workers[i], NULL, prime_thread, &args[i])) {
            perror("Thread creation failed");
            break;
        }

        current = gadd(args[i].n2, gen_1);
    }

    for (int i = 0; i < threads; i++) {
        if (workers[i]) {
            DEBUG_PRINT("Waiting for thread %d\n", i);
            pthread_join(workers[i], NULL);
        }
    }

    free(workers);
    free(args);
    DEBUG_PRINT("Parallel prime generation completed\n");
}

int main(int argc, char** argv) {
    // Set up crash handler
    signal(SIGSEGV, crash_handler);
    signal(SIGABRT, crash_handler);

    size_t pari_stack = DEFAULT_STACK;
    size_t max_stack = 2 * DEFAULT_STACK;
    int threads = get_cpu_cores();
    size_t thread_stack = THREAD_STACK;
    FILE* input = NULL;
    pthread_mutex_t file_mutex;

    // Initialize mutex with error checking
    if (pthread_mutex_init(&file_mutex, NULL)) {
        perror("Mutex init failed");
        return 1;
    }

    // Parse command line options
    int opt;
    while ((opt = getopt_long(argc, argv, "f:p:P:t:s:dhv", options, NULL)) != -1) {
        switch (opt) {
            case 'f': 
                input = fopen(optarg, "r");
                if (!input) {
                    fprintf(stderr, "Failed to open input file '%s': %s\n", 
                            optarg, strerror(errno));
                    pthread_mutex_destroy(&file_mutex);
                    return 1;
                }
                break;
            case 'p': 
                pari_stack = (size_t)atol(optarg) * 1024UL * 1024UL;
                break;
            case 'P':
                max_stack = (size_t)atol(optarg) * 1024UL * 1024UL;
                break;
            case 't':
                threads = atoi(optarg);
                if (threads < 1) {
                    fprintf(stderr, "Invalid thread count: %d\n", threads);
                    threads = 1;
                }
                break;
            case 's':
                thread_stack = (size_t)atol(optarg) * 1024UL * 1024UL;
                break;
            case 'd':
                debug_mode = 1;
                break;
            case 'h':
                printf("Usage: %s -f file.csv [options]\n", argv[0]);
                printf("Options:\n");
                printf("  -f, --file FILE          Input CSV file\n");
                printf("  -p, --pari-size SIZE     Initial PARI stack in MB\n");
                printf("  -P, --pari-max SIZE      Max PARI stack in MB\n");
                printf("  -t, --threads N          Number of threads\n");
                printf("  -s, --thread-stack SIZE  Thread stack in MB\n");
                printf("  -d, --debug              Enable debug output\n");
                printf("  -h, --help               Show this help\n");
                printf("  -v, --version            Show version\n");
                pthread_mutex_destroy(&file_mutex);
                return 0;
            case 'v':
                printf("Version: %s\n", VERSION);
                pthread_mutex_destroy(&file_mutex);
                return 0;
            default:
                fprintf(stderr, "Invalid option\n");
                pthread_mutex_destroy(&file_mutex);
                return 1;
        }
    }

    if (!input) {
        fprintf(stderr, "No input file specified (use -f option)\n");
        pthread_mutex_destroy(&file_mutex);
        return 1;
    }

    DEBUG_PRINT("Initializing PARI with stack %zu MB, max %zu MB\n", 
               pari_stack/(1024*1024), max_stack/(1024*1024));
    pari_init(pari_stack, max_stack);

    char line[MAX_LINE_LENGTH];
    char* saveptr;  // For strtok_r
    while (fgets(line, sizeof(line), input)) {
        pari_sp ltop = avma;
        line[strcspn(line, "\n")] = 0;
        DEBUG_PRINT("Processing line: %s\n", line);

        int n1, n2, n3, mult;
        char filename[MAX_FILENAME_LENGTH];
        char* token = strtok_r(line, ",", &saveptr);
        if (!token) {
            fprintf(stderr, "Skipping malformed line: %s\n", line);
            continue;
        }
        n1 = atoi(token);

        token = strtok_r(NULL, ",", &saveptr);
        if (!token) {
            fprintf(stderr, "Missing n2 in line: %s\n", line);
            continue;
        }
        n2 = atoi(token);

        token = strtok_r(NULL, ",", &saveptr);
        if (!token) {
            fprintf(stderr, "Missing multiplier in line: %s\n", line);
            continue;
        }
        mult = atoi(token);

        token = strtok_r(NULL, ",", &saveptr);
        if (!token) {
            fprintf(stderr, "Missing n3 in line: %s\n", line);
            continue;
        }
        n3 = atoi(token);

        token = strtok_r(NULL, ",", &saveptr);
        if (!token) {
            fprintf(stderr, "Missing filename in line: %s\n", line);
            continue;
        }
        strncpy(filename, token, MAX_FILENAME_LENGTH-1);
        filename[MAX_FILENAME_LENGTH-1] = '\0';

        DEBUG_PRINT("Generating primes for 10^%d to 10^%d + %d*10^%d -> %s\n",
                   n1, n2, mult, n3, filename);

        GEN base = stoi(10);
        GEN range_start = gpow(base, stoi(n1), -1);
        GEN part1 = gpow(base, stoi(n2), -1);
        GEN part2 = gmul(stoi(mult), gpow(base, stoi(n3), -1));
        GEN range_end = gadd(part1, part2);

                // In main(), after calculating range_start and range_end:
        if (gexpo(range_start) > 50) {  // Warn about very large numbers
            fprintf(stderr, "WARNING: Calculating primes for extremely large numbers (10^%ld)\n",
                gexpo(range_start));
            fprintf(stderr, "This may take prohibitively long or exhaust memory\n");
            // Consider adding: if (gexpo(range_start) > 50) continue;
}


        parallel_primes(range_start, range_end, filename, threads, &file_mutex, thread_stack);
        avma = ltop;
    }

    fclose(input);
    pari_close();
    pthread_mutex_destroy(&file_mutex);
    return 0;
}