Blog Entry Tuesday, July 16, 2024, Prime Number Zeta Function Values for s = 2 and Primes < 1,000,000 (c) James Pate Williams, Jr.

/*
* PrimeZetaFunction.c (c) Monday, July 15, 2024
* by James Pate Williams, Jr.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define MAX_SIEVE 100000000
#define MAX_PRIME_INDEX 5761454
#define SIEVE_SIZE (MAX_SIEVE / BITS_PER_LONG + 1)

long prime[MAX_PRIME_INDEX], sieve[SIEVE_SIZE];

long get_bit(long i, long* sieve)
{
    long b = i % BITS_PER_LONG;
    long c = i / BITS_PER_LONG;

    return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long* sieve)
{
    long b = i % BITS_PER_LONG;
    long c = i / BITS_PER_LONG;
    long mask = 1 << (BITS_PER_LONG_1 - b);

    if (v == 1)
        sieve[c] |= mask;
    else
        sieve[c] &= ~mask;
}

void Sieve(long n, long* sieve)
{
    long c, i, inc;

    set_bit(0, 0, sieve);
    set_bit(1, 0, sieve);
    set_bit(2, 1, sieve);
    for (i = 3; i <= n; i++)
        set_bit(i, i & 1, sieve);
    c = 3;
    do {
        i = c * c, inc = c + c;
        while (i <= n) {
            set_bit(i, 0, sieve);
            i += inc;
        }
        c += 2;
        while (!get_bit(c, sieve)) c++;
    } while (c * c <= n);
}

long double primeZetaFunction(long N, long s)
{
    long n, p;
    long double sum = 0.0L;

    for (n = N; n >= 0; n--)
    {
        p = prime[n];
        sum += 1.0 / powl((long double)p, (long double)s);
    }
        
    return sum;
}

int main()
{
    long N = 0, i = 0, p = 2, s = 0;
    double runtime = 0.0;
    clock_t time0 = clock(), time1 = 0;

    Sieve(MAX_SIEVE, sieve);

    for (i = 0; i <= MAX_PRIME_INDEX; i++) {
        while (!get_bit(p, sieve)) p++;
        prime[i] = p++;
    }

    time1 = clock();
    runtime = ((double)time1 - time0) / CLOCKS_PER_SEC;
    printf_s("Runtime in seconds to build sieve: %Lf\n", runtime);
    for (;;) {
        printf_s("N = <= %ld ", 1000000);
        scanf_s("%ld", &N);

        if (N == 0)
            break;

        printf_s("s > 1 = ");
        scanf_s("%ld", &s);
        printf_s("prime zeta = %16.15Lf\n", primeZetaFunction(N, s));
    }

    return 0;
}