1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
/*
* Licensed to the Apache Software Foundation (ASF) under one
* or more contributor license agreements. See the NOTICE file
* distributed with this work for additional information
* regarding copyright ownership. The ASF licenses this file
* to you under the Apache License, Version 2.0 (the
* "License"); you may not use this file except in compliance
* with the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing,
* software distributed under the License is distributed on an
* "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
* KIND, either express or implied. See the License for the
* specific language governing permissions and limitations
* under the License.
*/
#ifndef _HARMONICNUMBERS_INTERNAL_HPP_
#define _HARMONICNUMBERS_INTERNAL_HPP_
#include "HarmonicNumbers.hpp"
#include <cmath>
namespace datasketches {
template<typename A>
double HarmonicNumbers<A>::getBitMapEstimate(const int bitVectorLength, const int numBitsSet) {
return (bitVectorLength * (harmonicNumber(bitVectorLength) - harmonicNumber(bitVectorLength - numBitsSet)));
}
static const int NUM_EXACT_HARMONIC_NUMBERS = 25;
static double tableOfExactHarmonicNumbers[] = {
0.0, // 0
1.0, // 1
1.5, // 2
11.0 / 6.0, // 3
25.0 / 12.0, // 4
137.0 / 60.0, // 5
49.0 / 20.0, // 6
363.0 / 140.0, // 7
761.0 / 280.0, // 8
7129.0 / 2520.0, // 9
7381.0 / 2520.0, // 10
83711.0 / 27720.0, // 11
86021.0 / 27720.0, // 12
1145993.0 / 360360.0, // 13
1171733.0 / 360360.0, // 14
1195757.0 / 360360.0, // 15
2436559.0 / 720720.0, // 16
42142223.0 / 12252240.0, // 17
14274301.0 / 4084080.0, // 18
275295799.0 / 77597520.0, // 19
55835135.0 / 15519504.0, // 20
18858053.0 / 5173168.0, // 21
19093197.0 / 5173168.0, // 22
444316699.0 / 118982864.0, // 23
1347822955.0 / 356948592.0 // 24
};
static const double EULER_MASCHERONI_CONSTANT = 0.577215664901532860606512090082;
template<typename A>
double HarmonicNumbers<A>::harmonicNumber(const uint64_t x_i) {
if (x_i < NUM_EXACT_HARMONIC_NUMBERS) {
return tableOfExactHarmonicNumbers[x_i];
} else {
double x = static_cast<double>(x_i);
double invSq = 1.0 / (x * x);
double sum = log(x) + EULER_MASCHERONI_CONSTANT + (1.0 / (2.0 * x));
/* note: the number of terms included from this series expansion is appropriate
for the size of the exact table (25) and the precision of doubles */
double pow = invSq; // now n^-2
sum -= pow * (1.0 / 12.0);
pow *= invSq; // now n^-4
sum += pow * (1.0 / 120.0);
pow *= invSq; /* now n^-6 */
sum -= pow * (1.0 / 252.0);
pow *= invSq; /* now n^-8 */
sum += pow * (1.0 / 240.0);
return sum;
}
}
}
#endif // _HARMONICNUMBERS_INTERNAL_HPP_