Skip to content

Commit 8635a41

Browse files
leerhoclaude
andcommitted
Finish repairing the deferred KxQ rebuild, and pin the estimator constants
Three independent defects in the HLL estimator state. 1. The deferred KxQ rebuild from #364 still leaves curMin/numAtCurMin merge order dependent. #512 repaired the visible damage, but not the representation the rebuild writes: check_rebuild_kxq_cur_min() stores the true minimum register value and the count at that minimum, while the rest of the HLL_8 code maintains curMin == 0 with numAtCurMin counting the zero registers. When the merged array has no zero register the rebuild leaves curMin > 0, and from then on numAtCurMin -= (curVal == 0) never fires, so the stored pair freezes and drifts away from the registers. Its value depends on when the rebuild fired, hence on merge order. Over 4800 randomized merge-order comparisons 92 differed only in these two fields; the same sweep before #364 gives 0. Estimates and bounds are unaffected because both consumers branch on curMin == 0, so this is a serialization determinism defect. Emit the canonical form instead, and add the HLL_8 / HLL-mode guard that datasketches-java has, so the rebuild can never rewrite curMin on an HLL_4 array whose nibbles are stored relative to it. 2. HLL_HIP_RSE_FACTOR and HLL_NON_HIP_RSE_FACTOR were rounded to seven digits. For lg_k > 12 the bounds use the closed form rather than the interpolation table, so this is observable: the HIP bounds differ from Java by ~4.3e-11 relative and the non-HIP factor is off by 1.7e-6 relative, skewing the bounds of every union result above lg_k 12. Use full precision literals taken from Java's Double.toString output, rather than computing sqrt(log(2.0)) at runtime, so the value cannot vary with the platform libm. 3. harmonicNumber() and getHllBitMapEstimate() called the platform std::log. getBitMapEstimate computes K * (H(K) - H(K - numHit)), whose cancellation amplifies a 1 ULP log difference by about 17x. Apple libm differs from fdlibm on 0.62% of a 400000 sample sweep. Java's StrictMath.log is specified to be fdlibm, so add common/include/fdlibm_log.hpp (FDLIBM 5.3 __ieee754_log, Sun notice preserved) and use it in both places. fdlibm needs strict IEEE evaluation, and the clang pragma is overridden by an explicit -ffp-contract=fast, so the contraction sensitive expressions go through a volatile round trip; verified bit-identical to Java under -O2, -O3, -ffp-contract=on and =fast. This changes serialized bytes, unlike the preceding PR. Union results carry the canonical curMin/numAtCurMin, linear-counting estimates shift ~1e-15, and bounds above lg_k 12 shift as described. Reading is unaffected: older images still deserialize and yield identical estimates and bounds. Adds HllKxqRebuildTest.cpp. Three of the four assertions that also compile against the parent branch fail there and pass here; the fourth already held and is kept as a guard. With the companion datasketches-java change, a 1041 record corpus spanning lg_k 4..21, all three target types, 17 sizes across LIST/SET/HLL, round trips and 80 union scenarios goes from 428 differing records to 0, comparing every serialized byte, estimate, composite estimate and bound as raw IEEE bits. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
1 parent 637d748 commit 8635a41

7 files changed

Lines changed: 274 additions & 16 deletions

File tree

‎common/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ install(FILES
4040
include/conditional_back_inserter.hpp
4141
include/conditional_forward.hpp
4242
include/count_zeros.hpp
43+
include/fdlibm_log.hpp
4344
include/inv_pow2_table.hpp
4445
include/kolmogorov_smirnov_impl.hpp
4546
include/kolmogorov_smirnov.hpp

‎common/include/fdlibm_log.hpp‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
// fdlibm __ieee754_log, used by Java's StrictMath.log (and by Math.log on most JVMs).
2+
// Derived from FDLIBM 5.3, Copyright (C) 1993 by Sun Microsystems, Inc.
3+
// "Permission to use, copy, modify, and distribute this software is freely granted,
4+
// provided that this notice is preserved."
5+
#ifndef FDLIBM_LOG_HPP
6+
#define FDLIBM_LOG_HPP
7+
#include <cstdint>
8+
#include <cstring>
9+
#include <limits>
10+
11+
namespace fdlibm {
12+
13+
inline int32_t hi_word(double x){ uint64_t u; std::memcpy(&u,&x,8); return (int32_t)(uint32_t)(u>>32); }
14+
inline uint32_t lo_word(double x){ uint64_t u; std::memcpy(&u,&x,8); return (uint32_t)u; }
15+
inline void set_hi_word(double& x, uint32_t hi){ uint64_t u; std::memcpy(&u,&x,8);
16+
u = (u & 0x00000000ffffffffULL) | ((uint64_t)hi<<32); std::memcpy(&x,&u,8); }
17+
18+
// Forces a value to be rounded to a double before it is used again. fdlibm needs strict
19+
// IEEE-754 evaluation: a fused multiply-add anywhere in the polynomial below changes the
20+
// result. Pragmas are overridden by an explicit -ffp-contract=fast, so use a volatile
21+
// round-trip, which the standard requires the compiler to honour.
22+
inline double rnd(double v) { volatile double t = v; return t; }
23+
24+
inline double log(double x) {
25+
// fdlibm depends on strict IEEE-754 evaluation: a fused multiply-add would change the result
26+
// of the polynomial evaluation below, so contraction must be off for this function.
27+
#if defined(__clang__)
28+
#pragma clang fp contract(off)
29+
#endif
30+
static const double
31+
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
32+
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
33+
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
34+
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
35+
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
36+
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
37+
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
38+
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
39+
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
40+
Lg7 = 1.479819860511658591e-01, /* 3FC2F112 DF3E5244 */
41+
zero = 0.0;
42+
43+
double hfsq,f,s,z,R,w,t1,t2,dk;
44+
int32_t k,hx,i,j;
45+
uint32_t lx;
46+
47+
hx = hi_word(x);
48+
lx = lo_word(x);
49+
50+
k = 0;
51+
if (hx < 0x00100000) { /* x < 2**-1022 */
52+
// fdlibm writes these as -two54/zero and (x-x)/zero, which also raise the divide-by-zero
53+
// and invalid flags. MSVC rejects a compile-time division by a zero constant (C2124), so
54+
// return the same values directly. The estimators never call log() with these inputs.
55+
if (((hx & 0x7fffffff) | lx) == 0) { /* log(+-0) = -inf */
56+
return -std::numeric_limits<double>::infinity();
57+
}
58+
if (hx < 0) { /* log(-#) = NaN */
59+
return std::numeric_limits<double>::quiet_NaN();
60+
}
61+
k -= 54; x *= two54; /* subnormal: scale up */
62+
hx = hi_word(x);
63+
}
64+
if (hx >= 0x7ff00000) { return x+x; }
65+
k += (hx>>20) - 1023;
66+
hx &= 0x000fffff;
67+
i = (hx + 0x95f64) & 0x100000;
68+
set_hi_word(x, (uint32_t)(hx | (i ^ 0x3ff00000))); /* normalize x or x/2 */
69+
k += (i>>20);
70+
f = x - 1.0;
71+
if ((0x000fffff & (2+hx)) < 3) { /* |f| < 2**-20 */
72+
if (f == zero) {
73+
if (k == 0) { return zero; }
74+
dk = (double)k; return rnd(dk*ln2_hi) + rnd(dk*ln2_lo);
75+
}
76+
R = rnd(rnd(f*f)*rnd(0.5 - rnd(0.33333333333333333*f)));
77+
if (k == 0) { return f-R; }
78+
dk = (double)k; return rnd(dk*ln2_hi) - (rnd(R - rnd(dk*ln2_lo)) - f);
79+
}
80+
s = f/(2.0+f);
81+
dk = (double)k;
82+
z = s*s;
83+
i = hx - 0x6147a;
84+
w = z*z;
85+
j = 0x6b851 - hx;
86+
t1 = rnd(w*rnd(Lg2 + rnd(w*rnd(Lg4 + rnd(w*Lg6)))));
87+
t2 = rnd(z*rnd(Lg1 + rnd(w*rnd(Lg3 + rnd(w*rnd(Lg5 + rnd(w*Lg7)))))));
88+
i |= j;
89+
R = t2 + t1;
90+
if (i > 0) {
91+
hfsq = rnd(0.5*f)*f;
92+
if (k == 0) { return f - rnd(hfsq - rnd(s*(hfsq+R))); }
93+
return rnd(dk*ln2_hi) - (rnd(hfsq - rnd(rnd(s*(hfsq+R)) + rnd(dk*ln2_lo))) - f);
94+
} else {
95+
if (k == 0) { return f - rnd(s*(f-R)); }
96+
return rnd(dk*ln2_hi) - (rnd(rnd(s*(f-R)) - rnd(dk*ln2_lo)) - f);
97+
}
98+
}
99+
100+
} // namespace fdlibm
101+
#endif

‎hll/include/HarmonicNumbers-internal.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
#define _HARMONICNUMBERS_INTERNAL_HPP_
2222

2323
#include "HarmonicNumbers.hpp"
24+
#include "fdlibm_log.hpp"
2425

2526
#include <cmath>
2627

@@ -70,7 +71,7 @@ double HarmonicNumbers<A>::harmonicNumber(const uint64_t x_i) {
7071
} else {
7172
double x = static_cast<double>(x_i);
7273
double invSq = 1.0 / (x * x);
73-
double sum = log(x) + EULER_MASCHERONI_CONSTANT + (1.0 / (2.0 * x));
74+
double sum = fdlibm::log(x) + EULER_MASCHERONI_CONSTANT + (1.0 / (2.0 * x));
7475
/* note: the number of terms included from this series expansion is appropriate
7576
for the size of the exact table (25) and the precision of doubles */
7677
double pow = invSq; // now n^-2

‎hll/include/HllArray-internal.hpp‎

Lines changed: 17 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
#define _HLLARRAY_INTERNAL_HPP_
2222

2323
#include "HllArray.hpp"
24+
#include "fdlibm_log.hpp"
2425
#include "HllUtil.hpp"
2526
#include "HarmonicNumbers.hpp"
2627
#include "CubicInterpolation.hpp"
@@ -568,7 +569,7 @@ double HllArray<A>::getHllBitMapEstimate() const {
568569

569570
//This will eventually go away.
570571
if (numUnhitBuckets == 0) {
571-
return configK * log(configK / 0.5);
572+
return configK * fdlibm::log(configK / 0.5);
572573
}
573574

574575
const uint32_t numHitBuckets = configK - numUnhitBuckets;
@@ -601,34 +602,37 @@ bool HllArray<A>::isRebuildKxqCurminFlag() const {
601602
template<typename A>
602603
void HllArray<A>::check_rebuild_kxq_cur_min() {
603604
if (!rebuild_kxq_curmin_) { return; }
605+
// the deferred rebuild is only ever set on an HLL_8 union gadget. Guarding here also keeps
606+
// this from ever rewriting curMin_ on an HLL_4 array, whose nibbles are stored relative to it.
607+
if (this->getCurMode() != hll_mode::HLL || this->getTgtHllType() != target_hll_type::HLL_8) {
608+
rebuild_kxq_curmin_ = false;
609+
return;
610+
}
604611

605-
uint8_t cur_min = 64;
606-
uint32_t num_at_cur_min = 0;
612+
uint32_t num_zeros = 0;
607613
double kxq0 = 1 << this->lgConfigK_;
608614
double kxq1 = 0;
609615

610-
auto it = this->begin(true); // want all points to adjust cur_min
616+
auto it = this->begin(true); // want all slots, including the empty ones
611617
const auto end = this->end();
612618
while (it != end) {
613619
uint8_t v = HllUtil<A>::getValue(*it);
614620
if (v > 0) {
615621
if (v < 32) { kxq0 += INVERSE_POWERS_OF_2[v] - 1.0; }
616622
else { kxq1 += INVERSE_POWERS_OF_2[v] - 1.0; }
617-
}
618-
if (v > cur_min) { ++it; continue; }
619-
if (v < cur_min) {
620-
cur_min = v;
621-
num_at_cur_min = 1;
622623
} else {
623-
++num_at_cur_min;
624-
}
624+
++num_zeros;
625+
}
625626
++it;
626627
}
627628

628629
kxq0_ = kxq0;
629630
kxq1_ = kxq1;
630-
curMin_ = cur_min;
631-
numAtCurMin_ = num_at_cur_min;
631+
// HLL_8 convention: curMin is always 0 and numAtCurMin is the number of zero registers.
632+
// That is the representation the incremental update path maintains, so the rebuilt state is
633+
// indistinguishable from it and the timing of this rebuild is not observable in the image.
634+
curMin_ = 0;
635+
numAtCurMin_ = num_zeros;
632636
rebuild_kxq_curmin_ = false;
633637
// HipAccum is not affected
634638

‎hll/include/HllUtil.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -87,8 +87,8 @@ static const uint32_t EMPTY = 0;
8787
static const uint8_t MIN_LOG_K = 4;
8888
static const uint8_t MAX_LOG_K = 21;
8989

90-
static const double HLL_HIP_RSE_FACTOR = 0.8325546; // sqrt(ln(2))
91-
static const double HLL_NON_HIP_RSE_FACTOR = 1.03896; // sqrt((3 * ln(2)) - 1)
90+
static const double HLL_HIP_RSE_FACTOR = 0.8325546111576977; // sqrt(ln(2))
91+
static const double HLL_NON_HIP_RSE_FACTOR = 1.0389617614136892; // sqrt((3 * ln(2)) - 1)
9292
static const double COUPON_RSE_FACTOR = 0.409; // at transition point not the asymptote
9393
static const double COUPON_RSE = COUPON_RSE_FACTOR / (1 << 13);
9494

‎hll/test/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ target_sources(hll_test
4343
CrossCountingTest.cpp
4444
HllArrayTest.cpp
4545
HllFullSizeTest.cpp
46+
HllKxqRebuildTest.cpp
4647
HllSketchTest.cpp
4748
HllUnionTest.cpp
4849
TablesTest.cpp

‎hll/test/HllKxqRebuildTest.cpp‎

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
/*
2+
* Licensed to the Apache Software Foundation (ASF) under one
3+
* or more contributor license agreements. See the NOTICE file
4+
* distributed with this work for additional information
5+
* regarding copyright ownership. The ASF licenses this file
6+
* to you under the Apache License, Version 2.0 (the
7+
* "License"); you may not use this file except in compliance
8+
* with the License. You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing,
13+
* software distributed under the License is distributed on an
14+
* "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
15+
* KIND, either express or implied. See the License for the
16+
* specific language governing permissions and limitations
17+
* under the License.
18+
*/
19+
20+
#include <catch2/catch.hpp>
21+
#include <cmath>
22+
#include <cstring>
23+
#include <limits>
24+
25+
#include "hll.hpp"
26+
#include "fdlibm_log.hpp"
27+
28+
namespace datasketches {
29+
30+
// offsets into the HLL updatable image
31+
static const size_t CUR_MIN_BYTE = 6;
32+
static const size_t NUM_AT_CUR_MIN_INT = 32;
33+
static const size_t HLL_BYTE_ARR_START = 40;
34+
35+
static hll_sketch make(uint8_t lg_k, target_hll_type type, uint64_t lo, uint64_t hi) {
36+
hll_sketch sk(lg_k, type);
37+
for (uint64_t i = lo; i < hi; ++i) sk.update(i);
38+
return sk;
39+
}
40+
static uint32_t num_at_cur_min_of(const hll_sketch::vector_bytes& img) {
41+
uint32_t v; std::memcpy(&v, img.data() + NUM_AT_CUR_MIN_INT, sizeof(v)); return v;
42+
}
43+
44+
TEST_CASE("hll kxq rebuild: union result is merge-order independent", "[hll_kxq]") {
45+
// C stays in SET mode, so this exercises the coupon-update path into a gadget with a
46+
// deferred rebuild pending, which is where the stored curMin/numAtCurMin used to drift
47+
const hll_sketch a = make(12, HLL_4, 20000, 30364);
48+
const hll_sketch b = make(10, HLL_8, 5000, 14699);
49+
const hll_sketch c = make(17, HLL_4, 70000, 72598);
50+
const hll_sketch* in[3] = {&a, &b, &c};
51+
52+
auto run = [&](int i, int j, int k) {
53+
hll_union u(7);
54+
u.update(*in[i]); u.update(*in[j]); u.update(*in[k]);
55+
return u.get_result(HLL_8).serialize_updatable();
56+
};
57+
const auto ref = run(0, 1, 2);
58+
REQUIRE(run(0, 2, 1) == ref);
59+
REQUIRE(run(1, 0, 2) == ref);
60+
REQUIRE(run(1, 2, 0) == ref);
61+
REQUIRE(run(2, 0, 1) == ref);
62+
REQUIRE(run(2, 1, 0) == ref);
63+
}
64+
65+
TEST_CASE("hll kxq rebuild: reading an estimate does not change a later result",
66+
"[hll_kxq]") {
67+
// the rebuild is lazy; when it fires must not be observable in the serialized image
68+
const hll_sketch p = make(13, HLL_8, 0, 50000);
69+
const hll_sketch q = make(13, HLL_8, 50000, 100000);
70+
71+
for (uint8_t lg_max_k: {uint8_t(7), uint8_t(8), uint8_t(9)}) {
72+
hll_union peeked(lg_max_k);
73+
peeked.update(p); peeked.update(q);
74+
(void) peeked.get_estimate(); // forces the rebuild here
75+
for (uint64_t v = 9000000; v < 9400000; ++v) peeked.update(v);
76+
77+
hll_union plain(lg_max_k);
78+
plain.update(p); plain.update(q);
79+
for (uint64_t v = 9000000; v < 9400000; ++v) plain.update(v);
80+
81+
REQUIRE(peeked.get_result(HLL_8).serialize_updatable()
82+
== plain.get_result(HLL_8).serialize_updatable());
83+
}
84+
}
85+
86+
TEST_CASE("hll kxq rebuild: stored curMin and numAtCurMin agree with the registers",
87+
"[hll_kxq]") {
88+
// HLL_8 convention: curMin is always 0 and numAtCurMin counts the zero registers
89+
const hll_sketch a = make(15, HLL_8, 0, 100000);
90+
const hll_sketch b = make(8, HLL_8, 100000, 200000);
91+
hll_union u(8);
92+
u.update(a); u.update(b);
93+
const auto img = u.get_result(HLL_8).serialize_updatable();
94+
95+
uint32_t zeros = 0;
96+
for (size_t i = HLL_BYTE_ARR_START; i < img.size(); ++i) if (img[i] == 0) ++zeros;
97+
98+
REQUIRE(img[CUR_MIN_BYTE] == 0);
99+
REQUIRE(num_at_cur_min_of(img) == zeros);
100+
}
101+
102+
TEST_CASE("hll kxq rebuild: relative error constants are full precision", "[hll_kxq]") {
103+
// lg_k > 12 uses the closed form rather than the interpolation table, so a constant
104+
// truncated to seven digits is directly observable in the bounds
105+
const double hip = std::sqrt(std::log(2.0)); // sqrt(ln 2)
106+
const double non_hip = std::sqrt((3.0 * std::log(2.0)) - 1.0); // sqrt(3 ln 2 - 1)
107+
108+
for (uint8_t lg_k: {uint8_t(13), uint8_t(16), uint8_t(21)}) {
109+
const double k = static_cast<double>(1 << lg_k);
110+
for (uint8_t sd = 1; sd <= 3; ++sd) {
111+
const double got_hip = hll_union::get_rel_err(false, false, lg_k, sd);
112+
const double got_non = hll_union::get_rel_err(false, true, lg_k, sd);
113+
REQUIRE(got_hip == Approx(sd * hip / std::sqrt(k)).epsilon(1e-15));
114+
REQUIRE(got_non == Approx(sd * non_hip / std::sqrt(k)).epsilon(1e-15));
115+
}
116+
}
117+
}
118+
119+
TEST_CASE("hll kxq rebuild: log matches the fdlibm reference", "[hll_kxq]") {
120+
// the linear counting estimator subtracts two nearby harmonic numbers, which amplifies a
121+
// 1 ULP difference in log() by more than an order of magnitude. Pin log() to fdlibm, the
122+
// function datasketches-java's StrictMath.log is specified to be. These expected values
123+
// were taken from Java; the platform libm differs from every one of them.
124+
struct { int x; uint64_t bits; } expected[] = {
125+
{ 48, 0x400ef8383c50bb74ULL},
126+
{ 74, 0x4011375cd6fcab1cULL},
127+
{ 185, 0x4014e1a4f518c72cULL},
128+
{ 196, 0x40151cca16d7bba8ULL},
129+
{ 299, 0x4016cd411481a020ULL},
130+
{ 308, 0x4016eb9f470ac0b8ULL},
131+
{ 334, 0x40173e9bbe951e9cULL},
132+
{ 343, 0x401759d602a5c3c2ULL},
133+
{ 1261, 0x401c8f031e7e1220ULL},
134+
};
135+
for (const auto& e: expected) {
136+
const double got = fdlibm::log(static_cast<double>(e.x));
137+
uint64_t bits; std::memcpy(&bits, &got, sizeof(bits));
138+
REQUIRE(bits == e.bits);
139+
}
140+
141+
// edge cases: fdlibm produces these by dividing by a zero constant, which MSVC rejects at
142+
// compile time, so they are returned directly. Pin the values.
143+
REQUIRE(fdlibm::log(0.0) == -std::numeric_limits<double>::infinity());
144+
REQUIRE(fdlibm::log(-0.0) == -std::numeric_limits<double>::infinity());
145+
REQUIRE(std::isnan(fdlibm::log(-1.0)));
146+
REQUIRE(fdlibm::log(1.0) == 0.0);
147+
REQUIRE(fdlibm::log(5e-320) == Approx(-735.2178).epsilon(1e-6)); // subnormal scaling path
148+
}
149+
150+
} /* namespace datasketches */

0 commit comments

Comments
 (0)