Reciprocal Fibonacci constant and Dottie number. (boostorg#407)

doc/constants/constants.qbk

@@ -382,6 +382,8 @@ This section lists the mathematical constants, their use(s) (and sometimes ratio
 [[first_feigenbaum] [] [4.6692016] [[@https://en.wikipedia.org/wiki/Feigenbaum_constants First Feigenbaum constant]] ]
 [[plastic] [Real solution of x[super 3] = x + 1] [1.324717957] [[@https://en.wikipedia.org/wiki/Plastic_number Plastic constant]] ]
 [[gauss] [Reciprocal of agm(1, [radic]2)] [0.8346268] [[@https://en.wikipedia.org/wiki/Gauss%27s_constant Gauss's constant]] ]
+[[dottie] [Solution of cos(x) = x] [0.739085] [[@https://en.wikipedia.org/wiki/Dottie_number Dottie's number]] ]
+[[reciprocal_fibonacci] [Sum of reciprocals of Fibonacci numbers] [3.359885666] [[@https://en.wikipedia.org/wiki/Reciprocal_Fibonacci_constant Reciprocal Fibonacci constant]] ]
 ] [/table]
 
 

include/boost/math/constants/calculate_constants.hpp

@@ -1037,6 +1037,45 @@ inline T constant_gauss<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((boos
    return 2/(a + g);
 }
 
+template <class T>
+template<int N>
+inline T constant_dottie<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((boost::integral_constant<int, N>)))
+{
+  // Error analysis: cos(x(1+d)) - x(1+d) = -(sin(x)+1)xd; plug in x = 0.739 gives -1.236d; take d as half an ulp gives the termination criteria we want.
+  using std::cos;
+  using std::abs;
+  using std::sin;
+  T x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"};
+  T residual = cos(x) - x;
+  do {
+    x += residual/(sin(x)+1);
+    residual = cos(x) - x;
+  } while(abs(residual) > std::numeric_limits<T>::epsilon());
+  return x;
+}
+
+
+template <class T>
+template<int N>
+inline T constant_reciprocal_fibonacci<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((boost::integral_constant<int, N>)))
+{
+  // Wikipedia says Gosper has deviced a faster algorithm for this, but I read the linked paper and couldn't see it!
+  // In any case, k bits per iteration is fine, though it would be better to sum from smallest to largest.
+  // That said, the condition number is unity, so it should be fine.
+  T x0 = 1;
+  T x1 = 1;
+  T sum = 2;
+  T diff = 1;
+  while (diff > std::numeric_limits<T>::epsilon()) {
+    T tmp = x1 + x0;
+    diff = 1/tmp;
+    sum += diff;
+    x0 = x1;
+    x1 = tmp;
+  }
+  return sum;
+}
+
 #endif
 
 }

include/boost/math/constants/constants.hpp

@@ -325,6 +325,8 @@ namespace boost{ namespace math
   BOOST_DEFINE_MATH_CONSTANT(first_feigenbaum, 4.66920160910299067185320382046620161725818557747576863274,  "4.6692016091029906718532038204662016172581855774757686327456513430041343302113147371386897440239480138171");
   BOOST_DEFINE_MATH_CONSTANT(plastic, 1.324717957244746025960908854478097340734404056901733364534, "1.32471795724474602596090885447809734073440405690173336453401505030282785124554759405469934798178728032991");
   BOOST_DEFINE_MATH_CONSTANT(gauss, 0.834626841674073186281429732799046808993993013490347002449, "0.83462684167407318628142973279904680899399301349034700244982737010368199270952641186969116035127532412906785");
+  BOOST_DEFINE_MATH_CONSTANT(dottie, 0.739085133215160641655312087673873404013411758900757464965, "0.739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246");
+  BOOST_DEFINE_MATH_CONSTANT(reciprocal_fibonacci, 3.35988566624317755317201130291892717968890513, "3.35988566624317755317201130291892717968890513373196848649555381532513031899668338361541621645679008729704");
 #endif
 
 } // namespace constants

reporting/performance/constants_performance.cpp

@@ -10,6 +10,31 @@
 using namespace boost::math::constants;
 using boost::multiprecision::mpfr_float;
 
+void Dottie(benchmark::State& state)
+{
+    mpfr_float::default_precision(state.range(0));
+    for (auto _ : state)
+    {
+        benchmark::DoNotOptimize(dottie<mpfr_float>());
+    }
+    state.SetComplexityN(state.range(0));
+}
+
+BENCHMARK(Dottie)->RangeMultiplier(2)->Range(512, 1<<20)->Complexity()->Unit(benchmark::kMicrosecond);
+
+void ReciprocalFibonacci(benchmark::State& state)
+{
+    mpfr_float::default_precision(state.range(0));
+    for (auto _ : state)
+    {
+        benchmark::DoNotOptimize(reciprocal_fibonacci<mpfr_float>());
+    }
+    state.SetComplexityN(state.range(0));
+}
+
+BENCHMARK(ReciprocalFibonacci)->RangeMultiplier(2)->Range(512, 1<<20)->Complexity()->Unit(benchmark::kMicrosecond);
+
+
 void Pi(benchmark::State& state)
 {
     mpfr_float::default_precision(state.range(0));

test/test_constants.cpp

@@ -794,6 +794,22 @@ void test_gauss()
     CHECK_LE(Real(0.8346), G_computed);
 }
 
+template<typename Real>
+void test_dottie()
+{
+   using boost::math::constants::dottie;
+   using std::cos;
+   CHECK_ULP_CLOSE(dottie<Real>(), cos(dottie<Real>()), 1);
+}
+
+template<typename Real>
+void test_reciprocal_fibonacci()
+{
+   using boost::math::constants::reciprocal_fibonacci;
+   CHECK_LE(reciprocal_fibonacci<Real>(), Real(3.36));
+   CHECK_LE(Real(3.35), reciprocal_fibonacci<Real>());
+}
+
 #endif
 
 int main()
@@ -821,7 +837,10 @@ int main()
    test_gauss<double>();
    test_gauss<long double>();
    test_gauss<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<400>>>();
-
+   test_dottie<float>();
+   test_dottie<double>();
+   test_dottie<long double>();
+   test_dottie<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<400>>>();
 #endif
    return boost::math::test::report_errors();
 }