boost/multiprecision/detail/functions/constants.hpp

   1 // Copyright 2011 John Maddock. Distributed under the Boost
   2 // Distributed under the Boost Software License, Version 1.0.
   3 //    (See accompanying file LICENSE_1_0.txt or copy at
   4 //          http://www.boost.org/LICENSE_1_0.txt)
   5 //
   6 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
   7 //
   8
   9 template <class T>
  10 void calc_log2(T& num, unsigned digits)
  11 {
  12    typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
  13    typedef typename mpl::front<typename T::signed_types>::type                         si_type;
  14
  15    //
  16    // String value with 1100 digits:
  17    //
  18    static const char* string_val = "0."
  19                                    "6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875"
  20                                    "4200148102057068573368552023575813055703267075163507596193072757082837143519030703862389167347112335"
  21                                    "0115364497955239120475172681574932065155524734139525882950453007095326366642654104239157814952043740"
  22                                    "4303855008019441706416715186447128399681717845469570262716310645461502572074024816377733896385506952"
  23                                    "6066834113727387372292895649354702576265209885969320196505855476470330679365443254763274495125040606"
  24                                    "9438147104689946506220167720424524529612687946546193165174681392672504103802546259656869144192871608"
  25                                    "2938031727143677826548775664850856740776484514644399404614226031930967354025744460703080960850474866"
  26                                    "3852313818167675143866747664789088143714198549423151997354880375165861275352916610007105355824987941"
  27                                    "4729509293113897155998205654392871700072180857610252368892132449713893203784393530887748259701715591"
  28                                    "0708823683627589842589185353024363421436706118923678919237231467232172053401649256872747782344535347"
  29                                    "6481149418642386776774406069562657379600867076257199184734022651462837904883062033061144630073719489";
  30    //
  31    // Check if we can just construct from string:
  32    //
  33    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
  34    {
  35       num = string_val;
  36       return;
  37    }
  38    //
  39    // We calculate log2 from using the formula:
  40    //
  41    // ln(2) = 3/4 SUM[n>=0] ((-1)^n * N!^2 / (2^n(2n+1)!))
  42    //
  43    // Numerator and denominator are calculated separately and then
  44    // divided at the end, we also precalculate the terms up to n = 5
  45    // since these fit in a 32-bit integer anyway.
  46    //
  47    // See Gourdon, X., and Sebah, P. The logarithmic constant: log 2, Jan. 2004.
  48    // Also http://www.mpfr.org/algorithms.pdf.
  49    //
  50    num = static_cast<ui_type>(1180509120uL);
  51    T denom, next_term, temp;
  52    denom        = static_cast<ui_type>(1277337600uL);
  53    next_term    = static_cast<ui_type>(120uL);
  54    si_type sign = -1;
  55
  56    ui_type limit = digits / 3 + 1;
  57
  58    for (ui_type n = 6; n < limit; ++n)
  59    {
  60       temp = static_cast<ui_type>(2);
  61       eval_multiply(temp, ui_type(2 * n));
  62       eval_multiply(temp, ui_type(2 * n + 1));
  63       eval_multiply(num, temp);
  64       eval_multiply(denom, temp);
  65       sign = -sign;
  66       eval_multiply(next_term, n);
  67       eval_multiply(temp, next_term, next_term);
  68       if (sign < 0)
  69          temp.negate();
  70       eval_add(num, temp);
  71    }
  72    eval_multiply(denom, ui_type(4));
  73    eval_multiply(num, ui_type(3));
  74    INSTRUMENT_BACKEND(denom);
  75    INSTRUMENT_BACKEND(num);
  76    eval_divide(num, denom);
  77    INSTRUMENT_BACKEND(num);
  78 }
  79
  80 template <class T>
  81 void calc_e(T& result, unsigned digits)
  82 {
  83    typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
  84    //
  85    // 1100 digits in string form:
  86    //
  87    const char* string_val = "2."
  88                             "7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274"
  89                             "2746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901"
  90                             "1573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069"
  91                             "5517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416"
  92                             "9283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312"
  93                             "7736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117"
  94                             "3012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509"
  95                             "9618188159304169035159888851934580727386673858942287922849989208680582574927961048419844436346324496"
  96                             "8487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016"
  97                             "7683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354"
  98                             "0212340784981933432106817012100562788023519303322474501585390473041995777709350366041699732972508869";
  99    //
 100    // Check if we can just construct from string:
 101    //
 102    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
 103    {
 104       result = string_val;
 105       return;
 106    }
 107
 108    T lim;
 109    lim = ui_type(1);
 110    eval_ldexp(lim, lim, digits);
 111
 112    //
 113    // Standard evaluation from the definition of e: http://functions.wolfram.com/Constants/E/02/
 114    //
 115    result = ui_type(2);
 116    T denom;
 117    denom     = ui_type(1);
 118    ui_type i = 2;
 119    do
 120    {
 121       eval_multiply(denom, i);
 122       eval_multiply(result, i);
 123       eval_add(result, ui_type(1));
 124       ++i;
 125    } while (denom.compare(lim) <= 0);
 126    eval_divide(result, denom);
 127 }
 128
 129 template <class T>
 130 void calc_pi(T& result, unsigned digits)
 131 {
 132    typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
 133    typedef typename mpl::front<typename T::float_types>::type    real_type;
 134    //
 135    // 1100 digits in string form:
 136    //
 137    const char* string_val = "3."
 138                             "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
 139                             "8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"
 140                             "4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273"
 141                             "7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094"
 142                             "3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912"
 143                             "9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132"
 144                             "0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235"
 145                             "4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859"
 146                             "5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303"
 147                             "5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
 148                             "3809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913152";
 149    //
 150    // Check if we can just construct from string:
 151    //
 152    if (digits < 3640) // 3640 binary digits ~ 1100 decimal digits
 153    {
 154       result = string_val;
 155       return;
 156    }
 157
 158    T a;
 159    a = ui_type(1);
 160    T b;
 161    T A(a);
 162    T B;
 163    B = real_type(0.5f);
 164    T D;
 165    D = real_type(0.25f);
 166
 167    T lim;
 168    lim = ui_type(1);
 169    eval_ldexp(lim, lim, -(int)digits);
 170
 171    //
 172    // This algorithm is from:
 173    // Schonhage, A., Grotefeld, A. F. W., and Vetter, E. Fast Algorithms: A Multitape Turing
 174    // Machine Implementation. BI Wissenschaftverlag, 1994.
 175    // Also described in MPFR's algorithm guide: http://www.mpfr.org/algorithms.pdf.
 176    //
 177    // Let:
 178    // a[0] = A[0] = 1
 179    // B[0] = 1/2
 180    // D[0] = 1/4
 181    // Then:
 182    // S[k+1] = (A[k]+B[k]) / 4
 183    // b[k] = sqrt(B[k])
 184    // a[k+1] = a[k]^2
 185    // B[k+1] = 2(A[k+1]-S[k+1])
 186    // D[k+1] = D[k] - 2^k(A[k+1]-B[k+1])
 187    // Stop when |A[k]-B[k]| <= 2^(k-p)
 188    // and PI = B[k]/D[k]
 189
 190    unsigned k = 1;
 191
 192    do
 193    {
 194       eval_add(result, A, B);
 195       eval_ldexp(result, result, -2);
 196       eval_sqrt(b, B);
 197       eval_add(a, b);
 198       eval_ldexp(a, a, -1);
 199       eval_multiply(A, a, a);
 200       eval_subtract(B, A, result);
 201       eval_ldexp(B, B, 1);
 202       eval_subtract(result, A, B);
 203       bool neg = eval_get_sign(result) < 0;
 204       if (neg)
 205          result.negate();
 206       if (result.compare(lim) <= 0)
 207          break;
 208       if (neg)
 209          result.negate();
 210       eval_ldexp(result, result, k - 1);
 211       eval_subtract(D, result);
 212       ++k;
 213       eval_ldexp(lim, lim, 1);
 214    } while (true);
 215
 216    eval_divide(result, B, D);
 217 }
 218
 219 template <class T, const T& (*F)(void)>
 220 struct constant_initializer
 221 {
 222    static void do_nothing()
 223    {
 224       init.do_nothing();
 225    }
 226
 227  private:
 228    struct initializer
 229    {
 230       initializer()
 231       {
 232          F();
 233       }
 234       void do_nothing() const {}
 235    };
 236    static const initializer init;
 237 };
 238
 239 template <class T, const T& (*F)(void)>
 240 typename constant_initializer<T, F>::initializer const constant_initializer<T, F>::init;
 241
 242 template <class T>
 243 const T& get_constant_ln2()
 244 {
 245    static BOOST_MP_THREAD_LOCAL T    result;
 246    static BOOST_MP_THREAD_LOCAL long digits = 0;
 247 #ifndef BOOST_MP_USING_THREAD_LOCAL
 248    static BOOST_MP_THREAD_LOCAL bool b = false;
 249    constant_initializer<T, &get_constant_ln2<T> >::do_nothing();
 250
 251    if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 252    {
 253       b = true;
 254 #else
 255    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 256    {
 257 #endif
 258       calc_log2(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 259       digits = boost::multiprecision::detail::digits2<number<T> >::value();
 260    }
 261
 262    return result;
 263 }
 264
 265 template <class T>
 266 const T& get_constant_e()
 267 {
 268    static BOOST_MP_THREAD_LOCAL T    result;
 269    static BOOST_MP_THREAD_LOCAL long digits = 0;
 270 #ifndef BOOST_MP_USING_THREAD_LOCAL
 271    static BOOST_MP_THREAD_LOCAL bool b = false;
 272    constant_initializer<T, &get_constant_e<T> >::do_nothing();
 273
 274    if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 275    {
 276       b = true;
 277 #else
 278    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 279    {
 280 #endif
 281       calc_e(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 282       digits = boost::multiprecision::detail::digits2<number<T> >::value();
 283    }
 284
 285    return result;
 286 }
 287
 288 template <class T>
 289 const T& get_constant_pi()
 290 {
 291    static BOOST_MP_THREAD_LOCAL T    result;
 292    static BOOST_MP_THREAD_LOCAL long digits = 0;
 293 #ifndef BOOST_MP_USING_THREAD_LOCAL
 294    static BOOST_MP_THREAD_LOCAL bool b = false;
 295    constant_initializer<T, &get_constant_pi<T> >::do_nothing();
 296
 297    if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 298    {
 299       b = true;
 300 #else
 301    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 302    {
 303 #endif
 304       calc_pi(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
 305       digits = boost::multiprecision::detail::digits2<number<T> >::value();
 306    }
 307
 308    return result;
 309 }
 310
 311 template <class T>
 312 const T& get_constant_one_over_epsilon()
 313 {
 314    static BOOST_MP_THREAD_LOCAL T    result;
 315    static BOOST_MP_THREAD_LOCAL long digits = 0;
 316 #ifndef BOOST_MP_USING_THREAD_LOCAL
 317    static BOOST_MP_THREAD_LOCAL bool b = false;
 318    constant_initializer<T, &get_constant_one_over_epsilon<T> >::do_nothing();
 319
 320    if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 321    {
 322       b = true;
 323 #else
 324    if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
 325    {
 326 #endif
 327       typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
 328       result = static_cast<ui_type>(1u);
 329       eval_divide(result, std::numeric_limits<number<T> >::epsilon().backend());
 330    }
 331
 332    return result;
 333 }