Sum of square roots の近似関数

以下は sum of square roots (√1 + √2 + … + √N) を計算する関数です。 f32 か f64 ならほぼ正確ですが、任意精度実装と比較したときの相対誤差がマシンイプシロンを超える箇所がいくらかあります。 x^2 mod 1.0 の逆微分を計算するときに使えます。

#include <array>
#include <cmath>
#include <cstdint>
#include <type_traits>

// Only works up to N = 2^(significand_bits + 1).
template<typename Real> inline Real sum_of_sqrt(uint64_t N)
{
  alignas(64) static constexpr std::array<Real, 64> table{
    Real(0.00000000000000000e00), Real(1.00000000000000000e00),
    Real(2.41421356237309492e00), Real(4.14626436994197256e00),
    Real(6.14626436994197256e00), Real(8.38233234744176237e00),
    Real(1.08318220902249394e01), Real(1.34775734012895310e01),
    Real(1.63060005260357208e01), Real(1.93060005260357208e01),
    Real(2.24682781862040990e01), Real(2.57849029765595006e01),
    Real(2.92490045916972541e01), Real(3.28545558671612454e01),
    Real(3.65962132539351828e01), Real(4.04691966001426024e01),
    Real(4.44691966001426024e01), Real(4.85923022257602639e01),
    Real(5.28349429128795478e01), Real(5.71938418564202209e01),
    Real(6.16659778114198005e01), Real(6.62485535063756430e01),
    Real(7.09389692661990665e01), Real(7.57348007895117945e01),
    Real(8.06337802750781520e01), Real(8.56337802750781520e01),
    Real(9.07327997886709312e01), Real(9.59289522113775632e01),
    Real(1.01220454833506750e02), Real(1.06605619640641251e02),
    Real(1.12082845215692913e02), Real(1.17650609578522932e02),
    Real(1.23307463828015315e02), Real(1.29052026474553344e02),
    Real(1.34882978369398643e02), Real(1.40799058152498247e02),
    Real(1.46799058152498247e02), Real(1.52881820682796473e02),
    Real(1.59046234685765455e02), Real(1.65291232684163845e02),
    Real(1.71615788004500615e02), Real(1.78018912241933464e02),
    Real(1.84499652940341321e02), Real(1.91057091464643321e02),
    Real(1.97690341045354131e02), Real(2.04398544977853476e02),
    Real(2.11180874960978770e02), Real(2.18036529561379808e02),
    Real(2.24964732791655308e02), Real(2.31964732791655308e02),
    Real(2.39035800603520784e02), Real(2.46177229032063622e02),
    Real(2.53388331582991611e02), Real(2.60668441472272150e02),
    Real(2.68016910700621679e02), Real(2.75433109187717321e02),
    Real(2.82916423961265195e02), Real(2.90466258396535977e02),
    Real(2.98082031502399843e02), Real(3.05763177250268484e02),
    Real(3.13509143942683295e02), Real(3.21319393618589970e02),
    Real(3.29193401492601765e02), Real(3.37130655425795567e02),
  };
  if (N < static_cast<uint64_t>(table.size())) return table[N];

  Real n = static_cast<Real>(N);
  Real sqrt_n = std::sqrt(n);

  constexpr Real zeta_term = Real(-0.207886224977354566); // zeta(-0.5).
  Real base = ((Real(2) / Real(3)) * n + Real(0.5)) * sqrt_n + zeta_term;

  constexpr Real C3 = Real(1) / Real(9216);
  constexpr Real C2 = Real(-1) / Real(1920);
  constexpr Real C1 = Real(1) / Real(24);
  if constexpr (std::is_same_v<Real, float>) {
    return base + C1 / sqrt_n;
  } else {
    Real m = Real(1) / (n * n);
    return base + (((C3 * m + C2) * m + C1) / sqrt_n);
  }
}

参考文献: Ramanujan, Srinivasa. “On the sum of the square roots of the first n natural numbers.” J. Indian Math. Soc 7 (1915): 173-175.