/*
 * 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_