Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)

Subject: Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)
From: John Maddock (boost.regex_at_[hidden])
Date: 2012-06-25 13:44:14

> * 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.