93 lines · plain
1// Copyright 2020, Madhur Chauhan2 3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_SPECIAL_FIBO_HPP9#define BOOST_MATH_SPECIAL_FIBO_HPP10 11#include <boost/math/constants/constants.hpp>12#include <boost/math/policies/error_handling.hpp>13#include <cmath>14#include <limits>15 16#ifdef _MSC_VER17#pragma once18#endif19 20namespace boost {21namespace math {22 23namespace detail {24 constexpr double fib_bits_phi = 0.69424191363061730173879026;25 constexpr double fib_bits_deno = 1.1609640474436811739351597;26} // namespace detail27 28template <typename T>29inline BOOST_MATH_CXX14_CONSTEXPR T unchecked_fibonacci(unsigned long long n) noexcept(std::is_fundamental<T>::value) {30 // This function is called by the rest and computes the actual nth fibonacci number31 // First few fibonacci numbers: 0 (0th), 1 (1st), 1 (2nd), 2 (3rd), ...32 if (n <= 2) return n == 0 ? 0 : 1;33 /* 34 * This is based on the following identities by Dijkstra:35 * F(2*n-1) = F(n-1)^2 + F(n)^236 * F(2*n) = (2*F(n-1) + F(n)) * F(n)37 * The implementation is iterative and is unrolled version of trivial recursive implementation.38 */39 unsigned long long mask = 1;40 for (int ct = 1; ct != std::numeric_limits<unsigned long long>::digits && (mask << 1) <= n; ++ct, mask <<= 1)41 ;42 T a{1}, b{1};43 for (mask >>= 1; mask; mask >>= 1) {44 T t1 = a * a;45 a = 2 * a * b - t1, b = b * b + t1;46 if (mask & n) 47 t1 = b, b = b + a, a = t1; // equivalent to: swap(a,b), b += a;48 }49 return a;50}51 52template <typename T, class Policy>53T inline BOOST_MATH_CXX14_CONSTEXPR fibonacci(unsigned long long n, const Policy &pol) {54 // check for overflow using approximation to binet's formula: F_n ~ phi^n / sqrt(5)55 if (n > 20 && n * detail::fib_bits_phi - detail::fib_bits_deno > std::numeric_limits<T>::digits)56 return policies::raise_overflow_error<T>("boost::math::fibonacci<%1%>(unsigned long long)", "Possible overflow detected.", pol);57 return unchecked_fibonacci<T>(n);58}59 60template <typename T>61T inline BOOST_MATH_CXX14_CONSTEXPR fibonacci(unsigned long long n) {62 return fibonacci<T>(n, policies::policy<>());63}64 65// generator for next fibonacci number (see examples/reciprocal_fibonacci_constant.hpp)66template <typename T>67class fibonacci_generator {68 public:69 // return next fibonacci number70 T operator()() noexcept(std::is_fundamental<T>::value) {71 T ret = a;72 a = b, b = b + ret; // could've simply: swap(a, b), b += a;73 return ret;74 }75 76 // after set(nth), subsequent calls to the generator returns consecutive77 // fibonacci numbers starting with the nth fibonacci number78 void set(unsigned long long nth) noexcept(std::is_fundamental<T>::value) {79 n = nth;80 a = unchecked_fibonacci<T>(n);81 b = unchecked_fibonacci<T>(n + 1);82 }83 84 private:85 unsigned long long n = 0;86 T a = 0, b = 1;87};88 89} // namespace math90} // namespace boost91 92#endif93