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>
calc_log2(T & num,unsigned digits)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>
calc_e(T & result,unsigned digits)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 eval_multiply(denom, i);
121 eval_multiply(result, i);
122 eval_add(result, ui_type(1));
123 ++i;
124 }while(denom.compare(lim) <= 0);
125 eval_divide(result, denom);
126 }
127
128 template <class T>
calc_pi(T & result,unsigned digits)129 void calc_pi(T& result, unsigned digits)
130 {
131 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
132 typedef typename mpl::front<typename T::float_types>::type real_type;
133 //
134 // 1100 digits in string form:
135 //
136 const char* string_val = "3."
137 "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
138 "8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"
139 "4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273"
140 "7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094"
141 "3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912"
142 "9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132"
143 "0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235"
144 "4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859"
145 "5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303"
146 "5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
147 "3809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913152";
148 //
149 // Check if we can just construct from string:
150 //
151 if(digits < 3640) // 3640 binary digits ~ 1100 decimal digits
152 {
153 result = string_val;
154 return;
155 }
156
157 T a;
158 a = ui_type(1);
159 T b;
160 T A(a);
161 T B;
162 B = real_type(0.5f);
163 T D;
164 D = real_type(0.25f);
165
166 T lim;
167 lim = ui_type(1);
168 eval_ldexp(lim, lim, -(int)digits);
169
170 //
171 // This algorithm is from:
172 // Schonhage, A., Grotefeld, A. F. W., and Vetter, E. Fast Algorithms: A Multitape Turing
173 // Machine Implementation. BI Wissenschaftverlag, 1994.
174 // Also described in MPFR's algorithm guide: http://www.mpfr.org/algorithms.pdf.
175 //
176 // Let:
177 // a[0] = A[0] = 1
178 // B[0] = 1/2
179 // D[0] = 1/4
180 // Then:
181 // S[k+1] = (A[k]+B[k]) / 4
182 // b[k] = sqrt(B[k])
183 // a[k+1] = a[k]^2
184 // B[k+1] = 2(A[k+1]-S[k+1])
185 // D[k+1] = D[k] - 2^k(A[k+1]-B[k+1])
186 // Stop when |A[k]-B[k]| <= 2^(k-p)
187 // and PI = B[k]/D[k]
188
189 unsigned k = 1;
190
191 do
192 {
193 eval_add(result, A, B);
194 eval_ldexp(result, result, -2);
195 eval_sqrt(b, B);
196 eval_add(a, b);
197 eval_ldexp(a, a, -1);
198 eval_multiply(A, a, a);
199 eval_subtract(B, A, result);
200 eval_ldexp(B, B, 1);
201 eval_subtract(result, A, B);
202 bool neg = eval_get_sign(result) < 0;
203 if(neg)
204 result.negate();
205 if(result.compare(lim) <= 0)
206 break;
207 if(neg)
208 result.negate();
209 eval_ldexp(result, result, k - 1);
210 eval_subtract(D, result);
211 ++k;
212 eval_ldexp(lim, lim, 1);
213 }
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 {
do_nothingconstant_initializer222 static void do_nothing()
223 {
224 init.do_nothing();
225 }
226 private:
227 struct initializer
228 {
initializerconstant_initializer::initializer229 initializer()
230 {
231 F();
232 }
do_nothingconstant_initializer::initializer233 void do_nothing()const{}
234 };
235 static const initializer init;
236 };
237
238 template <class T, const T& (*F)(void)>
239 typename constant_initializer<T, F>::initializer const constant_initializer<T, F>::init;
240
241 template <class T>
get_constant_ln2()242 const T& get_constant_ln2()
243 {
244 static BOOST_MP_THREAD_LOCAL T result;
245 static BOOST_MP_THREAD_LOCAL bool b = false;
246 static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
247 if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
248 {
249 calc_log2(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
250 b = true;
251 digits = boost::multiprecision::detail::digits2<number<T> >::value();
252 }
253
254 constant_initializer<T, &get_constant_ln2<T> >::do_nothing();
255
256 return result;
257 }
258 #ifndef BOOST_MP_THREAD_LOCAL
259 #error 1
260 #endif
261
262 template <class T>
get_constant_e()263 const T& get_constant_e()
264 {
265 static BOOST_MP_THREAD_LOCAL T result;
266 static BOOST_MP_THREAD_LOCAL bool b = false;
267 static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
268 if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
269 {
270 calc_e(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
271 b = true;
272 digits = boost::multiprecision::detail::digits2<number<T> >::value();
273 }
274
275 constant_initializer<T, &get_constant_e<T> >::do_nothing();
276
277 return result;
278 }
279
280 template <class T>
get_constant_pi()281 const T& get_constant_pi()
282 {
283 static BOOST_MP_THREAD_LOCAL T result;
284 static BOOST_MP_THREAD_LOCAL bool b = false;
285 static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
286 if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
287 {
288 calc_pi(result, boost::multiprecision::detail::digits2<number<T, et_on> >::value());
289 b = true;
290 digits = boost::multiprecision::detail::digits2<number<T> >::value();
291 }
292
293 constant_initializer<T, &get_constant_pi<T> >::do_nothing();
294
295 return result;
296 }
297
298 template <class T>
get_constant_one_over_epsilon()299 const T& get_constant_one_over_epsilon()
300 {
301 static const bool is_init = false;
302 static T result;
303 if (is_init == false)
304 {
305 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
306 result = static_cast<ui_type>(1u);
307 eval_divide(result, std::numeric_limits<number<T> >::epsilon().backend());
308 }
309
310 constant_initializer<T, &get_constant_one_over_epsilon<T> >::do_nothing();
311
312 return result;
313 }
314
315