* line 251: calc_pi(result, ...digits2<...>::value + 10); Any particular reason for +10? I'd prefer calc_pi to do whatever adjustments are needed to get the right number of digits.
It looks to be an "implementation artifact", I'm investigating removing it.
detail/functions/pow.hpp:
* line 18 ff: pow_imp I don't like the fact that this implementation creates log2(p) temporaries on the stack and multiplies them all together on the way out of the recursion. I'd prefer to only use a constant number of temporaries.
* line 31: I don't see how the cases > 1 buy you anything. They don't change the number of calls to eval_multiply. (At least, it wouldn't if the algorithm were implemented correctly.) * line 47: ???: The work of this loop is repeated at every level of recursion. Please fix this function. It's a well-known algorithm. It shouldn't be hard to get it right.
Chris?
* line 76: pow_imp(denom, t, -p, mpl::false_()); This is broken for p == INT_MIN with twos-complement arithmetic.
Fixed locally.
* line 185: "The ldexp function..." I think you mean the "exp" function.
Nod, fixed locally.
* line 273: if(x.compare(ll) == 0) This is only legal if long long is in T::int_types. Better to use intmax_t which is guaranteed to be in T::int_types.
Nod. Fixed locally, there were a couple of other cases as well.
* line 287: const bool b_scale = (xx.compare(float_type(1e-4)) > 0); Unless I'm missing something, this will always be true because xx > 1 > 1e-4 as a result of the test on line 240.
Yep. Fixed locally. My fault for "improving" Chris's original algorithm.
* line 332: if(&arg == &result) This is totally unnecessary, since you already copy arg into a temporary using frexp on line 348.
Nod, fixed locally.
* line 351: if(t.compare(fp_type(2) / fp_type(3)) <= 0) Is the rounding error evaluating 2/3 important? I've convinced myself that it works, but it took a little thought. Maybe add a comment to the effect that the exact value of the boundary doesn't matter as long as its about 2/3?
Not sure, Chris?
* line 413: "The fabs function " s/fabs/log10/.
Fixed locally.
* line 464: eval_convert_to(&an, a); This makes the assumption that conversion to long long is always okay even if the value is out of the range of long long. This assumption makes me nervous.
Needs a try{}catch{} at the very least. Fixed locally.
* line 465: a.compare(an) compare(long long)
Yep.
* line 474: eval_subtract(T, long long) Also, this value of da is never used. It gets overwritten on line 480 or 485.
Nod. Fixed locally.
* line 493: (eval_get_sign(x) > 0) This expression should always be true at this point.
Nod. Fixed locally.
* line 614: *p_cosh = x; cosh is an even function, so this is wrong when x = -\inf
Nod. Fixed locally.
* line 625: T e_px, e_mx; These are only used in the generic implementation. It's better to restrict the scope of variables to the region that they're actually used in.
Nod. Fixed locally.
* line 643: eval_divide(*p_sinh, ui_type(2)); I notice that in other places you use ldexp for division by powers of 2. Which is better?
Good question, this is translated from Chris's original code which was specific to his decimal floating point code where division by 2 is (slightly) cheaper 'cos there's no simple ldexp. In the general case though, ldexp is usually a better bet as it's often a simple bit-fiddle. Changed locally.
detail/mp_number_base.hpp:
* line 423: (std::numeric_limits<T>::digits + 1) * 1000L This limits the number of bits to LONG_MAX/1000. I'd like there to be some way to know that as long as I use less than k bits, the library won't have any integer overflow problems. k can depend on INT_MAX, LONG_MAX, etc, but the implementation limits need to be clearly documented somewhere. I really hate the normal ad hoc situation where you just ignore the problem and hope that all the values are small enough to avoid overflow.
If there are really more digits than fit in a long then we'll likely have other issues - such as not enough memory to store them all in - however, I've added a static_assert and a note to this effect.
Other notes:
Backend Requirements:
* Do we always require that all number types be assignable from intmax_t, uintmax_t, and long double? Why must an integer be assignable from long double?
Good question. It feels right to me that wherever possible all types should have uniform requirements, but given that long double -> builtin int conversions are always possible, why not long double -> extended int? In other words it follows the idea that multiprecision types should as far as possible behave like builtin types. The conversion should probably be explicit in mp_number though (that's one on the TODO list).
* The code seems to assume that it's okay to pass int as the exp parameter of ldexp. These requirements don't mandate that this will work.
That's a bug in the requirements, will fix. Many thanks for the detailed comments, John.