math/special_functions/fibonacci.hpp

0001 // Copyright 2020, Madhur Chauhan
0002
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_SPECIAL_FIBO_HPP
0009 #define BOOST_MATH_SPECIAL_FIBO_HPP
0010
0011 #include <boost/math/constants/constants.hpp>
0012 #include <boost/math/policies/error_handling.hpp>
0013 #include <cmath>
0014 #include <limits>
0015
0016 #ifdef _MSC_VER
0017 #pragma once
0018 #endif
0019
0020 namespace boost {
0021 namespace math {
0022
0023 namespace detail {
0024    constexpr double fib_bits_phi = 0.69424191363061730173879026;
0025    constexpr double fib_bits_deno = 1.1609640474436811739351597;
0026 } // namespace detail
0027
0028 template <typename T>
0029 inline BOOST_CXX14_CONSTEXPR T unchecked_fibonacci(unsigned long long n) noexcept(std::is_fundamental<T>::value) {
0030     // This function is called by the rest and computes the actual nth fibonacci number
0031     // First few fibonacci numbers: 0 (0th), 1 (1st), 1 (2nd), 2 (3rd), ...
0032     if (n <= 2) return n == 0 ? 0 : 1;
0033     /*
0034      * This is based on the following identities by Dijkstra:
0035      *   F(2*n-1) = F(n-1)^2 + F(n)^2
0036      *   F(2*n)   = (2*F(n-1) + F(n)) * F(n)
0037      * The implementation is iterative and is unrolled version of trivial recursive implementation.
0038      */
0039     unsigned long long mask = 1;
0040     for (int ct = 1; ct != std::numeric_limits<unsigned long long>::digits && (mask << 1) <= n; ++ct, mask <<= 1)
0041         ;
0042     T a{1}, b{1};
0043     for (mask >>= 1; mask; mask >>= 1) {
0044         T t1 = a * a;
0045         a = 2 * a * b - t1, b = b * b + t1;
0046         if (mask & n)
0047             t1 = b, b = b + a, a = t1; // equivalent to: swap(a,b), b += a;
0048     }
0049     return a;
0050 }
0051
0052 template <typename T, class Policy>
0053 T inline BOOST_CXX14_CONSTEXPR fibonacci(unsigned long long n, const Policy &pol) {
0054     // check for overflow using approximation to binet's formula: F_n ~ phi^n / sqrt(5)
0055     if (n > 20 && n * detail::fib_bits_phi - detail::fib_bits_deno > std::numeric_limits<T>::digits)
0056         return policies::raise_overflow_error<T>("boost::math::fibonacci<%1%>(unsigned long long)", "Possible overflow detected.", pol);
0057     return unchecked_fibonacci<T>(n);
0058 }
0059
0060 template <typename T>
0061 T inline BOOST_CXX14_CONSTEXPR fibonacci(unsigned long long n) {
0062     return fibonacci<T>(n, policies::policy<>());
0063 }
0064
0065 // generator for next fibonacci number (see examples/reciprocal_fibonacci_constant.hpp)
0066 template <typename T>
0067 class fibonacci_generator {
0068   public:
0069     // return next fibonacci number
0070     T operator()() noexcept(std::is_fundamental<T>::value) {
0071         T ret = a;
0072         a = b, b = b + ret; // could've simply: swap(a, b), b += a;
0073         return ret;
0074     }
0075
0076     // after set(nth), subsequent calls to the generator returns consecutive
0077     // fibonacci numbers starting with the nth fibonacci number
0078     void set(unsigned long long nth) noexcept(std::is_fundamental<T>::value) {
0079         n = nth;
0080         a = unchecked_fibonacci<T>(n);
0081         b = unchecked_fibonacci<T>(n + 1);
0082     }
0083
0084   private:
0085     unsigned long long n = 0;
0086     T a = 0, b = 1;
0087 };
0088
0089 } // namespace math
0090 } // namespace boost
0091
0092 #endif