Subject: Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)
From: John Maddock (boost.regex_at_[hidden])
Date: 2012-06-25 08:19:53
> Here are my comments on the code so far.
> It's not very much, since I only got
> through a few files. I'll try to expand
> on this as I have time, but no promises.
Many thanks Steven.
> concepts/mp_number_architypes.hpp:
>
> * The file name is misspelled. It should
> be "archetypes."
Indeed, fixed in my local copy.
> * The #includes aren't quite right. The
> file does not use trunc, and does use
> fpclassify.
>
> * line 28: mpl::list requires #include <boost/mpl/list.hpp>
> * line 79: stringstream requires #include <sstream>
Fixed in my local copy.
> depricated:
>
> * Should be deprecated.
Will go away if the library is accepted anyway (just contains bits and bobs
not in the submission that I didn't want to delete just yet).
> detail/functions/constants.hpp:
>
> * line 12: typedef typename has_enough_bits<...> pred_type;
> If you're going to go through all this trouble to make
> sure that you have a type with enough bits, you shouldn't
> assume that unsigned has enough. IIRC, the standard
> only guarantees 16 bits for unsigned.
Nod. There's some new code that simplifies all that as well, simplified and
fixed in my local copy.
> * line 38: if(digits < 3640)
> Ugh. I'm assuming that digits means binary digits
> and 3640 = floor(log_2 10^1100)
Comments to that effect throughout.
> * line 51: I would appreciate a brief description
> of the formula. This code isn't very readable.
Nod. Added locally:
//
// We calculate log2 from using the formula:
//
// ln(2) = 3/4 SUM[n>=0] ((-1)^n * N!^2 / (2^n(2n+1)!))
//
// Numerator and denominator are calculated separately and then
// divided at the end, we also precalculate the terms up to n = 5
// since these fit in a 32-bit integer anyway.
//
// See Gourdon, X., and Sebah, P. The logarithmic constant: log 2, Jan.
2004.
// Also http://www.mpfr.org/algorithms.pdf.
//
Also spotted a missing optimisation in the arithmetic.
> * line 105: It looks like this is using the formula
> e \approx \sum_0^n{\frac{1}{i!}} = (\sum_0^n\frac{n!}{i!})/n!
> and induction using the rule:
> \sum_0^n\frac{n!}{i!} =
> (\sum_0^(n-1)\frac{(n-1)!}{i!}) * n + 1
>
> Please explain the math or point to somewhere
> that explains it.
Nod adding locally:
//
// Standard evaluation from the definition of e:
http://functions.wolfram.com/Constants/E/02/
//
> * line 147: ditto for pi
Nod. Adding locally:
//
// This algorithm is from:
// Schonhage, A., Grotefeld, A. F. W., and Vetter, E. Fast Algorithms: A
Multitape Turing
// Machine Implementation. BI Wissenschaftverlag, 1994.
// Also described in MPFR's algorithm guide:
http://www.mpfr.org/algorithms.pdf.
//
// Let:
// a[0] = A[0] = 1
// B[0] = 1/2
// D[0] = 1/4
// Then:
// S[k+1] = (A[k]+B[k]) / 4
// b[k] = sqrt(B[k])
// a[k+1] = a[k]^2
// B[k+1] = 2(A[k+1]-S[k+1])
// D[k+1] = D[k] - 2^k(A[k+1]-B[k+1])
// Stop when |A[k]-B[k]| <= 2^(k-p)
// and PI = B[k]/D[k]
Currently testing those changes before moving on to the other comments,
Thanks again, John.