detail/functions/constants.hpp

0001 // Copyright 2011 John Maddock.
0002 // Distributed under the Boost Software License, Version 1.0.
0003 //    (See accompanying file LICENSE_1_0.txt or copy at
0004 //          http://www.boost.org/LICENSE_1_0.txt)
0005 //
0006 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
0007 //
0008
0009 template <class T>
0010 void calc_log2(T& num, unsigned digits)
0011 {
0012    using ui_type = typename boost::multiprecision::detail::canonical<std::uint32_t, T>::type;
0013    using si_type = typename std::tuple_element<0, typename T::signed_types>::type                        ;
0014
0015    //
0016    // String value with 1100 digits:
0017    //
0018    static const char* string_val = "0."
0019                                    "6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875"
0020                                    "4200148102057068573368552023575813055703267075163507596193072757082837143519030703862389167347112335"
0021                                    "0115364497955239120475172681574932065155524734139525882950453007095326366642654104239157814952043740"
0022                                    "4303855008019441706416715186447128399681717845469570262716310645461502572074024816377733896385506952"
0023                                    "6066834113727387372292895649354702576265209885969320196505855476470330679365443254763274495125040606"
0024                                    "9438147104689946506220167720424524529612687946546193165174681392672504103802546259656869144192871608"
0025                                    "2938031727143677826548775664850856740776484514644399404614226031930967354025744460703080960850474866"
0026                                    "3852313818167675143866747664789088143714198549423151997354880375165861275352916610007105355824987941"
0027                                    "4729509293113897155998205654392871700072180857610252368892132449713893203784393530887748259701715591"
0028                                    "0708823683627589842589185353024363421436706118923678919237231467232172053401649256872747782344535347"
0029                                    "6481149418642386776774406069562657379600867076257199184734022651462837904883062033061144630073719489";
0030    //
0031    // Check if we can just construct from string:
0032    //
0033    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
0034    {
0035       num = string_val;
0036       return;
0037    }
0038    //
0039    // We calculate log2 from using the formula:
0040    //
0041    // ln(2) = 3/4 SUM[n>=0] ((-1)^n * N!^2 / (2^n(2n+1)!))
0042    //
0043    // Numerator and denominator are calculated separately and then
0044    // divided at the end, we also precalculate the terms up to n = 5
0045    // since these fit in a 32-bit integer anyway.
0046    //
0047    // See Gourdon, X., and Sebah, P. The logarithmic constant: log 2, Jan. 2004.
0048    // Also http://www.mpfr.org/algorithms.pdf.
0049    //
0050    num = static_cast<ui_type>(1180509120uL);
0051    T denom, next_term, temp;
0052    denom        = static_cast<ui_type>(1277337600uL);
0053    next_term    = static_cast<ui_type>(120uL);
0054    si_type sign = -1;
0055
0056    ui_type limit = digits / 3 + 1;
0057
0058    for (ui_type n = 6; n < limit; ++n)
0059    {
0060       temp = static_cast<ui_type>(2);
0061       eval_multiply(temp, ui_type(2 * n));
0062       eval_multiply(temp, ui_type(2 * n + 1));
0063       eval_multiply(num, temp);
0064       eval_multiply(denom, temp);
0065       sign = -sign;
0066       eval_multiply(next_term, n);
0067       eval_multiply(temp, next_term, next_term);
0068       if (sign < 0)
0069          temp.negate();
0070       eval_add(num, temp);
0071    }
0072    eval_multiply(denom, ui_type(4));
0073    eval_multiply(num, ui_type(3));
0074    INSTRUMENT_BACKEND(denom);
0075    INSTRUMENT_BACKEND(num);
0076    eval_divide(num, denom);
0077    INSTRUMENT_BACKEND(num);
0078 }
0079
0080 template <class T>
0081 void calc_e(T& result, unsigned digits)
0082 {
0083    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0084    //
0085    // 1100 digits in string form:
0086    //
0087    const char* string_val = "2."
0088                             "7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274"
0089                             "2746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901"
0090                             "1573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069"
0091                             "5517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416"
0092                             "9283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312"
0093                             "7736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117"
0094                             "3012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509"
0095                             "9618188159304169035159888851934580727386673858942287922849989208680582574927961048419844436346324496"
0096                             "8487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016"
0097                             "7683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354"
0098                             "0212340784981933432106817012100562788023519303322474501585390473041995777709350366041699732972508869";
0099    //
0100    // Check if we can just construct from string:
0101    //
0102    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
0103    {
0104       result = string_val;
0105       return;
0106    }
0107
0108    T lim;
0109    lim = ui_type(1);
0110    eval_ldexp(lim, lim, digits);
0111
0112    //
0113    // Standard evaluation from the definition of e: http://functions.wolfram.com/Constants/E/02/
0114    //
0115    result = ui_type(2);
0116    T denom;
0117    denom     = ui_type(1);
0118    ui_type i = 2;
0119    do
0120    {
0121       eval_multiply(denom, i);
0122       eval_multiply(result, i);
0123       eval_add(result, ui_type(1));
0124       ++i;
0125    } while (denom.compare(lim) <= 0);
0126    eval_divide(result, denom);
0127 }
0128
0129 template <class T>
0130 void calc_pi(T& result, unsigned digits)
0131 {
0132    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0133    using real_type = typename std::tuple_element<0, typename T::float_types>::type   ;
0134    //
0135    // 1100 digits in string form:
0136    //
0137    const char* string_val = "3."
0138                             "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
0139                             "8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"
0140                             "4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273"
0141                             "7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094"
0142                             "3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912"
0143                             "9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132"
0144                             "0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235"
0145                             "4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859"
0146                             "5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303"
0147                             "5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
0148                             "3809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913152";
0149    //
0150    // Check if we can just construct from string:
0151    //
0152    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
0153    {
0154       result = string_val;
0155       return;
0156    }
0157
0158    T a;
0159    a = ui_type(1);
0160    T b;
0161    T A(a);
0162    T B;
0163    B = real_type(0.5f);
0164    T D;
0165    D = real_type(0.25f);
0166
0167    T lim;
0168    lim = ui_type(1);
0169    eval_ldexp(lim, lim, -static_cast<int>(digits));
0170
0171    //
0172    // This algorithm is from:
0173    // Schonhage, A., Grotefeld, A. F. W., and Vetter, E. Fast Algorithms: A Multitape Turing
0174    // Machine Implementation. BI Wissenschaftverlag, 1994.
0175    // Also described in MPFR's algorithm guide: http://www.mpfr.org/algorithms.pdf.
0176    //
0177    // Let:
0178    // a[0] = A[0] = 1
0179    // B[0] = 1/2
0180    // D[0] = 1/4
0181    // Then:
0182    // S[k+1] = (A[k]+B[k]) / 4
0183    // b[k] = sqrt(B[k])
0184    // a[k+1] = a[k]^2
0185    // B[k+1] = 2(A[k+1]-S[k+1])
0186    // D[k+1] = D[k] - 2^k(A[k+1]-B[k+1])
0187    // Stop when |A[k]-B[k]| <= 2^(k-p)
0188    // and PI = B[k]/D[k]
0189
0190    unsigned k = 1;
0191
0192    do
0193    {
0194       eval_add(result, A, B);
0195       eval_ldexp(result, result, -2);
0196       eval_sqrt(b, B);
0197       eval_add(a, b);
0198       eval_ldexp(a, a, -1);
0199       eval_multiply(A, a, a);
0200       eval_subtract(B, A, result);
0201       eval_ldexp(B, B, 1);
0202       eval_subtract(result, A, B);
0203       bool neg = eval_get_sign(result) < 0;
0204       if (neg)
0205          result.negate();
0206       if (result.compare(lim) <= 0)
0207          break;
0208       if (neg)
0209          result.negate();
0210       eval_ldexp(result, result, static_cast<int>(k - 1u));
0211       eval_subtract(D, result);
0212       ++k;
0213       eval_ldexp(lim, lim, 1);
0214    } while (true);
0215
0216    eval_divide(result, B, D);
0217 }
0218
0219 template <class T>
0220 const T& get_constant_ln2()
0221 {
0222    static BOOST_MP_THREAD_LOCAL T    result;
0223    static BOOST_MP_THREAD_LOCAL long digits = 0;
0224    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
0225    {
0226       boost::multiprecision::detail::maybe_promote_precision(&result);
0227       calc_log2(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0228       digits = boost::multiprecision::detail::digits2<number<T> >::value();
0229    }
0230
0231    return result;
0232 }
0233
0234 template <class T>
0235 const T& get_constant_e()
0236 {
0237    static BOOST_MP_THREAD_LOCAL T    result;
0238    static BOOST_MP_THREAD_LOCAL long digits = 0;
0239    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
0240    {
0241       boost::multiprecision::detail::maybe_promote_precision(&result);
0242       calc_e(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0243       digits = boost::multiprecision::detail::digits2<number<T> >::value();
0244    }
0245
0246    return result;
0247 }
0248
0249 template <class T>
0250 const T& get_constant_pi()
0251 {
0252    static BOOST_MP_THREAD_LOCAL T             result;
0253    static BOOST_MP_THREAD_LOCAL long digits = 0;
0254    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
0255    {
0256       boost::multiprecision::detail::maybe_promote_precision(&result);
0257       calc_pi(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
0258       digits = boost::multiprecision::detail::digits2<number<T> >::value();
0259    }
0260
0261    return result;
0262 }
0263 #ifdef BOOST_MSVC
0264 #pragma warning(push)
0265 #pragma warning(disable : 4127) // conditional expression is constant
0266 #endif
0267 template <class T>
0268 const T& get_constant_one_over_epsilon()
0269 {
0270    static BOOST_MP_THREAD_LOCAL T             result;
0271    static BOOST_MP_THREAD_LOCAL long digits = 0;
0272    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
0273    {
0274       using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0275       boost::multiprecision::detail::maybe_promote_precision(&result);
0276       result = static_cast<ui_type>(1u);
0277       BOOST_IF_CONSTEXPR(std::numeric_limits<number<T> >::is_specialized)
0278          eval_divide(result, std::numeric_limits<number<T> >::epsilon().backend());
0279       else
0280          eval_ldexp(result, result, boost::multiprecision::detail::digits2<number<T> >::value() - 1);
0281       digits = boost::multiprecision::detail::digits2<number<T> >::value();
0282    }
0283
0284    return result;
0285 }
0286 #ifdef BOOST_MSVC
0287 #pragma warning(pop)
0288 #endif