diff options
| author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-12-06 15:40:25 +0000 |
|---|---|---|
| committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-12-06 15:40:25 +0000 |
| commit | 2580a2b391cb47540b7d33dbd38ea7d301d2036c (patch) | |
| tree | 0d3b41cf73225c9cd96127a16734ad483e26b19e | |
| parent | 4ddd0dcfbf54fbed9e9c6117a7248090dd253ac2 (diff) | |
| download | mpfr-2580a2b391cb47540b7d33dbd38ea7d301d2036c.tar.gz | |
Merged the latest changes from the trunk except those related to
mpfr_zeta (it is buggy and not fixed yet in the trunk, but the bug
will not be triggered here, allowing the tests to pass).
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/faithful@10971 280ebfd0-de03-0410-8827-d642c229c3f4
46 files changed, 1145 insertions, 619 deletions
@@ -543,6 +543,18 @@ a. Using MinGW which is simpler and less demanding than Cygwin. Contrary to Cygwin, it also provides native Windows code. + If you also use MSYS, you should use "make" for MSYS instead of + the "make" utility from MinGW-W64 or from GCC, which causes the + following error: + + libtool: warning: libobj name 'extract.Tpo -c -o extract.lo extract.lo' + may not contain shell special characters. + rm: unknown option -- c + + References about this issue and solution: + https://sourceforge.net/p/msys2/tickets/223/ + https://sympa.inria.fr/sympa/arc/mpfr/2016-07/msg00009.html + 2 - If you just want to make a binary with gcc, there is nothing to do: GMP, MPFR and the program compile exactly as under Linux. (It is recommended to pass --build=xxx-yyy-mingw64 to the GMP configure command, @@ -611,17 +611,12 @@ Table of contents: (The one done in r10573 allowed us to find a bug even without assertion checking.) -- after the fix of mpfr_can_round in r9883, tzeta has a different trace - from r9882, which may be expected as mpfr_can_round has an influence - in tgeneric.c, but a closer look is needed: check only the first - difference (with various GMP_CHECK_RANDOMIZE values), since after it, - the random numbers do not match. - With the 32-bit ABI, tzeta r9955 is even 2.17 times as slow as r9954; - the corresponding change: - - test_generic (2, 70, 5); - + test_generic (MPFR_PREC_MIN, 70, 5); - Find the real cause, whether a wrong result could be obtained before - the fix r9883, and whether the slowdown can be avoided. +- tzeta has been much slower since r9848 (which increases the precision + of the input for the low output precisions), at least with the x86 + 32-bit ABI. This seems to come from the fact that the working precision + in the mpfr_zeta implementation depends on the precision of the input. + Once mpfr_zeta has improved, change the last argument of test_generic + in tzeta.c back to 5 (as it was before r10667). ############################################################################## diff --git a/acinclude.m4 b/acinclude.m4 index 1fefe840e..85fe826cd 100644 --- a/acinclude.m4 +++ b/acinclude.m4 @@ -213,7 +213,15 @@ dnl Check for POSIX Thread. Since the AX_PTHREAD macro is not standard dnl (it is provided by autoconf-archive), we need to detect whether it dnl is left unexpanded, otherwise the configure script won't fail and dnl "make distcheck" won't give any error, yielding buggy tarballs! -m4_pattern_forbid([AX_PTHREAD]) +dnl The \b is necessary to avoid an error with recent ax_pthread.m4 +dnl (such as with Debian's autoconf-archive 20160320-1), which contains +dnl AX_PTHREAD_ZOS_MISSING, etc. It is not documented, but see: +dnl https://lists.gnu.org/archive/html/autoconf/2015-03/msg00011.html +dnl +dnl Note: each time a change is done in m4_pattern_forbid, autogen.sh +dnl should be tested with and without ax_pthread.m4 availability (in +dnl the latter case, there should be an error). +m4_pattern_forbid([AX_PTHREAD\b]) AX_PTHREAD([]) dnl Check for ISO C11 Thread diff --git a/doc/README.dev b/doc/README.dev index 0476892fd..f81470cbd 100644 --- a/doc/README.dev +++ b/doc/README.dev @@ -218,8 +218,8 @@ To make a release (for the MPFR team): with and without -std=c11, with and without --enable-thread-safe, in various FPU precisions (double, double extended and single) if the platform supports that (e.g. under Linux/x86, with GCC and its -mpc64 - option to simulate the FreeBSD/NetBSD behavior, where by default, the - x87 FPU is configured to round on 53 bits), and in various locales + option to simulate the FreeBSD / NetBSD 6- behavior, where by default, + the x87 FPU is configured to round on 53 bits), and in various locales (LC_ALL=tr_TR in particular, if installed). On x86, test with -m96bit-long-double and -m128bit-long-double. Try also with gcc's -fno-common option. @@ -273,7 +273,11 @@ To make a release (for the MPFR team): diagnostic messages at runtime, not a failure; GCC 5 supports -fno-sanitize-recover like clang. - Test with GCC's AddressSanitizer (-fsanitize=address). + Test with GCC's AddressSanitizer (-fsanitize=address). One needs to + unset LD_PRELOAD to avoid failures. Alternatively, -static-libasan + could be used, but there are currently issues with it: + https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=836855 + https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=836864 Test with i586-mingw32msvc under Wine (see below). @@ -282,9 +286,10 @@ To make a release (for the MPFR team): Check various warnings, in particular for obsolescent features. With GCC: "-Wall -Wold-style-declaration -Wold-style-definition -Wmissing-parameter-type -Wmissing-prototypes -Wmissing-declarations - -Wmissing-field-initializers". They can easily be checked in - automatic tests by adding "-Werror -Wno-error=unused-function", - but this needs: + -Wmissing-field-initializers". The -Wint-in-bool-context option + could be added once available. These warnings can easily be checked + in automatic tests by adding "-Werror -Wno-error=unused-function + -Wno-error=deprecated-declarations", but this needs: * GCC 4.9+ * a patched autoconf: http://lists.gnu.org/archive/html/autoconf-patches/2014-01/msg00003.html diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 2d9081303..e43691f32 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -4253,7 +4253,7 @@ where \subsection{The Digamma Function} The Digamma function \texttt{mpfr\_digamma} is defined by: -\[ \psi(x) = \frac{d}{dx} \log\Gamma(x), \] +\[ \psi(x) = \frac{d}{{\mathrm d}x} \log\Gamma(x), \] and is computed from the asymptotic series \cite{Smith01} \[ \psi(x) \sim \log x - \frac{1}{2x} - \sum_{n=1}^{\infty} \frac{B_{2n}}{2n x^{2n}}. \] @@ -4629,7 +4629,7 @@ $|\epsilon''| \leq 2 \epsilon + \epsilon' \leq (26 + 2^{k+7}) 2^{k-p} Euler's constant is computed using the formula $\gamma = S(n) - R(n) - \log n$ where: \[ S(n) = \sum_{k=1}^{\infty} \frac{n^k (-1)^{k-1}}{k! k}, \quad - R(n) = \int_n^{\infty} \frac{\exp(-u)}{u} du \sim \frac{\exp(-n)}{n} + R(n) = \int_n^{\infty} \frac{\exp(-u)}{u} {\mathrm d}u \sim \frac{\exp(-n)}{n} \sum_{k=0}^{\infty} \frac{k!}{(-n)^k}. \] This identity is attributed to Sweeney by Brent \cite{Brent78}. (See also \cite{Smith01}.) @@ -4702,7 +4702,7 @@ as in \cite{Brent78}: % and $n^{n-1}$ within at most $n-2$ ulps. % The division by $n^{n-1}$ and the multiplication introduce one more ulp of % error, thus the total error on $R'(n)$ is at most $n+2$ ulps. -Let us introduce $I_k = \int_n^{\infty} \frac{e^{-u}}{u^k} du$. +Let us introduce $I_k = \int_n^{\infty} \frac{e^{-u}}{u^k} {\mathrm d}u$. We have $R(n) = I_1$ and the recurrence $I_k = \frac{e^{-n}}{n^k} - k I_{k+1}$, which gives \[ R(n) = \frac{e^{-n}}{n} \sum_{k=0}^{n-2} \frac{k!}{(-n)^k} @@ -4716,7 +4716,7 @@ $e^t < 1+2t$ for $0 < t < 1$, thus $(x!)^2 \le 2\pi x^{2x} e^{-2x} (x+\frac{1}{3})$.} for $n \ge 1$, we have: \[ |R(n) - R'(n)| = (n-1)! I_n - \le \frac{n!}{n} \int_n^{\infty} \frac{e^{-n}}{u^n} du + \le \frac{n!}{n} \int_n^{\infty} \frac{e^{-n}}{u^n} {\mathrm d}u \le \frac{n^{n-1}}{e^n} \sqrt{2 \pi (n+1)} \frac{e^{-n}}{(n-1) n^{n-1}} \] and since $\sqrt{2 \pi (n+1)}/(n-1) \le 1$ for $n \ge 9$: @@ -4763,16 +4763,16 @@ $\gamma = S(n) - R(n) - \log n$ with $n$ such that $e^{-2n} \le \ulp(1/2) = 2^{-{\rm prec}}$, i.e.\ $n \ge {\rm prec} \frac{\log 2}{2}$: \begin{quote} $s \leftarrow S'(n)$ \\ -$l \leftarrow \log(n)$ \\ +$\ell \leftarrow \log(n)$ \\ $r \leftarrow R'(n)$ \\ -return $(s - l) - r$ +return $(s - \ell) - r$ \end{quote} Since the final result is in $[\frac{1}{2}, 1]$, and $R'(n) \le \frac{e^{-n}}{n}$, then $S'(n)$ approximates $\log n$. If the target precision is $m$, and -we use $m + \lceil \log_2({\rm prec}) \rceil$ bits to evaluate $s$ and $l$, -then the error on $s-l$ will be at most $3 \ulp(1/2)$, -so the error on $(s - l) - r$ is at most $5 \ulp(1/2)$, +we use $m + \lceil \log_2({\rm prec}) \rceil$ bits to evaluate $s$ and $\ell$, +then the error on $s-\ell$ will be at most $3 \ulp(1/2)$, +so the error on $(s - \ell) - r$ is at most $5 \ulp(1/2)$, and adding the $3 e^{-2n}$ truncation error, we get a bound of $8 \ulp(1/2)$. [\textbf{FIXME: check with new method to compute S}] diff --git a/doc/sum.txt b/doc/sum.txt index 60fa5c52c..3dc0f048a 100644 --- a/doc/sum.txt +++ b/doc/sum.txt @@ -103,9 +103,6 @@ where sum will contain the result, x is an array of pointers to the inputs, n is the length of this array, and rnd is the rounding mode. The return value of type int will be the usual ternary value. -FIXME: It is currently assumed that no input pointers point to sum, -i.e., that the output sum is not shared with an input. - If n = 0, then the result is +0, whatever the rounding mode. This is equivalent to mpfr_set_ui and mpfr_set_si on the integer 0, and this choice is consistent with IEEE 754-2008's sum reduction operation of @@ -132,11 +129,11 @@ This is equivalent to apply the following steps: * if all the inputs have the same sign (i.e., all +0 or all −0), then the result has the same sign as the inputs; * otherwise, either because all inputs are zeros with at least - a +0 and a −0, or because some inputs are non-zero (but they + a +0 and a −0, or because some inputs are nonzero (but they globally cancel), the result is +0, except for the MPFR_RNDD rounding mode, where it is −0. -5. Otherwise the exact result is a non-zero real number, and the +5. Otherwise the exact result is a nonzero real number, and the conventional rounding function is applied. Concerning the exceptions, one generates a NaN exception in cases (1) @@ -192,7 +189,7 @@ following. * If n = 1, use mpfr_set. * If n = 2, use mpfr_add. -2. Iterate on the n input numbers to: +2. Iterate over the n input numbers to: (A) detect singular values (NaN, infinity, zero); (B) among the regular values, get the maximum exponent. @@ -211,7 +208,7 @@ following. Thus it is initialized to MPFR_EXP_MIN. * rn, which will contain the number n' of regular inputs, i.e., - those which are non-zero finite numbers. + those which are nonzero finite numbers. * sign_inf, which will contain the sign of an infinity result. It is initialized to 0, meaning no infinities yet. When the first @@ -253,7 +250,7 @@ following. sum_aux (sum, x, n, rnd, maxexp, rn); i.e., in addition to the parameters of mpfr_sum, we pass: * the maximum exponent; - * the number of regular inputs, i.e., the number of non-zero + * the number of regular inputs, i.e., the number of nonzero inputs. This number will be used instead of n to determine bounds on the sum (to avoid internal overflows) and error bounds. @@ -310,11 +307,12 @@ output precision, denoted sq, for two reasons: accumulator is a multiple of GMP_NUMB_BITS. Since logn ≥ 2, we have dq ≥ 4 as wanted. +Accumulator precision: wq = cq + maxexp − minexp = cq + sq + dq. + + ┌──────────────────── wq ────────────────────┐ Accumulator: [--------]-----------------------------------] cq └─ maxexp sq + dq minexp ─┘ -Accumulator precision: wq = cq + maxexp − minexp = cq + sq + dq. - This idea is not sufficient because if the Table Maker's Dilemma occurs, the needed precision for the truncated sum will grow. In particular, one can easily reach a huge precision with a few small-precision inputs. For @@ -356,10 +354,10 @@ let e and err denote the exponents of the computed result and of the error bound respectively. While e − err is less than some given bound denoted prec, shift the accumulator, update maxexp and minexp, and reiterate. For the caller, this bound must be large enough in order to -reach some wanted accuracy. However it cannot be too large since the +reach some wanted accuracy. However, it cannot be too large since the accumulator has a limited precision: we will need to make sure that if a reiteration is needed, then the cause is a partial cancellation, -so that the shift count is non-zero, otherwise minexp would not change +so that the shift count is nonzero, otherwise minexp would not change and one would get an infinite loop. Details and formal definitions are given later. @@ -387,13 +385,13 @@ Inputs of the sum_raw function (see the code for the types): * n: size of this array (number of inputs, regular or not). * minexp: exponent of the least significant bit of the first block. * maxexp: exponent of the first block (exponent of its MSB + 1). - * tp: pointer to a temporary area (pre-allocated). + * tp: pointer to a temporary area (preallocated). * ts: size of this temporary area. * logn: ceil(log2(rn)), where rn is the number of regular inputs. * prec: lower bound for e − err (as described above). * ep: pointer to mpfr_exp_t (see below), or a null pointer. - * minexpp: pointer to mpfr_exp_t (see below). - * maxexpp: pointer to mpfr_exp_t (see below). + * minexpp: pointer to mpfr_exp_t (see below), or a null pointer. + * maxexpp: pointer to mpfr_exp_t (see below), or a null pointer. We require as preconditions (explanations are given later): * prec ≥ 1 @@ -403,11 +401,14 @@ This function returns 0 if the accumulator is 0 (which implies that the exact sum for this sum_raw invocation is 0), otherwise the number of cancelled bits (>= 1), defined as the number of identical bits on the most significant part of the accumulator. In the latter case, it -also returns the following data in variables passed by reference: - * in ep: the exponent e of the computed result (unless ep is NULL); +also returns the following data in variables passed by reference, +unless ep, minexpp and maxexpp are null pointers (such data are useful +only after the first invocation of sum_raw): + * in ep: the exponent e of the computed result; * in minexpp: the last value of minexp; - * in maxexpp: the new value of maxexp (for the next iteration after - the first invocation of sum_raw in the main code). + * in maxexpp: the last value of maxexp2 (which would be the new value + of maxexp for the next iteration, i.e. the first iteration of the + second invocation of sum_raw in case of TMD resolution). Notes: * All the bits of the accumulator are meaningful; minexp is also @@ -438,12 +439,12 @@ Notes: Note: This definition can also be justified as follows. Negating a number in two's complement (in order to obtain its magnitude) can be done by complementing all the bits of the infinite series, and - the intervals would be the same if one considered non-canonical + the intervals would be the same if one considered noncanonical binary expansions, i.e. with an infinite trailing string of 1's: [−2^e,−2^(e−1)]. The difference between the open/closed bounds for [−2^e,−2^(e−1)[ and ]−2^e,−2^(e−1)] is due to the fact that a number represented by a finite number of bits (i.e., a canonical binary - expansion) is changed to a non-canonical binary expansion after + expansion) is changed to a noncanonical binary expansion after complement: the infinite trailing string of 0's is changed to an infinite trailing string of 1's, and the conversion to a canonical binary expansion could change the exponent. @@ -466,7 +467,7 @@ Then one has a loop over the inputs x[i]. Each input is processed with the following steps: 1. If the input is not regular, skip it. - Note: if there are many non-regular inputs, it may be more efficient + Note: if there are many singular inputs, it may be more efficient to have an array pointing to the regular inputs only, but such a case is assumed to be rare, and the number of iterations of this inner loop is also limited by the relatively small number of regular @@ -498,7 +499,7 @@ the following steps: added to the carry; carry propagation is done with mpn_add_1 if the size of the destination is larger than the size of the block. Note: There may be still be a carry out, but it is just ignored. This - occurs when a negative value in the accumulator becomes non-negative, + occurs when a negative value in the accumulator becomes nonnegative, and this fact is part of the usual two's complement arithmetic. If x[i] is negative, we do similar computations by using mpn_sub_n for the subtraction and mpn_sub_1 to propagate borrows. @@ -538,7 +539,7 @@ then we return 0, otherwise we reiterate at the beginning of sum_raw with: * minexp = cq + maxexp2 − wq * maxexp = maxexp2 -so that we still have in this subcase: wq = cq + maxexp - minexp. +so that we still have in this subcase: wq = cq + maxexp − minexp. We can now assume that the truncated sum is not 0. @@ -549,14 +550,14 @@ unlimited and the bits of exponent less than minexp are regarded as 0's: [--- accumulator ---]000000... - minexp ─┘└─ first non-represented bit = 0 + minexp ─┘└─ first nonrepresented bit = 0 So, we need to check that the computed value of cancel matches this mathematical point of view: * If the cancelled bits are 0's, then, since the truncated sum is not 0, there must be at least a bit 1 in the accumulator. * If the cancelled bits are 1's, then this sequence entirely fits - in the accumulator, as the first non-represented bit is a 0. + in the accumulator, as the first nonrepresented bit is a 0. The analysis below maps the truncated sum to the destination without considering rounding yet. @@ -588,15 +589,16 @@ inequality ensures that the value of the accumulator will not change with the updated minexp: shiftq is subtracted from minexp at the same time. The reiteration is done with maxexp = maxexp2, as said above. -Example. If there is a additional iteration with maxexp2 = minexp − 4 +Example. If there is an additional iteration with maxexp2 = minexp − 4 and a shift of 26 bits (due to cancellation), here is the accumulator before and after the shift: Before shift: [--------]-----------------------------------] - └─ maxexp minexp ─┘ - + cq ├─ maxexp sq + dq minexp ─┤ + 00000000000000000000000000001----------------] + <--- identical bits (0) ---> <------- 26 zeros -------> - After shift: [-------------------00000000000000000000000000 + After shift: 001-----------------00000000000000000000000000 This iteration: minexp ─┘ ├─ maxexp2 │ Next iteration: └─ maxexp minexp ─┘ @@ -684,7 +686,7 @@ any longer, but it can help to understand, with its examples): └─ maxexp i.e., an overflow in two's complement. But leaving at least - 2 + max(0, err - e) itentical bits in the most significant part, + 2 + max(0, err − e) identical bits in the most significant part, such as 0011111111111111111111111111110101000000 @@ -714,13 +716,13 @@ any longer, but it can help to understand, with its examples): value negative) just with additions of positive numbers. Said otherwise, leaving at least i identical bits allows one to - represent numbers in [-2^(e+i-1),2^(e+i-1)[. The current sum is in - [-2^e,2^e[, and taking into account the error terms, one wants to - be able to represent arbitrary numbers in ]-2^e-2^err,2^e+2^err[. - So, i must be chosen in such a way that 2^e + 2^err ≤ 2^(e+i-1), - i.e., 2^0 + 2^(err-e) ≤ 2^(i-1). The smallest power of two larger - than or equal to 2^0 + 2^(err-e) has exponent 1 + max(0, err - e). - Hence the choice i ≥ 2 + max(0, err - e). + represent numbers in [−2^(e+i−1),2^(e+i−1)[. The current sum is in + [−2^e,2^e[, and taking into account the error terms, one wants to + be able to represent arbitrary numbers in ]−2^e−2^err,2^e+2^err[. + So, i must be chosen in such a way that 2^e + 2^err ≤ 2^(e+i−1), + i.e., 2^0 + 2^(err−e) ≤ 2^(i−1). The smallest power of two larger + than or equal to 2^0 + 2^(err−e) has exponent 1 + max(0, err − e). + Hence the choice i ≥ 2 + max(0, err − e). Algorithm and Implementation: Back to the Generic Case ────────────────────────────────────────────────────── @@ -733,12 +735,25 @@ Accumulator: [--------]-----------------------------------] where cq = logn + 1, sq is the target precision, and dq ≥ logn + 2. -TODO: Update for code supporting reuse. - -Memory is allocated both for the accumulator and for the temporary area -needed by sum_raw. For performance reasons, the allocation is done in -the stack if the size is small enough (see MPFR_TMP_LIMBS_ALLOC macro). -No other memory allocation will be needed (except for auto variables). +A single chunk of memory is allocated both for the accumulator and for +the temporary area needed by sum_raw. We also chose to allocate memory +for the possible TMD resolution (as explained later) in the same chunk; +this second accumulator will be useless in most cases (it is necessary +only if the TMD occurs and some input is reused as the output), but in +the current implementation, it takes at most two limbs in practice, +so that this does not make a noticeable difference. For performance +reasons, the memory is allocated in the stack instead of the heap if +its size is small enough (see MPFR_TMP_LIMBS_ALLOC macro). No other +memory allocation will be needed (except for auto variables). + +Note: Having a small-size accumulator for sum_raw, either for the main +computation or for the TMD resolution, is not the best choice for the +worst-case complexity. For the time being, we focus on correctness and +make sure that the implementation is fast on almost all cases and not +too slow on corner cases. In the future, we may want to fix a minimal +size for the accumulator or allow it to grow dynamically, for instance +in a geometric progression after a few iterations (similarly to what is +done for Ziv loops in the TMD resolution for mathematical functions). The accumulator is zeroed and sum_raw is invoked to compute an accurate approximation of the sum. Among its parameters, maxexp was @@ -754,7 +769,7 @@ IEEE 754 rules (positive, except for MPFR_RNDD, where it is negative), and return with ternary value = 0. Now, the accumulator contains the significand of a good approximation -to the non-zero exact sum. The corresponding exponent is e and the sign +to the nonzero exact sum. The corresponding exponent is e and the sign is determined from one of the cancelled bits. The exponent of the ulp for the target precision is u = e − sq. We have: * err = maxexp + logn (where maxexp is the maxexp2 value returned by @@ -782,16 +797,17 @@ inputs x[i] whose exponent was very small compared to minexp. In such a case, the bit r and some of the least significant bits x may fall outside of the accumulator, in which case they are regarded as 0's (still due to truncation). In the following, we will make sure that -we do not try to read non-represented bits. +we do not try to read nonrepresented bits. -When maxexp (value returned by sum_raw) is different from MPFR_EXP_MIN, -i.e., when some bits of the inputs have still not been considered, we -will need to determine whether the TMD occurs. In this case, we will -compute d = u − err, which is ≥ 3, and can be very large if maxexp is -very small; nevertheless, d is representable in a mpfr_exp_t since: +When maxexp (corresponding to the maxexp2 value returned by sum_raw) +is different from MPFR_EXP_MIN, i.e., when some bits of the inputs +have still not been considered, we will need to determine whether the +TMD occurs. In this case, we will compute d = u − err, which is ≥ 3, +and can be very large if maxexp is very small; nevertheless, d is +representable in a mpfr_exp_t since: * If maxexp < minexp, then maxexp is the exponent of an input x[i], - so that maxpexp ≥ MPFR_EMIN_MIN; and since u ≤ MPFR_EMAX_MAX, we + so that maxexp ≥ MPFR_EMIN_MIN; and since u ≤ MPFR_EMAX_MAX, we have d ≤ MPFR_EMAX_MAX − MPFR_EMIN_MIN, which is representable in a mpfr_exp_t as per definition of MPFR_EMIN_MIN and MPFR_EMAX_MAX in MPFR (this choice was made to simplify internal computations). @@ -822,70 +838,116 @@ the TMD occurs when: identical. Several things need to be considered, in arbitrary order: - * the copy to the destination; - * a shift (for the alignment), if the shift count is non-zero + * the copy to the destination (if the destination is used by an input, + the TMD may need to be resolved first); + * a shift (for the alignment), if the shift count is nonzero (this is the most probable case); * a negation/complement if the value is negative (cancelled bits = 1), since the significand of MPFR numbers uses the conventional sign + absolute value representation; * rounding (the TMD needs to be resolved first if it occurs). -Some of these operations will be merged, i.e., done at the same time, -for performance reasons, and this depends on the operations provided -by the mpn layer of GMP. Ideally, all these operations should be -merged together, but this is not possible with the current version of -GMP (6.1.0). + +It is more efficient to merge some of these operations, i.e., do them at +the same time, and this possibility depends on the operations provided +by the mpn layer of GMP. Ideally, all these operations should be merged +together, but this is not possible with the current version of GMP +(6.1.1). For practical reasons, the shift should be done before the rounding, so that all the bits are represented for the rounding. The copy itself should be done together with the shift or the negation, because this is where most of the limbs are changed in general. We chose to do it -with the shift as it is assumed that the proportion of non-zero shift +with the shift as it is assumed that the proportion of nonzero shift counts is higher than the proportion of negations. -Moreover, when the value is negative, the difference between negation -and complement is similar to the difference between rounding directions +Moreover, for negative values, the difference between negation and +complement is similar to the difference between rounding directions (these operations are identical on the real numbers, i.e., in infinite precision), so that negation/complement and rounding can naturally be merged, as detailed later. -We start by doing the first two operations at the same time: the bits -of exponents in ⟦max(u,minexp),e⟦ are copied with a possible shift to -the most significant part of the destination, and the least significant -limbs (if any) are zeroed. +We will do the following: -Note: By definition of e, the most significant bit that is copied is -the complement of the value of the cancelled bits. A variable pos will -contain its value, i.e., pos = 1 for a positive number, pos = 0 for a -negative number. - -The values of 3 variables are also determined: - * inex: 0 if the final sum is known to be exact, otherwise 1. - * rbit: the rounding bit (0 or 1) of the truncated sum, corrected - in some halfway cases (see below). +1. Determine how the result will be rounded. If the TMD occurs, it is + resolved at this step. + +2. Copy the truncated accumulator (shifted) to the destination. For + simplicity, after this step, the trailing bits of the destination + (present when the precision is not a multiple of GMP_NUMB_BITS) + contain garbage. Since rounding needs a specific operation on the + least significant limb, these trailing bits (located in this limb) + will be zeroed in the next step. + +3. Take the complement if the result is negative, and at the same time, + do the rounding and zero the trailing bits. + +4. Set the exponent and handle a possible overflow/underflow. + +Details for each of these 4 steps are given below. + +Algorithm and Implementation: Rounding Information +────────────────────────────────────────────────── + +The values of 3 variables are determined: + * inex: 0 if the final sum is known to be exact (which can be the + case only if maxexp = MPFR_EXP_MIN), otherwise 1. + * rbit: the rounding bit (0 or 1) of the truncated sum, changed to 0 + for halfway cases that will round toward −infinity if the rounding + mode is to nearest (so that this bit gives the rounding direction), + as explained below. * tmd: three possible values: 0: the TMD does not occur. 1: the TMD occurs on a machine number. 2: the TMD occurs on a midpoint. -All this is done by considering two cases: +Note: The value of inex will be used only if the TMD does not occur +(i.e. tmd = 0). So, inex could be left indeterminate when tmd ≠ 0, +but this would not simplify the current code. + +This is done by considering two cases: - * u > minexp. One just has a shift (mpn_lshift) or a copy (MPN_COPY). - There are no least significant limbs to zero. In this case, the - rounding bit, which is represented, is read. Then there are two - subcases: + * u > minexp. The rounding bit, which is represented, is read. Then + there are two subcases: * maxexp = MPFR_EXP_MIN. The sum in the accumulator is exact. Thus inex will be the value of rounding bit || sticky bit. In round to nearest, rbit = 1 will mean that the value is - to be rounded toward +inf; therefore we need to detect and - handle halfway cases now (information on the trailing bits - will not be kept). More precisely we need to determine the - sticky bit if + to be rounded toward +infinity, even for halfway cases as + it is easier to handle these cases now. The variable rbit + is initially set to the value of the rounding bit. We need + to determine the sticky bit if * rbit = 0, or * rbit = 1 and rnd = MPFR_RNDN and the least significant bit - of the truncated sq-bit significand is 0. In this case, if - the sticky bit is 0, this halfway value will have to be - rounded toward −infinity; so, rbit is changed to 0. + of the truncated sq-bit significand (i.e., the bit before + the rounding bit) is 0; if the sticky bit is 0, this + halfway value will have to be rounded toward −infinity, + so that rbit is changed to 0. Note that for sq ≥ 2, the + parity of the rounded significand does not depend on the + representation (two's complement or sign + magnitude); + that is why, even though the significand is currently + represented in two's complement, we round to even. To + illustrate this point, let us give an example with a + negative value: + …1110.1100[100000] (two's complement) + …1110.1100 (rounded with rbit = 0, to even) + …0001.0100 (magnitude) + where the bits inside the brackets are those after the + truncated sq-bit significand. If we had converted the + accumulator first, we would have obtained: + …0001.0011[100000] (magnitude) + …0001.0100 (rounded with rbit = 1, to even) + i.e., the same result. + For sq = 1, the IEEE 754 rule for halfway cases is to choose + the value larger in magnitude, i.e., round away from zero. + Discussion about this rule and errata: + http://grouper.ieee.org/groups/754/email/msg03907.html + http://speleotrove.com/misc/IEEE754-errata.html + therefore, in this case, we want to keep rbit to 1 for + positive values, and set it to 0 for negative values; but + it happens that this corresponds to the rule chosen for + sq ≥ 2 (since the only bit of the truncated significand + is $1$ for positive values and $0$ for negative values), + so that there is no need to distinguish cases in the code. And tmd is set to 0 because one can round correctly, knowing the exact sum. @@ -895,24 +957,31 @@ All this is done by considering two cases: we need to take into account the fact that not all the bits are represented). - * u ≤ minexp. One just has a shift (bit shifted in are 0's, which is - what we want) or a copy, and if there are remaining low significant - limbs in the destination, they are zeroed with MPN_ZERO. We set - rbit to 0 because the rounding bit is not represented, thus its - value is 0. If maxexp = MPFR_EXP_MIN, then both inex and tmd are - set to 0; otherwise they are set to 1 (the bits following the ulp - bit are not represented, thus are all 0's, implying that the TMD - occurs on a machine number). + * u ≤ minexp. The rounding bit is not represented (its value is 0), + thus rbit is set to 0. If maxexp = MPFR_EXP_MIN, then both inex and + tmd are set to 0; otherwise they are set to 1 (the bits following + the ulp bit are not represented, thus are all 0's, implying that + the TMD occurs on a machine number). + +We also determine the sign of the result: a variable neg is set to the +value of the most significant bit of the accumulator, and a variable +sgn to the corresponding sign. In short: + ┌──────────┬─────┬─────┐ + │ number │ neg │ sgn │ + ├──────────┼─────┼─────┤ + │ positive │ 0 │ 1 │ + │ negative │ 1 │ −1 │ + └──────────┴─────┴─────┘ Now we seek to determine how the value will be rounded, more precisely, -what correction will be done to the significand that has been copied -just above. We currently have a significand, a trailing term t in the -accumulator such that 0 ≤ t < 1 ulp (non-negative thanks to the two's -complement representation), and an error on the trailing term bounded by +what correction will be done to the significand that will be copied. +We currently have a significand, a trailing term t in the accumulator +such that 0 ≤ t < 1 ulp (nonnegative thanks to the two's complement +representation), and an error on the trailing term bounded by t' ≤ 2^(u−3) = 2^(−3) ulp in absolute value, so that the error ε on the significand satisfies −t' ≤ ε < 1 ulp + t'. Thus one has 4 correction cases, denoted by an integer value corr between −1 and 2, which depends -on ε, the sign of the significand, rbit, and the rounding mode: +on ε, the sign of the result (neg / sgn), rbit, and the rounding mode: −1: equivalent to nextbelow; 0: no correction; +1: equivalent to nextabove. @@ -929,36 +998,45 @@ We distinguish two cases: absolute value of the ternary value. Therefore: * for MPFR_RNDD, corr ← 0; * for MPFR_RNDU, corr ← inex; - * for MPFR_RNDZ, corr ← inex && !pos; - * for MPFR_RNDA, corr ← inex && pos; + * for MPFR_RNDZ, corr ← inex && neg; + * for MPFR_RNDA, corr ← inex && !neg; * for MPFR_RNDN, corr ← rbit. If inex ≠ 0 and corr = 0, we set inex to −1. The variable inex now contains the ternary value. * tmd ≠ 0. The TMD occurs and will be solved by determining the sign (−1, 0 or +1) of a secondary term with a second sum_raw invocation - with a small-precision accumulator. - - TODO: Update for code supporting reuse. Still explain what was done - before. - - Since the uncorrected significand (which consists of the bits of - exponents in ⟦max(u,minexp),e⟦ as mentioned above) has already been - copied from the accumulator to the destination, we can reuse the - memory of the accumulator, but its precision is now set to cq + dq - rounded up to the next multiple of the word size (GMP_NUMB_BITS). - There may remain some bits that have not been copied, but they will - be taken into account as described below. + using a small-precision accumulator. + + Note: In the code written before the support of reused inputs as the + output, the accumulator had already been copied to the destination, + so that a part of the memory of this accumulator could be reused for + the small-size accumulator for the TMD resolution. This is no longer + possible, but currently not a problem since the accumulator for the + TMD resolution takes at most only 2 limbs in practice; however, in + the future, we might want the accumulators to grow dynamically, as + explained above. + + We set up a new accumulator of size cq + dq (= wq − sq) rounded up + to the next multiple of the word size (GMP_NUMB_BITS); let us call + this size zq (it will correspond to wq in sum_raw). From the old + accumulator, bits whose exponent is in ⟦minexp,u⟦ (when u > minexp) + will not be copied to the destination; these bits will be taken + into account as described below. Let us recall that the d−1 bits from exponent u−2 to u−d (= err) are identical. We distinguish two subcases: * err ≥ minexp. The last two over the d−1 identical bits and the - following bits, i.e., the bits from err+1 to minexp, are shifted - to the most significant part of the accumulator. - - After this shift, we have minexp = err + 2 − wq, so that: - maxexp − minexp = wq − logn − 2 ≤ wq − cq, + following bits, i.e., the bits from err+1 to minexp, are copied + (possibly with a shift) to the most significant part of the new + accumulator. + + The minexp value of this new accumulator is thus defined as: + minexp = err + 2 − zq, + so that: + maxexp − minexp = (err − logn) − (err + 2 − zq) + = zq − logn − 2 ≤ zq − cq, satisfying the condition for ts. * err < minexp. Here at least one of the identical bits is not @@ -967,8 +1045,8 @@ We distinguish two cases: from maxexp, with cq bits reserved to avoid overflows, just like in the main sum. - Once the accumulator is set up, sum_raw is called with prec = 1, - satisfying the first sum_raw precondition (prec ≥ 1). And we have: + Then sum_raw is called with prec = 1, satisfying the first sum_raw + precondition (prec ≥ 1). And we have: wq ≥ cq + dq ≥ logn + 3 = logn + prec + 2, @@ -982,20 +1060,20 @@ We distinguish two cases: negative if MSB = 1. Otherwise the secondary term is 0, and so is its sign; however, for the halfway cases (tmd = 2), we want to eliminate the ambiguity of their rounding due to the even-rounding - rule by choosing a non-zero value for the sign: −1 if the truncated + rule by choosing a nonzero value for the sign: −1 if the truncated significand is even, +1 if it is odd. Then, from the values of the variables rnd (rounding mode), tmd, rbit (rounding bit), sst (sign of the secondary term, corrected - for the halfway cases), and pos (sign of the sum), we determine: - * the correction case (between −1 and +2); - * the ternary value (negative, zero, or positive). + for the halfway cases), and sgn (sign of the sum), we determine: + * the ternary value (negative, zero, or positive); + * the correction case (between −1 and +2). The different cases are summarized in the table below. The two lines corresponding to halfway cases are not possible since sst has been - changed to an equivalent non-zero value as said above. The rounding + changed to an equivalent nonzero value as said above. The rounding modes MPFR_RNDZ and MPFR_RNDA are not in this table since they are - handled like MPFR_RNDD and MPFR_RNDU depending on the value of pos + handled like MPFR_RNDD and MPFR_RNDU depending on the value of sgn (see the MPFR_IS_LIKE_RNDD and MPFR_IS_LIKE_RNDU macros). rnd tmd rbit sst correction ternary @@ -1031,18 +1109,45 @@ We distinguish two cases: At this point, the variable inex contains the correct ternary value (before the overflow/underflow detection) and we know the correction -that needs to be applied to the significand. So, the next step is to -change the representation of the significand, i.e., do a negation or -complement if the sum is negative, and apply the correction at the -same time. +that needs to be applied to the significand. + +Algorithm and Implementation: Copy/Shift to the Destination +─────────────────────────────────────────────────────────── + +First, we can set the sign of the MPFR number from the value of sgn. + +The bits of the accumulator that need to be taken into account for the +destination are those of exponents in ⟦max(u,minexp),e⟦ (if u < minexp, +the nonrepresented bits are seen as 0's). We distinguish two cases: + + * u > minexp. We need to copy the bits of exponents in ⟦u,e⟦, i.e. + all the bits are represented in the accumulator. We just have a + left shift (mpn_lshift) or a copy (MPN_COPY). In the process, some + bits of exponent < u can have been copied to the trailing bits; they + are seen as garbage. Since rounding will need a specific operation + on the least significant limb, these trailing bits (located in this + limb) will be zeroed at the same time in the next step. + + * u ≤ minexp. We just have a left shift (bits that are shifted in are + 0's as specified by GMP, which is what we want) or a copy, and if + there are remaining low significant limbs in the destination, they + are zeroed with MPN_ZERO. + +Note: By definition of e, the most significant bit that is copied is +the first bit after the cancelled bits: 1 for a positive number, 0 for +a negative number. + +Algorithm and Implementation: Complement and Rounding +───────────────────────────────────────────────────── + +Note: The various cases below are not in the same order as in the code, +but this is better for the explanations. -First, we can set the sign of the MPFR number from the value of pos. For the moment, let us assume that sq ≥ 2. We distinguish two cases: - * pos = 1 (positive sum). We first set the sign of the MPFR number - to positive. Since the significand can contain garbage in the - trailing bits (present when the precision is not a multiple of - GMP_NUMB_BITS), we set these trailing bits to 0 as required by + * neg = 0 (positive sum). Since the significand can contain garbage + in the trailing bits (present when the precision is not a multiple + of GMP_NUMB_BITS), we set these trailing bits to 0 as required by the format of MPFR numbers. If corr > 0, we need to add corr to the significand (we can see that this works even if corr = 2 and the significand contains all 1's, which was not obvious). This is @@ -1056,17 +1161,16 @@ For the moment, let us assume that sq ≥ 2. We distinguish two cases: carry out, meaning a change of binade, the most significant bit of the significand is set to 1 without touching the other bits (this is important because if corr = 2 and the significand has - only one limb, the least significant non-trailing bit may be 1), + only one limb, the least significant nontrailing bit may be 1), and e is incremented. If corr < 0, then it is −1, so that we - subtract 1 from the significand with mpn_sub_1. In case of - borrow out, meaning a change of binade, the most significant bit - of the significand is set to 1 so that all the (non-trailing) - bits of the significand are 1's, and e is decremented. - - * pos = 0 (negative sum). We first set the sign of the MPFR number - to negative. Contrary to the positive case, where we could do - MPN number + limb (with mpn_add_1) or MPN number − limb (with - mpn_sub_1), GMP does not provide a function to do an operation + subtract 1 from the significand with mpn_sub_1. If the MSB of the + significand becomes 0, meaning a change of binade, then it is set + back to 1 so that all the (nontrailing) bits of the significand + are 1's, and e is decremented. + + * neg = 1 (negative sum). Contrary to the positive case, where we + could do MPN number + limb (with mpn_add_1) or MPN number − limb + (with mpn_sub_1), GMP does not provide a function to do an operation like: limb − MPN number. However, we will show that this can be emulated with mpn_neg, which does a negation, and mpn_com, which does a complement. This allows us to avoid separate mpn_com (or @@ -1091,7 +1195,7 @@ For the moment, let us assume that sq ≥ 2. We distinguish two cases: = com(x) + (1 − corr) where −1 ≤ 1 − corr ≤ 2. We consider two subcases, leading to a - non-negative case for the correction, and a negative case: + nonnegative case for the correction, and a negative case: * corr ≤ 1, i.e., 1 − corr ≥ 0. We first compute the least significant limb by setting the @@ -1105,7 +1209,7 @@ For the moment, let us assume that sq ≥ 2. We distinguish two cases: (if the significand has more than one limb). If there is still a carry, i.e., if the significand has exactly one limb or if there is no borrow out of the mpn_neg, then - we handle the change of binade just like in the pos = 1 + we handle the change of binade just like in the positive case for corr > 0. * If there is no carry, we apply mpn_com on the next limbs (if the significand has more than one limb). There cannot @@ -1116,8 +1220,8 @@ For the moment, let us assume that sq ≥ 2. We distinguish two cases: Here we want to compute com(x) − 1, but GMP does not provide an operation for that. The fact is that a sequence of low significant bits 1 is invariant, and we need to do the loop - ourselves in C instead of uses an optimized assembly version - from GMP. However this may not be a problem in practice, as + ourselves in C instead of using an optimized assembly version + from GMP. However, this may not be a problem in practice, as the difference is probably not noticeable. When a limb with a zero is reached (there is at least one since the most significant bit of the significand is a 0), we compute its @@ -1127,7 +1231,7 @@ For the moment, let us assume that sq ≥ 2. We distinguish two cases: the most significant limb can lead to a change of binade; more precisely, this happens on the significand 01111…111, whose complement is 10000…000 and com(x) − 1 is 01111…111. - The change of binade is handled like in the pos = 1 case for + The change of binade is handled like in the positive case for corr < 0. If sq = 1, the solution described above does not work when we need to @@ -1143,6 +1247,9 @@ In the negative case, as forcing the only limb to MPFR_LIMB_HIGHBIT corresponds to the computation of com(x), we just add (1 − corr) to the exponent e, following the formula given in the case sq > 1. +Algorithm and Implementation: Exponent Consideration +──────────────────────────────────────────────────── + Finally, we set the exponent of the MPFR number to e (not with MPFR_SET_EXP because this exponent may not be in the current exponent range) and call mpfr_check_range as usual to check whether e is in @@ -1240,7 +1347,7 @@ is known, and the accumulator lies just after the temporary area. But the gain would probably not be noticeable in practice. Note 2: It was considered to determine the number of cancelled limbs -instead of the number of cancelled bits in order to avoid a non-trivial +instead of the number of cancelled bits in order to avoid a nontrivial shift between iterations. Moreover, on platforms with a MMU, thus most platforms where GNU MPFR would be used, we could arrange to reserve a specific address space for the accumulator and choose the shift count @@ -1263,7 +1370,7 @@ have the minimal precision 2, and the output precision is p: Here long carry propagation will occur for each addition of the initial iteration, so that the complexity will be O(n·p) instead of O(n+p) if we -choose to delay carry propagation (however such a choice may slower the +choose to delay carry propagation (however, such a choice may slower the average case and take more memory, such as around 3p instead of 2p). When a new iteration is needed due to cancellation, a second accumulator diff --git a/examples/version.c b/examples/version.c index 3459027b1..00d651dcd 100644 --- a/examples/version.c +++ b/examples/version.c @@ -29,9 +29,6 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include <gmp.h> #include <mpfr.h> -#define SIGNED_STR(V) ((V) < 0 ? "signed" : "unsigned") -#define SIGNED(I) SIGNED_STR((I) - (I) - 1) - /* The following failure can occur if GMP has been rebuilt with * a different ABI, e.g. * 1. GMP built with ABI=mode32. @@ -49,37 +46,156 @@ static void failure_test (void) mpfr_clear (x); } +static void patches (void) +{ + const char *p = mpfr_get_patches (); + printf ("MPFR patches: %s\n", p[0] ? p : "[none]"); +} + +#define STRINGIZE(S) #S +#define MAKE_STR(S) STRINGIZE(S) + +#define SIGNED_STR(V) ((V) < 0 ? "signed" : "unsigned") +#define SIGNED(I) SIGNED_STR((I) - (I) - 1) + int main (void) { unsigned long c; - mp_limb_t t[4] = { -1, -1, -1, -1 }; + mp_limb_t t[4]; + int i; + + /* Casts are for C++ compilers. */ + for (i = 0; i < (int) (sizeof (t) / sizeof (mp_limb_t)); i++) + t[i] = (mp_limb_t) -1; + /**************** Information about the C implementation ****************/ + + /* This is useful, as this can affect the behavior of MPFR. */ + +#define COMP "Compiler: " +#ifdef __INTEL_COMPILER +# ifdef __VERSION__ +# define ICCV " [" __VERSION__ "]" +# else +# define ICCV "" +# endif + printf (COMP "ICC %d.%d.%d" ICCV "\n", __INTEL_COMPILER / 100, + __INTEL_COMPILER % 100, __INTEL_COMPILER_UPDATE); +#elif (defined(__GNUC__) || defined(__clang__)) && defined(__VERSION__) +# ifdef __clang__ +# define COMP2 COMP +# else +# define COMP2 COMP "GCC " +# endif + printf (COMP2 "%s\n", __VERSION__); +#endif + +#if defined(__STDC__) || defined(__STDC_VERSION__) + printf ("C/C++: __STDC__ = " +#if defined(__STDC__) + MAKE_STR(__STDC__) +#else + "undef" +#endif + ", __STDC_VERSION__ = " +#if defined(__STDC_VERSION__) + MAKE_STR(__STDC_VERSION__) +#else + "undef" +#endif #if defined(__cplusplus) - printf ("A C++ compiler is used.\n"); + ", C++" +#endif + "\n"); +#endif + +#if defined(__GNUC__) + printf ("GNU compatibility: __GNUC__ = " MAKE_STR(__GNUC__) + ", __GNUC_MINOR__ = " +#if defined(__GNUC_MINOR__) + MAKE_STR(__GNUC_MINOR__) +#else + "undef" #endif + "\n"); +#endif + +#if defined(__ICC) || defined(__INTEL_COMPILER) + printf ("Intel compiler: __ICC = " +#if defined(__ICC) + MAKE_STR(__ICC) +#else + "undef" +#endif + ", __INTEL_COMPILER = " +#if defined(__INTEL_COMPILER) + MAKE_STR(__INTEL_COMPILER) +#else + "undef" +#endif + "\n"); +#endif + +#if defined(_WIN32) || defined(_MSC_VER) + printf ("MS Windows: _WIN32 = " +#if defined(_WIN32) + MAKE_STR(_WIN32) +#else + "undef" +#endif + ", _MSC_VER = " +#if defined(_MSC_VER) + MAKE_STR(_MSC_VER) +#else + "undef" +#endif + "\n"); +#endif + +#if defined(__GLIBC__) + printf ("GNU C library: __GLIBC__ = " MAKE_STR(__GLIBC__) + ", __GLIBC_MINOR__ = " +#if defined(__GLIBC_MINOR__) + MAKE_STR(__GLIBC_MINOR__) +#else + "undef" +#endif + "\n"); +#endif + + printf ("\n"); + /************************************************************************/ + +#if defined(__MPIR_VERSION) + printf ("MPIR .... Library: %-12s Header: %d.%d.%d\n", + mpir_version, __MPIR_VERSION, __MPIR_VERSION_MINOR, + __MPIR_VERSION_PATCHLEVEL); +#else printf ("GMP ..... Library: %-12s Header: %d.%d.%d\n", gmp_version, __GNU_MP_VERSION, __GNU_MP_VERSION_MINOR, __GNU_MP_VERSION_PATCHLEVEL); +#endif printf ("MPFR .... Library: %-12s Header: %s (based on %d.%d.%d)\n", mpfr_get_version (), MPFR_VERSION_STRING, MPFR_VERSION_MAJOR, MPFR_VERSION_MINOR, MPFR_VERSION_PATCHLEVEL); -#if MPFR_VERSION_MAJOR >= 3 - printf ("MPFR features: TLS = %s, decimal = %s", - mpfr_buildopt_tls_p () ? "yes" : "no", - mpfr_buildopt_decimal_p () ? "yes" : "no"); -# if MPFR_VERSION_MAJOR > 3 || MPFR_VERSION_MINOR >= 1 + printf ("MPFR features: TLS = %s", mpfr_buildopt_tls_p () ? "yes" : "no"); +#if MPFR_VERSION_MAJOR >= 4 + printf (", float128 = %s", mpfr_buildopt_float128_p () ? "yes" : "no"); +#endif + printf (", decimal = %s", mpfr_buildopt_decimal_p () ? "yes" : "no"); +#if MPFR_VERSION_MAJOR > 3 || MPFR_VERSION_MINOR >= 1 printf (", GMP internals = %s\nMPFR tuning: %s", mpfr_buildopt_gmpinternals_p () ? "yes" : "no", mpfr_buildopt_tune_case ()); -# endif +#endif /* 3.1 */ printf ("\n"); -#endif - printf ("MPFR patches: %s\n\n", mpfr_get_patches ()); + patches (); + printf ("\n"); #ifdef __GMP_CC printf ("__GMP_CC = \"%s\"\n", __GMP_CC); #endif @@ -105,10 +221,8 @@ int main (void) printf ("\n"); printf ("sizeof(mpfr_prec_t) = %d (%s type)\n", (int) sizeof(mpfr_prec_t), SIGNED_STR((mpfr_prec_t) -1)); -#if MPFR_VERSION_MAJOR >= 3 printf ("sizeof(mpfr_exp_t) = %d (%s type)\n", (int) sizeof(mpfr_exp_t), SIGNED_STR((mpfr_exp_t) -1)); -#endif #ifdef _MPFR_PREC_FORMAT printf ("_MPFR_PREC_FORMAT = %d\n", (int) _MPFR_PREC_FORMAT); #endif @@ -122,7 +236,7 @@ int main (void) printf ("_MPFR_EXP_FORMAT = %d\n", (int) _MPFR_EXP_FORMAT); #endif printf ("sizeof(mpfr_t) = %d\n", (int) sizeof(mpfr_t)); - + printf ("sizeof(mpfr_ptr) = %d\n", (int) sizeof(mpfr_ptr)); failure_test (); return 0; diff --git a/src/atan.c b/src/atan.c index 335144481..7c0269f65 100644 --- a/src/atan.c +++ b/src/atan.c @@ -38,9 +38,9 @@ static void mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab) { mpz_t *S, *Q, *ptoj; - mp_bitcnt_t n, i, k, j, l; /* unsigned type, which is >= unsigned long */ + mp_bitcnt_t n, h, j; /* unsigned type, which is >= unsigned long */ mpfr_exp_t diff, expo; - int im, done; + int im, i, k, l, done; mpfr_prec_t mult; mpfr_prec_t accu[MPFR_PREC_BITS], log2_nb_terms[MPFR_PREC_BITS]; mpfr_prec_t precy = MPFR_PREC(y); @@ -151,18 +151,17 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab) } /* we need to combine S[0]/Q[0]...S[k-1]/Q[k-1] */ - l = 0; /* number of terms accumulated in S[k]/Q[k] */ + h = 0; /* number of terms accumulated in S[k]/Q[k] */ while (k > 1) { k --; /* combine S[k-1]/Q[k-1] and S[k]/Q[k] */ - j = log2_nb_terms[k-1]; mpz_mul (S[k], S[k], Q[k-1]); if (mpz_cmp_ui (p, 1) != 0) - mpz_mul (S[k], S[k], ptoj[j]); + mpz_mul (S[k], S[k], ptoj[log2_nb_terms[k-1]]); mpz_mul (S[k-1], S[k-1], Q[k]); - l += 1 << log2_nb_terms[k]; - mpz_mul_2exp (S[k-1], S[k-1], r * l); + h += (mp_bitcnt_t) 1 << log2_nb_terms[k]; + mpz_mul_2exp (S[k-1], S[k-1], r * h); mpz_add (S[k-1], S[k-1], S[k]); mpz_mul (Q[k-1], Q[k-1], Q[k]); } diff --git a/src/const_catalan.c b/src/const_catalan.c index 41b6bd97d..d7daa36e1 100644 --- a/src/const_catalan.c +++ b/src/const_catalan.c @@ -24,15 +24,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" /* Declare the cache */ -MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_catalan, mpfr_const_catalan_internal); - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -mpfr_cache_t * -__gmpfr_cache_const_catalan_f() -{ - return &__gmpfr_cache_const_catalan; -} -#endif +MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_catalan, mpfr_const_catalan_internal) /* Set User Interface */ #undef mpfr_const_catalan diff --git a/src/const_euler.c b/src/const_euler.c index 51b6a2c09..d49d4a68d 100644 --- a/src/const_euler.c +++ b/src/const_euler.c @@ -27,15 +27,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" /* Declare the cache */ -MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_euler, mpfr_const_euler_internal); - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -mpfr_cache_t * -__gmpfr_cache_const_euler_f() -{ - return &__gmpfr_cache_const_euler; -} -#endif +MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_euler, mpfr_const_euler_internal) /* Set User Interface */ #undef mpfr_const_euler diff --git a/src/const_log2.c b/src/const_log2.c index 91c9afe72..be11981fb 100644 --- a/src/const_log2.c +++ b/src/const_log2.c @@ -25,21 +25,11 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., /* Declare the cache */ #ifndef MPFR_USE_LOGGING -MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2, mpfr_const_log2_internal); +MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_log2, mpfr_const_log2_internal) #else -MPFR_DECL_INIT_CACHE(__gmpfr_normal_log2, mpfr_const_log2_internal); -MPFR_DECL_INIT_CACHE(__gmpfr_logging_log2, mpfr_const_log2_internal); -MPFR_THREAD_ATTR mpfr_cache_ptr __gmpfr_cache_const_log2 = __gmpfr_normal_log2; -#endif - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -# ifndef MPFR_USE_LOGGING -mpfr_cache_t * __gmpfr_cache_const_log2_f() { return &__gmpfr_cache_const_log2; } -# else -mpfr_cache_t * __gmpfr_normal_log2_f() { return &__gmpfr_normal_log2; } -mpfr_cache_t * __gmpfr_logging_log2_f() { return &__gmpfr_logging_log2; } -mpfr_cache_ptr * __gmpfr_cache_const_log2_f() { return &__gmpfr_cache_const_log2; } -# endif +MPFR_DECL_INIT_CACHE (__gmpfr_normal_log2, mpfr_const_log2_internal) +MPFR_DECL_INIT_CACHE (__gmpfr_logging_log2, mpfr_const_log2_internal) +MPFR_THREAD_VAR (mpfr_cache_ptr, __gmpfr_cache_const_log2, __gmpfr_normal_log2) #endif /* Set User interface */ diff --git a/src/const_pi.c b/src/const_pi.c index 5945cc9e0..e9682ed97 100644 --- a/src/const_pi.c +++ b/src/const_pi.c @@ -24,21 +24,11 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., /* Declare the cache */ #ifndef MPFR_USE_LOGGING -MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_pi, mpfr_const_pi_internal); +MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_pi, mpfr_const_pi_internal) #else -MPFR_DECL_INIT_CACHE(__gmpfr_normal_pi, mpfr_const_pi_internal); -MPFR_DECL_INIT_CACHE(__gmpfr_logging_pi, mpfr_const_pi_internal); -MPFR_THREAD_ATTR mpfr_cache_ptr __gmpfr_cache_const_pi = __gmpfr_normal_pi; -#endif - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -# ifndef MPFR_USE_LOGGING -mpfr_cache_t * __gmpfr_cache_const_pi_f() { return &__gmpfr_cache_const_pi; } -# else -mpfr_cache_t * __gmpfr_normal_pi_f() { return &__gmpfr_normal_pi; } -mpfr_cache_t * __gmpfr_logging_pi_f() { return &__gmpfr_logging_pi; } -mpfr_cache_ptr * __gmpfr_cache_const_pi_f() { return &__gmpfr_cache_const_pi; } -# endif +MPFR_DECL_INIT_CACHE (__gmpfr_normal_pi, mpfr_const_pi_internal) +MPFR_DECL_INIT_CACHE (__gmpfr_logging_pi, mpfr_const_pi_internal) +MPFR_THREAD_VAR (mpfr_cache_ptr, __gmpfr_cache_const_pi, __gmpfr_normal_pi) #endif /* Set User Interface */ @@ -31,123 +31,10 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #if !defined(MPFR_GENERIC_ABI) -#if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_INVERT_LIMB) -/* The following macro is copied from GMP-6.1.1, file gmp-impl.h, - macro udiv_qrnnd_preinv. - It computes q and r such that nh*2^GMP_NUMB_BITS + nl = q*d + r, - with 0 <= r < d, assuming di = __gmpn_invert_limb (d). */ -#define __udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ - do { \ - mp_limb_t _qh, _ql, _r, _mask; \ - umul_ppmm (_qh, _ql, (nh), (di)); \ - if (__builtin_constant_p (nl) && (nl) == 0) \ - { \ - _qh += (nh) + 1; \ - _r = - _qh * (d); \ - _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ - _qh += _mask; \ - _r += _mask & (d); \ - } \ - else \ - { \ - add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \ - _r = (nl) - _qh * (d); \ - _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ - _qh += _mask; \ - _r += _mask & (d); \ - if (MPFR_UNLIKELY (_r >= (d))) \ - { \ - _r -= (d); \ - _qh++; \ - } \ - } \ - (r) = _r; \ - (q) = _qh; \ - } while (0) - -/* specialized version for nl = 0, with di computed inside */ -#define __udiv_qrnd_preinv(q, r, nh, d) \ - do { \ - mp_limb_t _di; \ - \ - MPFR_ASSERTD ((d) != 0); \ - MPFR_ASSERTD ((nh) < (d)); \ - MPFR_ASSERTD ((d) & MPFR_LIMB_HIGHBIT); \ - \ - _di = __gmpn_invert_limb (d); \ - __udiv_qrnnd_preinv (q, r, nh, 0, d, _di); \ - } while (0) -#else -/* Same as __udiv_qrnnd_c from longlong.h, but using a single UWType/UWtype - division instead of two, and with n0=0. The analysis below assumes a limb - has 64 bits for simplicity. */ -#define __udiv_qrnd_preinv(q, r, n1, d) \ - do { \ - UWtype __d1, __d0, __q1, __q0, __r1, __r0, __i; \ - \ - MPFR_ASSERTD ((d) != 0); \ - MPFR_ASSERTD ((n1) < (d)); \ - MPFR_ASSERTD ((d) & MPFR_LIMB_HIGHBIT); \ - \ - __d1 = __ll_highpart (d); \ - /* 2^31 <= d1 < 2^32 */ \ - __d0 = __ll_lowpart (d); \ - /* 0 <= d0 < 2^32 */ \ - __i = ~(UWtype) 0 / __d1; \ - /* 2^32 < i < 2^33 with i < 2^64/d1 */ \ - \ - __q1 = (((n1) / __ll_B) * __i) / __ll_B; \ - /* Since n1 < d, high(n1) <= d1 = high(d), thus */ \ - /* q1 <= high(n1) * (2^64/d1) / 2^32 <= 2^32 */ \ - /* and also q1 <= n1/d1 thus r1 >= 0 below */ \ - __r1 = (n1) - __q1 * __d1; \ - while (__r1 >= __d1) \ - __q1 ++, __r1 -= __d1; \ - /* by construction, we have n1 = q1*d1+r1, and 0 <= r1 < d1 */ \ - /* thus q1 <= n1/d1 < 2^32+2 */ \ - /* q1*d0 <= (2^32+1)*(2^32-1) <= 2^64-1 thus it fits in a limb. */ \ - __r0 = __r1 * __ll_B; \ - __r1 = __r0 - __q1 * __d0; \ - /* At most two corrections are needed like in __udiv_qrnnd_c. */ \ - if (__r1 > __r0) /* borrow when subtracting q1*d0 */ \ - { \ - __q1--, __r1 += (d); \ - if (__r1 > (d)) /* no carry when adding d */ \ - __q1--, __r1 += (d); \ - } \ - /* We can have r1 < m here, but in this case the true remainder */ \ - /* is 2^64+r1, which is > m (same analysis as below for r0). */ \ - /* Now we have n1*2^32 = q1*d + r1, with 0 <= r1 < d. */ \ - MPFR_ASSERTD(__r1 < (d)); \ - \ - /* The same analysis as above applies, with n1 replaced by r1, */ \ - /* q1 by q0, r1 by r0. */ \ - __q0 = ((__r1 / __ll_B) * __i) / __ll_B; \ - __r0 = __r1 - __q0 * __d1; \ - while (__r0 >= __d1) \ - __q0 ++, __r0 -= __d1; \ - __r1 = __r0 * __ll_B; \ - __r0 = __r1 - __q0 * __d0; \ - /* We know the quotient floor(r1*2^64/d) is q0, q0-1 or q0-2.*/ \ - if (__r0 > __r1) /* borrow when subtracting q0*d0 */ \ - { \ - /* The quotient is either q0-1 or q0-2. */ \ - __q0--, __r0 += (d); \ - if (__r0 > (d)) /* no carry when adding d */ \ - __q0--, __r0 += (d); \ - } \ - /* Now we have n1*2^64 = (q1*2^32+q0)*d + r0, with 0 <= r0 < d. */ \ - MPFR_ASSERTD(__r0 < (d)); \ - \ - (q) = __q1 * __ll_B | __q0; \ - (r) = __r0; \ - } while (0) -#endif - /* special code for p=PREC(q) < GMP_NUMB_BITS, and PREC(u), PREC(v) <= GMP_NUMB_BITS */ static int -mpfr_divsp1 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) +mpfr_div_1 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { mpfr_prec_t p = MPFR_GET_PREC(q); mpfr_limb_ptr up = MPFR_MANT(u); @@ -267,7 +154,7 @@ mpfr_divsp1 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) /* special code for GMP_NUMB_BITS < PREC(q) < 2*GMP_NUMB_BITS and GMP_NUMB_BITS < PREC(u), PREC(v) <= 2*GMP_NUMB_BITS */ static int -mpfr_divsp2 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) +mpfr_div_2 (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) { mpfr_prec_t p = MPFR_GET_PREC(q); mpfr_limb_ptr up = MPFR_MANT(u); @@ -794,12 +681,12 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) #if !defined(MPFR_GENERIC_ABI) if (MPFR_GET_PREC(q) < GMP_NUMB_BITS && usize == 1 && vsize == 1) - return mpfr_divsp1 (q, u, v, rnd_mode); + return mpfr_div_1 (q, u, v, rnd_mode); #if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_INVERT_LIMB) if (GMP_NUMB_BITS < MPFR_GET_PREC(q) && MPFR_GET_PREC(q) < 2 * GMP_NUMB_BITS && usize == 2 && vsize == 2) - return mpfr_divsp2 (q, u, v, rnd_mode); + return mpfr_div_2 (q, u, v, rnd_mode); #endif #endif /* !defined(MPFR_GENERIC_ABI) */ diff --git a/src/exceptions.c b/src/exceptions.c index 2d94a7385..97026f81a 100644 --- a/src/exceptions.c +++ b/src/exceptions.c @@ -22,16 +22,9 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" -MPFR_THREAD_ATTR mpfr_flags_t __gmpfr_flags = 0; - -MPFR_THREAD_ATTR mpfr_exp_t __gmpfr_emin = MPFR_EMIN_DEFAULT; -MPFR_THREAD_ATTR mpfr_exp_t __gmpfr_emax = MPFR_EMAX_DEFAULT; - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -mpfr_flags_t * __gmpfr_flags_f() { return &__gmpfr_flags; } -mpfr_exp_t * __gmpfr_emin_f() { return &__gmpfr_emin; } -mpfr_exp_t * __gmpfr_emax_f() { return &__gmpfr_emax; } -#endif +MPFR_THREAD_VAR (mpfr_flags_t, __gmpfr_flags, 0) +MPFR_THREAD_VAR (mpfr_exp_t, __gmpfr_emin, MPFR_EMIN_DEFAULT) +MPFR_THREAD_VAR (mpfr_exp_t, __gmpfr_emax, MPFR_EMAX_DEFAULT) #undef mpfr_get_emin diff --git a/src/exp3.c b/src/exp3.c index 560099c9b..6b8998c12 100644 --- a/src/exp3.c +++ b/src/exp3.c @@ -42,7 +42,7 @@ static void mpfr_exp_rational (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *Q, mpfr_prec_t *mult) { - mp_bitcnt_t n, i, j; /* unsigned type, which is >= unsigned long */ + mp_bitcnt_t n, h, i, j; /* unsigned type, which is >= unsigned long */ mpz_t *S, *ptoj; mpfr_prec_t *log2_nb_terms; mpfr_exp_t diff, expo; @@ -114,14 +114,14 @@ mpfr_exp_rational (mpfr_ptr y, mpz_ptr p, long r, int m, /* accumulate all products in S[0] and Q[0]. Warning: contrary to above, here we do not have log2_nb_terms[k-1] = log2_nb_terms[k]+1. */ - l = 0; /* number of accumulated terms in the right part S[k]/Q[k] */ + h = 0; /* number of accumulated terms in the right part S[k]/Q[k] */ while (k > 0) { j = log2_nb_terms[k-1]; mpz_mul (S[k], S[k], ptoj[j]); mpz_mul (S[k-1], S[k-1], Q[k]); - l += 1 << log2_nb_terms[k]; - mpz_mul_2exp (S[k-1], S[k-1], r * l); + h += (mp_bitcnt_t) 1 << log2_nb_terms[k]; + mpz_mul_2exp (S[k-1], S[k-1], r * h); mpz_add (S[k-1], S[k-1], S[k]); mpz_mul (Q[k-1], Q[k-1], Q[k]); k--; diff --git a/src/get_d.c b/src/get_d.c index aa5751e8a..9108576be 100644 --- a/src/get_d.c +++ b/src/get_d.c @@ -151,8 +151,7 @@ mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode) return negative ? DBL_NEG_ZERO : 0.0; } - tmp[0] = *src; /* Hack copy mpfr_t */ - MPFR_SET_EXP (tmp, 0); + MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0); ret = mpfr_get_d (tmp, rnd_mode); if (MPFR_IS_PURE_FP(src)) diff --git a/src/get_ld.c b/src/get_ld.c index 777fb40ec..02f4dbe54 100644 --- a/src/get_ld.c +++ b/src/get_ld.c @@ -36,6 +36,12 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) #elif defined(HAVE_LDOUBLE_IEEE_EXT_LITTLE) +/* Note: The code will return a result with a 64-bit precision, even + if the rounding precision is only 53 bits like on FreeBSD and + NetBSD 6- (or with GCC's -mpc64 option to simulate this on other + platforms). This is consistent with how strtold behaves in these + cases, for instance. */ + /* special code for IEEE 754 little-endian extended format */ long double mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) @@ -50,7 +56,7 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG); inex = mpfr_set (tmp, x, rnd_mode); - mpfr_set_emin (-16382-63); + mpfr_set_emin (-16381-63); /* emin=-16444, see below */ mpfr_set_emax (16384); mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode); mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */ @@ -113,17 +119,12 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) long double mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) { - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) return (long double) mpfr_get_d (x, rnd_mode); else /* now x is a normal non-zero number */ { long double r; /* result */ - long double m; double s; /* part of result */ - mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */ - mpfr_t y, z; - int sign; #if defined(HAVE_LDOUBLE_MAYBE_DOUBLE_DOUBLE) if (MPFR_LDBL_MANT_DIG == 106) @@ -131,21 +132,36 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) /* Assume double-double format (as found with the PowerPC ABI). The generic code below isn't used because numbers with precision > 106 would not be supported. */ - sh = 0; /* force sh to 0 otherwise if say x = 2^1023 + 2^(-1074) - then after shifting mpfr_get_d (y, rnd_mode) will - underflow to 0 */ - mpfr_init2 (y, mpfr_get_prec (x)); - mpfr_init2 (z, IEEE_DBL_MANT_DIG); /* keep the precision small */ - mpfr_set (y, x, rnd_mode); /* exact */ s = mpfr_get_d (x, MPFR_RNDN); /* high part of x */ - mpfr_set_d (z, s, MPFR_RNDN); /* exact */ - mpfr_sub (y, x, z, MPFR_RNDN); /* exact */ - /* Add the second part of y (in the correct rounding mode). */ - r = (long double) s + (long double) mpfr_get_d (y, rnd_mode); + /* Let's first consider special cases separately. The test for + infinity is really needed to avoid a NaN result. The test + for NaN is mainly for optimization. The test for 0 is useful + to get the correct sign (assuming mpfr_get_d supports signed + zeros on the implementation). */ + if (s == 0 || DOUBLE_ISNAN (s) || DOUBLE_ISINF (s)) + r = (long double) s; + else + { + mpfr_t y, z; + + mpfr_init2 (y, mpfr_get_prec (x)); + mpfr_init2 (z, IEEE_DBL_MANT_DIG); /* keep the precision small */ + mpfr_set_d (z, s, MPFR_RNDN); /* exact */ + mpfr_sub (y, x, z, MPFR_RNDN); /* exact */ + /* Add the second part of y (in the correct rounding mode). */ + r = (long double) s + (long double) mpfr_get_d (y, rnd_mode); + mpfr_clear (z); + mpfr_clear (y); + } } else #endif { + long double m; + mpfr_exp_t sh; /* exponent shift -> x/2^sh is in the double range */ + mpfr_t y, z; + int sign; + /* First round x to the target long double precision, so that all subsequent operations are exact (this avoids double rounding problems). However if the format contains numbers that have more @@ -175,42 +191,42 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_sub (y, y, z, MPFR_RNDN); /* exact */ } while (!MPFR_IS_ZERO (y)); - } - mpfr_clear (z); - mpfr_clear (y); - - /* we now have to multiply back by 2^sh */ - MPFR_ASSERTD (r > 0); - if (sh != 0) - { - /* An overflow may occur (example: 0.5*2^1024) */ - while (r < 1.0) - { - r += r; - sh--; - } - - if (sh > 0) - m = 2.0; - else - { - m = 0.5; - sh = -sh; - } + mpfr_clear (z); + mpfr_clear (y); - for (;;) + /* we now have to multiply back by 2^sh */ + MPFR_ASSERTD (r > 0); + if (sh != 0) { - if (sh % 2) - r = r * m; - sh >>= 1; - if (sh == 0) - break; - m = m * m; + /* An overflow may occur (example: 0.5*2^1024) */ + while (r < 1.0) + { + r += r; + sh--; + } + + if (sh > 0) + m = 2.0; + else + { + m = 0.5; + sh = -sh; + } + + for (;;) + { + if (sh % 2) + r = r * m; + sh >>= 1; + if (sh == 0) + break; + m = m * m; + } } + if (sign < 0) + r = -r; } - if (sign < 0) - r = -r; return r; } } @@ -228,8 +244,7 @@ mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode) if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src))) return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode); - tmp[0] = *src; /* Hack copy mpfr_t */ - MPFR_SET_EXP (tmp, 0); + MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0); ret = mpfr_get_ld (tmp, rnd_mode); if (MPFR_IS_PURE_FP(src)) diff --git a/src/init2.c b/src/init2.c index b519e379c..ca114618b 100644 --- a/src/init2.c +++ b/src/init2.c @@ -26,7 +26,7 @@ MPFR_HOT_FUNCTION_ATTR void mpfr_init2 (mpfr_ptr x, mpfr_prec_t p) { mp_size_t xsize; - mpfr_limb_ptr tmp; + mpfr_size_limb_t *tmp; /* Check if we can represent the number of limbs * associated to the maximum of mpfr_prec_t*/ @@ -52,7 +52,7 @@ mpfr_init2 (mpfr_ptr x, mpfr_prec_t p) MPFR_ASSERTN (MPFR_PREC_COND (p)); xsize = MPFR_PREC2LIMBS (p); - tmp = (mpfr_limb_ptr) (*__gmp_allocate_func)(MPFR_MALLOC_SIZE(xsize)); + tmp = (mpfr_size_limb_t *) (*__gmp_allocate_func)(MPFR_MALLOC_SIZE(xsize)); MPFR_PREC(x) = p; /* Set prec */ MPFR_EXP (x) = MPFR_EXP_INVALID; /* make sure that the exp field has a diff --git a/src/log_ui.c b/src/log_ui.c index b544e0d64..4acaf20a6 100644 --- a/src/log_ui.c +++ b/src/log_ui.c @@ -120,8 +120,16 @@ mpfr_log_ui (mpfr_ptr x, unsigned long n, mpfr_rnd_t rnd_mode) /* here n >= 3 */ - /* argument reduction: compute k such that 2/3 <= n/2^k < 4/3, - i.e., 2^(k+1) <= 3n < 2^(k+2) */ + /* Argument reduction: compute k such that 2/3 <= n/2^k < 4/3, + i.e., 2^(k+1) <= 3n < 2^(k+2). + + FIXME: we could do better by considering n/(2^k*3^i*5^j), + which reduces the maximal distance to 1 from 1/3 to 1/8, + thus needing about 1.89 less terms in the Taylor expansion of + the reduced argument. Then log(2^k*3^i*5^j) can be computed + using a combination of log(16/15), log(25/24) and log(81/80), + see Section 6.5 of "A Fortran Multiple-Precision Arithmetic Package", + Richard P. Brent, ACM Transactions on Mathematical Software, 1978. */ mpz_init_set_ui (three_n, n); mpz_mul_ui (three_n, three_n, 3); diff --git a/src/mpfr-gmp.h b/src/mpfr-gmp.h index 3088338fa..4e8c0a120 100644 --- a/src/mpfr-gmp.h +++ b/src/mpfr-gmp.h @@ -246,11 +246,34 @@ __MPFR_DECLSPEC extern const struct bases mpfr_bases[257]; #undef ABS #undef MIN #undef MAX -#undef numberof #define ABS(x) ((x) >= 0 ? (x) : -(x)) #define MIN(l,o) ((l) < (o) ? (l) : (o)) #define MAX(h,i) ((h) > (i) ? (h) : (i)) -#define numberof(x) (sizeof (x) / sizeof ((x)[0])) + +/* Size of an array, as a constant expression. */ +#define numberof_const(x) (sizeof (x) / sizeof ((x)[0])) + +/* Size of an array, safe version but not a constant expression: + Since an array can silently be converted to a pointer, we check + that this macro is applied on an array, not a pointer. */ +#undef numberof +#if 0 +/* The following should work with GCC as documented in its manual, + but fails: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=38377#c10 + Thus disabled for now. */ +# define numberof(x) \ + ( __extension__ ({ \ + int is_array = (void *) &(x) == (void *) &(x)[0]; \ + MPFR_STAT_STATIC_ASSERT (__builtin_constant_p (is_array) ? \ + is_array : 1); \ + MPFR_ASSERTN (is_array); \ + numberof_const (x); \ + }) ) +#else +# define numberof(x) \ + (MPFR_ASSERTN ((void *) &(x) == (void *) &(x)[0]), \ + numberof_const (x)) +#endif /* Allocate func are defined in gmp-impl.h. Warning: the code below fetches the GMP memory allocation functions the first @@ -348,13 +371,20 @@ __MPFR_DECLSPEC void mpfr_tmp_free (struct tmp_marker *); It returns invxl = floor((B^2-1)/xl)-B, where B=2^BITS_PER_LIMB, assuming the most significant bit of xl is set. */ #ifndef invert_limb +#if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_INVERT_LIMB) +#define invert_limb(invxl,xl) \ + do { \ + invxl = __gmpn_invert_limb (xl); \ + } while (0) +#else #define invert_limb(invxl,xl) \ do { \ mp_limb_t dummy MPFR_MAYBE_UNUSED; \ MPFR_ASSERTD ((xl) != 0); \ udiv_qrnnd (invxl, dummy, ~(xl), MPFR_LIMB_MAX, xl); \ } while (0) -#endif +#endif /* defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_INVERT_LIMB) */ +#endif /* ifndef invert_limb */ /* udiv_qr_3by2 macro, adapted from GMP 5.0.2, file gmp-impl.h. Compute quotient the quotient and remainder for n / d. Requires d @@ -435,7 +465,118 @@ typedef struct {mp_limb_t inv32;} mpfr_pi1_t; # define mpn_copyd MPN_COPY #endif +#if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_INVERT_LIMB) +/* The following macro is copied from GMP-6.1.1, file gmp-impl.h, + macro udiv_qrnnd_preinv. + It computes q and r such that nh*2^GMP_NUMB_BITS + nl = q*d + r, + with 0 <= r < d, assuming di = __gmpn_invert_limb (d). */ +#define __udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ + do { \ + mp_limb_t _qh, _ql, _r, _mask; \ + umul_ppmm (_qh, _ql, (nh), (di)); \ + if (__builtin_constant_p (nl) && (nl) == 0) \ + { \ + _qh += (nh) + 1; \ + _r = - _qh * (d); \ + _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ + _qh += _mask; \ + _r += _mask & (d); \ + } \ + else \ + { \ + add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \ + _r = (nl) - _qh * (d); \ + _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ + _qh += _mask; \ + _r += _mask & (d); \ + if (MPFR_UNLIKELY (_r >= (d))) \ + { \ + _r -= (d); \ + _qh++; \ + } \ + } \ + (r) = _r; \ + (q) = _qh; \ + } while (0) +/* specialized version for nl = 0, with di computed inside */ +#define __udiv_qrnd_preinv(q, r, nh, d) \ + do { \ + mp_limb_t _di; \ + \ + MPFR_ASSERTD ((d) != 0); \ + MPFR_ASSERTD ((nh) < (d)); \ + MPFR_ASSERTD ((d) & MPFR_LIMB_HIGHBIT); \ + \ + _di = __gmpn_invert_limb (d); \ + __udiv_qrnnd_preinv (q, r, nh, 0, d, _di); \ + } while (0) +#else +/* Same as __udiv_qrnnd_c from longlong.h, but using a single UWType/UWtype + division instead of two, and with n0=0. The analysis below assumes a limb + has 64 bits for simplicity. */ +#define __udiv_qrnd_preinv(q, r, n1, d) \ + do { \ + UWtype __d1, __d0, __q1, __q0, __r1, __r0, __i; \ + \ + MPFR_ASSERTD ((d) != 0); \ + MPFR_ASSERTD ((n1) < (d)); \ + MPFR_ASSERTD ((d) & MPFR_LIMB_HIGHBIT); \ + \ + __d1 = __ll_highpart (d); \ + /* 2^31 <= d1 < 2^32 */ \ + __d0 = __ll_lowpart (d); \ + /* 0 <= d0 < 2^32 */ \ + __i = ~(UWtype) 0 / __d1; \ + /* 2^32 < i < 2^33 with i < 2^64/d1 */ \ + \ + __q1 = (((n1) / __ll_B) * __i) / __ll_B; \ + /* Since n1 < d, high(n1) <= d1 = high(d), thus */ \ + /* q1 <= high(n1) * (2^64/d1) / 2^32 <= 2^32 */ \ + /* and also q1 <= n1/d1 thus r1 >= 0 below */ \ + __r1 = (n1) - __q1 * __d1; \ + while (__r1 >= __d1) \ + __q1 ++, __r1 -= __d1; \ + /* by construction, we have n1 = q1*d1+r1, and 0 <= r1 < d1 */ \ + /* thus q1 <= n1/d1 < 2^32+2 */ \ + /* q1*d0 <= (2^32+1)*(2^32-1) <= 2^64-1 thus it fits in a limb. */ \ + __r0 = __r1 * __ll_B; \ + __r1 = __r0 - __q1 * __d0; \ + /* At most two corrections are needed like in __udiv_qrnnd_c. */ \ + if (__r1 > __r0) /* borrow when subtracting q1*d0 */ \ + { \ + __q1--, __r1 += (d); \ + if (__r1 > (d)) /* no carry when adding d */ \ + __q1--, __r1 += (d); \ + } \ + /* We can have r1 < m here, but in this case the true remainder */ \ + /* is 2^64+r1, which is > m (same analysis as below for r0). */ \ + /* Now we have n1*2^32 = q1*d + r1, with 0 <= r1 < d. */ \ + MPFR_ASSERTD(__r1 < (d)); \ + \ + /* The same analysis as above applies, with n1 replaced by r1, */ \ + /* q1 by q0, r1 by r0. */ \ + __q0 = ((__r1 / __ll_B) * __i) / __ll_B; \ + __r0 = __r1 - __q0 * __d1; \ + while (__r0 >= __d1) \ + __q0 ++, __r0 -= __d1; \ + __r1 = __r0 * __ll_B; \ + __r0 = __r1 - __q0 * __d0; \ + /* We know the quotient floor(r1*2^64/d) is q0, q0-1 or q0-2.*/ \ + if (__r0 > __r1) /* borrow when subtracting q0*d0 */ \ + { \ + /* The quotient is either q0-1 or q0-2. */ \ + __q0--, __r0 += (d); \ + if (__r0 > (d)) /* no carry when adding d */ \ + __q0--, __r0 += (d); \ + } \ + /* Now we have n1*2^64 = (q1*2^32+q0)*d + r0, with 0 <= r0 < d. */ \ + MPFR_ASSERTD(__r0 < (d)); \ + \ + (q) = __q1 * __ll_B | __q0; \ + (r) = __r0; \ + } while (0) +#endif /****************************************************** ************* GMP Basic Pointer Types **************** @@ -476,7 +617,11 @@ typedef const mp_limb_t *mpfr_limb_srcptr; of _Bool, signed int, unsigned int, or some other implementation-defined type. - So, unsigned int should be better. + So, unsigned int should be better. This will fail with implementations + having 16-bit int's, but such implementations are not required to + support bit-fields of size > 16 anyway; if ever an implementation with + 16-bit int's is found, the appropriate minimal changes could still be + done in the future. */ #ifndef _MPFR_IEEE_FLOATS diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 7e9fd60bd..64678e5f0 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -248,23 +248,24 @@ extern MPFR_CACHE_ATTR mpfr_cache_ptr __gmpfr_cache_const_log2; #endif #ifdef MPFR_WIN_THREAD_SAFE_DLL -__MPFR_DECLSPEC mpfr_flags_t * __gmpfr_flags_f(); -__MPFR_DECLSPEC mpfr_exp_t * __gmpfr_emin_f(); -__MPFR_DECLSPEC mpfr_exp_t * __gmpfr_emax_f(); -__MPFR_DECLSPEC mpfr_prec_t * __gmpfr_default_fp_bit_precision_f(); -__MPFR_DECLSPEC mpfr_rnd_t * __gmpfr_default_rounding_mode_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_euler_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_catalan_f(); +# define MPFR_MAKE_VARFCT(T,N) T * N ## _f (void) { return &N; } +__MPFR_DECLSPEC mpfr_flags_t * __gmpfr_flags_f (void); +__MPFR_DECLSPEC mpfr_exp_t * __gmpfr_emin_f (void); +__MPFR_DECLSPEC mpfr_exp_t * __gmpfr_emax_f (void); +__MPFR_DECLSPEC mpfr_prec_t * __gmpfr_default_fp_bit_precision_f (void); +__MPFR_DECLSPEC mpfr_rnd_t * __gmpfr_default_rounding_mode_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_euler_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_catalan_f (void); # ifndef MPFR_USE_LOGGING -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_pi_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_log2_f(); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_pi_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_log2_f (void); # else -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_normal_pi_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_normal_log2_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_logging_pi_f(); -__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_logging_log2_f(); -__MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_pi_f(); -__MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_log2_f(); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_normal_pi_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_normal_log2_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_logging_pi_f (void); +__MPFR_DECLSPEC mpfr_cache_t * __gmpfr_logging_log2_f (void); +__MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_pi_f (void); +__MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_log2_f (void); # endif # ifndef __MPFR_WITHIN_MPFR # define __gmpfr_flags (*__gmpfr_flags_f()) @@ -285,8 +286,14 @@ __MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_log2_f(); # define __gmpfr_cache_const_log2 (*__gmpfr_cache_const_log2_f()) # endif # endif +#else +# define MPFR_MAKE_VARFCT(T,N) #endif +# define MPFR_THREAD_VAR(T,N,V) \ + MPFR_THREAD_ATTR T N = (V); \ + MPFR_MAKE_VARFCT (T,N) + #define BASE_MAX 62 __MPFR_DECLSPEC extern const __mpfr_struct __gmpfr_l2b[BASE_MAX-1][2]; @@ -991,7 +998,7 @@ typedef uintmax_t mpfr_ueexp_t; following two macros, unless the flag comes from another function returning the ternary inexact value */ #define MPFR_RET(I) return \ - (I) ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0 + (I) != 0 ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0 #define MPFR_RET_NAN return (__gmpfr_flags |= MPFR_FLAGS_NAN), 0 #define SIGN(I) ((I) < 0 ? -1 : (I) > 0) @@ -1072,12 +1079,28 @@ typedef uintmax_t mpfr_ueexp_t; #define MPFR_BYTES_PER_MP_LIMB (GMP_NUMB_BITS/CHAR_BIT) -/* Heap memory handling */ +/* Heap memory handling + -------------------- + Memory allocated for a significand (mantissa) has the following + format: + * A mp_size_t in a mpfr_size_limb_t union (see below). + * An array of mp_limb_t (not all of them are necessarily used, + as the precision can change without a reallocation). + The goal of the mpfr_size_limb_t union is to make sure that + size and alignment requirements are satisfied if mp_size_t and + mp_limb_t have different sizes and/or alignment requirements. + Moreover, pointer conversions are not fully specified by the + C standard, and the use of a union (and the double casts below) + might help even if mp_size_t and mp_limb_t have the same size + and the same alignment requirements. Still, there is currently + no guarantee that this code is portable. Note that union members + are not used at all. +*/ typedef union { mp_size_t s; mp_limb_t l; } mpfr_size_limb_t; #define MPFR_GET_ALLOC_SIZE(x) \ - (((mp_size_t *) MPFR_MANT(x))[-1] + 0) + (((mp_size_t *) (mpfr_size_limb_t *) MPFR_MANT(x))[-1] + 0) #define MPFR_SET_ALLOC_SIZE(x, n) \ - (((mp_size_t *) MPFR_MANT(x))[-1] = (n)) + (((mp_size_t *) (mpfr_size_limb_t *) MPFR_MANT(x))[-1] = (n)) #define MPFR_MALLOC_SIZE(s) \ (sizeof(mpfr_size_limb_t) + MPFR_BYTES_PER_MP_LIMB * (size_t) (s)) #define MPFR_SET_MANT_PTR(x,p) \ @@ -1161,7 +1184,8 @@ typedef union { mp_size_t s; mp_limb_t l; } mpfr_size_limb_t; MPFR_CACHE_ATTR mpfr_cache_t _cache = {{ \ {{ 0, MPFR_SIGN_POS, 0, (mp_limb_t *) 0 }}, 0, _func \ MPFR_DEFERRED_INIT_SLAVE_VALUE(_func) \ - }} + }}; \ + MPFR_MAKE_VARFCT (mpfr_cache_t,_cache) /****************************************************** *************** Threshold parameters *************** @@ -211,7 +211,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) /* special code for prec(a) < GMP_NUMB_BITS and prec(b), prec(c) <= GMP_NUMB_BITS */ static int -mpfr_mulsp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, +mpfr_mul_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, mpfr_prec_t p) { mp_limb_t h; @@ -308,7 +308,7 @@ mpfr_mulsp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, /* special code for GMP_NUMB_BITS < prec(a) < 2*GMP_NUMB_BITS and GMP_NUMB_BITS < prec(b), prec(c) <= 2*GMP_NUMB_BITS */ static int -mpfr_mulsp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, +mpfr_mul_2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, mpfr_prec_t p) { mp_limb_t h, l, u, v; @@ -486,12 +486,12 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) cq = MPFR_GET_PREC (c); if (MPFR_GET_PREC(a) < GMP_NUMB_BITS && bq <= GMP_NUMB_BITS && cq <= GMP_NUMB_BITS) - return mpfr_mulsp1 (a, b, c, rnd_mode, MPFR_GET_PREC(a)); + return mpfr_mul_1 (a, b, c, rnd_mode, MPFR_GET_PREC(a)); if (GMP_NUMB_BITS < MPFR_GET_PREC(a) && MPFR_GET_PREC(a) < 2 * GMP_NUMB_BITS && GMP_NUMB_BITS < bq && bq <= 2 * GMP_NUMB_BITS && GMP_NUMB_BITS < cq && cq <= 2 * GMP_NUMB_BITS) - return mpfr_mulsp2 (a, b, c, rnd_mode, MPFR_GET_PREC(a)); + return mpfr_mul_2 (a, b, c, rnd_mode, MPFR_GET_PREC(a)); sign = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); diff --git a/src/mulders.c b/src/mulders.c index 4cc94cb31..d1be67541 100644 --- a/src/mulders.c +++ b/src/mulders.c @@ -41,7 +41,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., static short mulhigh_ktab[MPFR_MULHIGH_TAB_SIZE]; #else static short mulhigh_ktab[] = {MPFR_MULHIGH_TAB}; -#define MPFR_MULHIGH_TAB_SIZE ((mp_size_t) numberof (mulhigh_ktab)) +#define MPFR_MULHIGH_TAB_SIZE (numberof_const (mulhigh_ktab)) #endif /* Put in rp[n..2n-1] an approximation of the n high limbs @@ -151,7 +151,7 @@ mpfr_mullow_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp, static short sqrhigh_ktab[MPFR_SQRHIGH_TAB_SIZE]; #else static short sqrhigh_ktab[] = {MPFR_SQRHIGH_TAB}; -#define MPFR_SQRHIGH_TAB_SIZE (numberof (sqrhigh_ktab)) +#define MPFR_SQRHIGH_TAB_SIZE (numberof_const (sqrhigh_ktab)) #endif /* Put in rp[n..2n-1] an approximation of the n high limbs @@ -192,7 +192,7 @@ mpfr_sqrhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mp_size_t n) static short divhigh_ktab[MPFR_DIVHIGH_TAB_SIZE]; #else static short divhigh_ktab[] = {MPFR_DIVHIGH_TAB}; -#define MPFR_DIVHIGH_TAB_SIZE (numberof (divhigh_ktab)) +#define MPFR_DIVHIGH_TAB_SIZE (numberof_const (divhigh_ktab)) #endif #if !(defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_SBPI1_DIVAPPR_Q)) diff --git a/src/round_prec.c b/src/round_prec.c index 10baa55c1..bb4fa695b 100644 --- a/src/round_prec.c +++ b/src/round_prec.c @@ -71,8 +71,10 @@ mpfr_prec_round (mpfr_ptr x, mpfr_prec_t prec, mpfr_rnd_t rnd_mode) ow = MPFR_GET_ALLOC_SIZE(x); if (nw > ow) { + mpfr_size_limb_t *tmpx; + /* Realloc significand */ - mpfr_limb_ptr tmpx = (mpfr_limb_ptr) (*__gmp_reallocate_func) + tmpx = (mpfr_size_limb_t *) (*__gmp_reallocate_func) (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw)); MPFR_SET_MANT_PTR(x, tmpx); /* mant ptr must be set before alloc size */ diff --git a/src/set_dfl_prec.c b/src/set_dfl_prec.c index 748bf7fa4..36dd8bf23 100644 --- a/src/set_dfl_prec.c +++ b/src/set_dfl_prec.c @@ -23,15 +23,8 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" /* default is IEEE double precision, i.e. 53 bits */ -MPFR_THREAD_ATTR mpfr_prec_t __gmpfr_default_fp_bit_precision \ - = IEEE_DBL_MANT_DIG; - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -mpfr_prec_t * __gmpfr_default_fp_bit_precision_f() -{ - return &__gmpfr_default_fp_bit_precision; -} -#endif +MPFR_THREAD_VAR (mpfr_prec_t, __gmpfr_default_fp_bit_precision, + IEEE_DBL_MANT_DIG) void mpfr_set_default_prec (mpfr_prec_t prec) diff --git a/src/set_ld.c b/src/set_ld.c index 41256140b..9f63e5138 100644 --- a/src/set_ld.c +++ b/src/set_ld.c @@ -45,6 +45,12 @@ static const union { #define MPFR_LDBL_MAX LDBL_MAX #endif +/* To check for +inf, one can use the test x > MPFR_LDBL_MAX, as LDBL_MAX + is the maximum finite number representable in a long double, according + to DR 467; see + http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2092.htm + If this fails on some platform, a test x - x != 0 might be used. */ + #if defined(HAVE_LDOUBLE_IS_DOUBLE) /* the "long double" format is the same as "double" */ @@ -165,10 +171,6 @@ mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode) LONGDOUBLE_NAN_ACTION (d, goto nan); /* Check for INF */ - /* Note: according to the ISO C standard, there may be finite numbers - larger than LDBL_MAX, among the values that are not floating-point - numbers. If the following fails on some platform, a test d - d != 0 - could be used. */ if (d > MPFR_LDBL_MAX) { mpfr_set_inf (r, 1); @@ -227,10 +229,6 @@ mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode) LONGDOUBLE_NAN_ACTION (d, goto nan); /* Check for INF */ - /* Note: according to the ISO C standard, there may be finite numbers - larger than LDBL_MAX, among the values that are not floating-point - numbers. If the following fails on some platform, a test d - d != 0 - could be used. */ if (d > MPFR_LDBL_MAX) { mpfr_set_inf (r, 1); diff --git a/src/set_prec.c b/src/set_prec.c index bd2864dd9..eb2a721c7 100644 --- a/src/set_prec.c +++ b/src/set_prec.c @@ -26,7 +26,6 @@ MPFR_HOT_FUNCTION_ATTR void mpfr_set_prec (mpfr_ptr x, mpfr_prec_t p) { mp_size_t xsize, xoldsize; - mpfr_limb_ptr tmp; /* first, check if p is correct */ MPFR_ASSERTN (MPFR_PREC_COND (p)); @@ -38,8 +37,12 @@ mpfr_set_prec (mpfr_ptr x, mpfr_prec_t p) xoldsize = MPFR_GET_ALLOC_SIZE (x); if (xsize > xoldsize) { - tmp = (mpfr_limb_ptr) (*__gmp_reallocate_func) - (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(xoldsize), MPFR_MALLOC_SIZE(xsize)); + mpfr_size_limb_t *tmp; + + tmp = (mpfr_size_limb_t *) (*__gmp_reallocate_func) + (MPFR_GET_REAL_PTR(x), + MPFR_MALLOC_SIZE(xoldsize), + MPFR_MALLOC_SIZE(xsize)); MPFR_SET_MANT_PTR(x, tmp); MPFR_SET_ALLOC_SIZE(x, xsize); } diff --git a/src/set_rnd.c b/src/set_rnd.c index 10fc23231..bb96930da 100644 --- a/src/set_rnd.c +++ b/src/set_rnd.c @@ -23,15 +23,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" -MPFR_THREAD_ATTR mpfr_rnd_t __gmpfr_default_rounding_mode = MPFR_RNDN; - -#ifdef MPFR_WIN_THREAD_SAFE_DLL -mpfr_rnd_t * -__gmpfr_default_rounding_mode_f() -{ - return &__gmpfr_default_rounding_mode; -} -#endif +MPFR_THREAD_VAR (mpfr_rnd_t, __gmpfr_default_rounding_mode, MPFR_RNDN) void mpfr_set_default_rounding_mode (mpfr_rnd_t rnd_mode) diff --git a/src/sin_cos.c b/src/sin_cos.c index 37c113a25..0ecdcf98e 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -289,9 +289,9 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r, mpz_t T[GMP_NUMB_BITS], Q[GMP_NUMB_BITS], ptoj[GMP_NUMB_BITS], pp; mpfr_prec_t log2_nb_terms[GMP_NUMB_BITS], mult[GMP_NUMB_BITS]; mpfr_prec_t accu[GMP_NUMB_BITS], size_ptoj[GMP_NUMB_BITS]; - mpfr_prec_t prec_i_have, r0 = r, pp_s, p_s; - unsigned long alloc, i, j, k; - mpfr_prec_t l; + mpfr_prec_t prec_i_have, h, r0 = r, pp_s, p_s; + unsigned long i, j, m; + int alloc, k, l; if (MPFR_UNLIKELY(mpz_cmp_ui (p, 0) == 0)) /* sin(x)/x -> 1 */ { @@ -307,10 +307,10 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r, mpz_init (pp); /* normalize p (non-zero here) */ - l = mpz_scan1 (p, 0); - mpz_fdiv_q_2exp (pp, p, l); /* p = pp * 2^l */ + h = mpz_scan1 (p, 0); + mpz_fdiv_q_2exp (pp, p, h); /* p = pp * 2^h */ mpz_mul (pp, pp, pp); - r = 2 * (r - l); /* x^2 = (p/2^r0)^2 = pp / 2^r */ + r = 2 * (r - h); /* x^2 = (p/2^r0)^2 = pp / 2^r */ /* now p is odd */ alloc = 2; @@ -388,46 +388,45 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r, /* accumulate all products in T[0] and Q[0]. Warning: contrary to above, here we do not have log2_nb_terms[k-1] = log2_nb_terms[k]+1. */ - l = 0; /* number of accumulated terms in the right part T[k]/Q[k] */ + h = 0; /* number of accumulated terms in the right part T[k]/Q[k] */ while (k > 0) { - j = log2_nb_terms[k-1]; - mpz_mul (T[k], T[k], ptoj[j]); + mpz_mul (T[k], T[k], ptoj[log2_nb_terms[k-1]]); mpz_mul (T[k-1], T[k-1], Q[k]); - l += 1 << log2_nb_terms[k]; - mpz_mul_2exp (T[k-1], T[k-1], r * l); + h += (mpfr_prec_t) 1 << log2_nb_terms[k]; + mpz_mul_2exp (T[k-1], T[k-1], r * h); mpz_add (T[k-1], T[k-1], T[k]); mpz_mul (Q[k-1], Q[k-1], Q[k]); k--; } - l = r0 + r * (i - 1); /* implicit multiplier 2^r for Q0 */ - /* at this point T[0]/(2^l*Q[0]) is an approximation of sin(x) where the 1st + m = r0 + r * (i - 1); /* implicit multiplier 2^r for Q0 */ + /* at this point T[0]/(2^m*Q[0]) is an approximation of sin(x) where the 1st neglected term has contribution < 1/2^prec, thus since the series has alternate signs, the error is < 1/2^prec */ /* we truncate Q0 to prec bits: the relative error is at most 2^(1-prec), which means that Q0 = Q[0] * (1+theta) with |theta| <= 2^(1-prec) [up to a power of two] */ - l += reduce (Q0, Q[0], prec); - l -= reduce (T[0], T[0], prec); - /* multiply by x = p/2^l */ + m += reduce (Q0, Q[0], prec); + m -= reduce (T[0], T[0], prec); + /* multiply by x = p/2^m */ mpz_mul (S0, T[0], p); - l -= reduce (S0, S0, prec); /* S0 = T[0] * (1 + theta)^2 up to power of 2 */ + m -= reduce (S0, S0, prec); /* S0 = T[0] * (1 + theta)^2 up to power of 2 */ /* sin(X) ~ S0/Q0*(1 + theta)^3 + err with |theta| <= 2^(1-prec) and |err| <= 2^(-prec), thus since |S0/Q0| <= 1: |sin(X) - S0/Q0| <= 4*|theta*S0/Q0| + |err| <= 9*2^(-prec) */ mpz_clear (pp); - for (j = 0; j < alloc; j ++) + for (k = 0; k < alloc; k ++) { - mpz_clear (T[j]); - mpz_clear (Q[j]); - mpz_clear (ptoj[j]); + mpz_clear (T[k]); + mpz_clear (Q[k]); + mpz_clear (ptoj[k]); } /* compute cos(X) from sin(X): sqrt(1-(S/Q)^2) = sqrt(Q^2-S^2)/Q - = sqrt(Q0^2*2^(2l)-S0^2)/Q0. + = sqrt(Q0^2*2^(2m)-S0^2)/Q0. Write S/Q = sin(X) + eps with |eps| <= 9*2^(-prec), then sqrt(Q^2-S^2) = sqrt(Q^2-Q^2*(sin(X)+eps)^2) = sqrt(Q^2*cos(X)^2-Q^2*(2*sin(X)*eps+eps^2)) @@ -441,14 +440,14 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r, = Q*cos(X)*(1+eps3+eps2/(Q*cos(X))) = Q*cos(X)*(1+eps4) with |eps4| <= 9*2^(-prec) since |Q| >= 2^(prec-1) */ - /* we assume that Q0*2^l >= 2^(prec-1) */ - MPFR_ASSERTN(l + mpz_sizeinbase (Q0, 2) >= prec); + /* we assume that Q0*2^m >= 2^(prec-1) */ + MPFR_ASSERTN(m + mpz_sizeinbase (Q0, 2) >= prec); mpz_mul (C0, Q0, Q0); - mpz_mul_2exp (C0, C0, 2 * l); + mpz_mul_2exp (C0, C0, 2 * m); mpz_submul (C0, S0, S0); mpz_sqrt (C0, C0); - return l; + return m; } /* Put in s and c approximations of sin(x) and cos(x) respectively. diff --git a/src/sqrt.c b/src/sqrt.c index 332720e3f..3dcf408da 100644 --- a/src/sqrt.c +++ b/src/sqrt.c @@ -216,7 +216,8 @@ mpn_rsqrtrem1 (mp_limb_t a0) mp_limb_t a = a0 >> (GMP_NUMB_BITS - 4); mp_limb_t b = (a0 >> (GMP_NUMB_BITS - 8)) & 0xf; mp_limb_t c = (a0 >> (GMP_NUMB_BITS - 12)) & 0xf; - mp_limb_t x0, a1, t; + mp_limb_t x0, t; + MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 32 || GMP_NUMB_BITS == 64); @@ -243,19 +244,23 @@ mpn_rsqrtrem1 (mp_limb_t a0) return x0 - 4; /* ensures x0 <= 2^36/sqrt(a0) */ #else /* GMP_NUMB_BITS = 64 */ - a1 = a0 >> (GMP_NUMB_BITS - 32); - /* a1 has 32 bits, thus a1*x0^2 has 64 bits */ - t = (mp_limb_signed_t) (-a1 * (x0 * x0)) >> 32; - /* t has 32 bits now */ + { + mp_limb_t a1 = a0 >> (GMP_NUMB_BITS - 32); + /* a1 has 32 bits, thus a1*x0^2 has 64 bits */ + t = (mp_limb_signed_t) (-a1 * (x0 * x0)) >> 32; + /* t has 32 bits now */ + } x0 = (x0 << 15) + ((mp_limb_signed_t) (x0 * t) >> (17+1)); /* now x0 is a 31-bit approximation (32 bits, 1 <= x0/2^31 <= 2), with maximal error 2^(-19.19) */ - a1 = a0 >> (GMP_NUMB_BITS - 40); /* a1 has 40 bits */ - t = (x0 * x0) >> 22; /* t has 40 bits */ - /* a1 * t has 80 bits, but we know the upper 19 bits cancel with 1 */ - t = (mp_limb_signed_t) (-a1 * t) >> 31; /* it remains 49 bits in theory, - but t has only 31 bits at most */ + { + mp_limb_t a1 = a0 >> (GMP_NUMB_BITS - 40); /* a1 has 40 bits */ + t = (x0 * x0) >> 22; /* t has 40 bits */ + /* a1 * t has 80 bits, but we know the upper 19 bits cancel with 1 */ + t = (mp_limb_signed_t) (-a1 * t) >> 31; + /* it remains 49 bits in theory, but t has only 31 bits at most */ + } x0 = (x0 << 9) + ((mp_limb_signed_t) (x0 * t) >> (31+1+49-40)); /* now x0 is a 1+40-bit approximation, @@ -413,7 +418,9 @@ mpn_sqrtrem2 (mpfr_limb_ptr sp, mpfr_limb_ptr rp, mpfr_limb_srcptr np) return x; } -/* Special code for prec(r), prec(u) < GMP_NUMB_BITS. */ +/* Special code for prec(r), prec(u) < GMP_NUMB_BITS. We cannot have + prec(u) = GMP_NUMB_BITS here, since when the exponent of u is odd, + we need to shift u by one bit to the right without losing any bit. */ static int mpfr_sqrt1 (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) { @@ -451,7 +458,7 @@ mpfr_sqrt1 (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) if (exp_r > __gmpfr_emax) return mpfr_overflow (r, rnd_mode, 1); - /* See comments in mpfr_divsp1 */ + /* See comments in mpfr_div_1 */ if (exp_r < __gmpfr_emin) { if (rnd_mode == MPFR_RNDN) @@ -557,6 +564,71 @@ mpn_rsqrtrem2 (mpfr_limb_ptr xp, mpfr_limb_srcptr ap) xp[1] = xp[1] >> 24; } +#if 0 +/* this code seems to be slightly slower than the one below, for precision + 113 bits on a 64-bit machine, and in addition it requires GMP internals + (for __gmpn_invert_limb) */ +static mp_limb_t +mpn_sqrtrem4 (mpfr_limb_ptr sp, mpfr_limb_ptr rp, mpfr_limb_srcptr ap) +{ + mp_limb_t x, r, t; + + x = mpn_sqrtrem2 (sp + 1, rp + 1, ap + 2); + + /* now a1 = s1^2 + x*B + r1, with a1 = {ap+2,2}, s1 = sp[1], r1 = rp[1], + with 0 <= x*B + r1 <= 2*s1 */ + + /* use Newton's iteration for the square root, x' = 1/2 * (x + a/x), which + rewrites as: + + x' = x + (a - x^2)/(2x) + + Since x' - sqrt(a) = (x-sqrt(a))^2/(2x), we have x' >= sqrt(a). + If we round (a - x^2)/(2x) downwards, then we will still get + an upper approximation of floor(sqrt(a)). */ + + t = rp[1] << (GMP_NUMB_BITS - 1); + rp[1] = (x << (GMP_NUMB_BITS - 1)) | (rp[1] >> 1); /* rp[1] = floor((a-x^2)/2) */ + /* Since x*B + r1 <= 2*s1, we now have rp[1] <= s1, and since __udiv_qrnnd_preinv + requires its 3rd argument to be smaller than its 5th argument, we must + distinguish the case rp[1] == s1. */ + if (MPFR_UNLIKELY(rp[1] == sp[1])) + { + /* Necessarily t=0 since x*B + r1 <= 2*s1 initially. + Taking sp[0] = 2^GMP_NUMB_BITS (if representable) would be too large + since it would mean sp[1] = sp[1]+1, and the remainder rp[1] would + become negative. */ + sp[0] = ~MPFR_LIMB_ZERO; + r = sp[1]; + } + else + { + r = __gmpn_invert_limb (sp[1]); + __udiv_qrnnd_preinv (sp[0], r, rp[1], t, sp[1], r); + } + + /* now sp[0] = floor((a-x^2)/(2x)) = floor((x*B+r1)/2/s1) */ + + /* r1*B = 2*s1*s0 + 2*r + a - s^2 = a1*B^2 + a0 - (s1*B+s0)^2 + = (s1^2+r1)*B^2 + a0 - s1^2*B^2 - 2*s1*s0*B - s0^2 + = r1*B^2 + a0 - 2*s1*s0*B - s0^2 + = 2*r*B + a0 - s0^2 */ + umul_ppmm (rp[1], rp[0], sp[0], sp[0]); /* s0^2 */ + t = -mpn_sub_n (rp, ap, rp, 2); /* a0 - s0^2 */ + t += 2 * (r >> (GMP_NUMB_BITS - 1)); + r <<= 1; + rp[1] += r; + t += rp[1] < r; + if ((mp_limb_signed_t) t < 0) + { + mpn_sub_1 (sp, sp, 2, 1); + t += mpn_addmul_1 (rp, sp, 2, 2); + t += mpn_add_1 (rp, rp, 2, 1); + } + return t; +} +#else /* Given as input ap[0-3], with B/4 <= ap[3] (where B = 2^GMP_NUMB_BITS), mpn_sqrtrem4 returns a value x, 0 <= x <= 1, and stores values s in sp[0-1] and r in rp[0-1] such that: @@ -565,7 +637,8 @@ mpn_rsqrtrem2 (mpfr_limb_ptr xp, mpfr_limb_srcptr ap) or equivalently x*B + r <= 2*s. - This code currently assumes GMP_NUMB_BITS = 64. */ + This code currently assumes GMP_NUMB_BITS = 64, and takes on average + 135.68 cycles on an Intel i5-6500 for precision 113 bits */ static mp_limb_t mpn_sqrtrem4 (mpfr_limb_ptr sp, mpfr_limb_ptr rp, mpfr_limb_srcptr ap) { @@ -668,12 +741,13 @@ mpn_sqrtrem4 (mpfr_limb_ptr sp, mpfr_limb_ptr rp, mpfr_limb_srcptr ap) mpn_add_1 (sp, sp, 2, 1); /* add 2 to t0*B + {b, 2}: t0 += mpn_add_1 (b, b, 2, 2) */ b[0] += 2; - if (b[0] < 2) - t0 += (b[1]++ == 0); + if (b[0] == 0) + t0 += (++b[1] == 0); } return b[2]; } +#endif /* Special code for GMP_NUMB_BITS < prec(r) < 2*GMP_NUMB_BITS, and GMP_NUMB_BITS < prec(u) <= 2*GMP_NUMB_BITS. @@ -690,18 +764,20 @@ mpfr_sqrt2 (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) if (((unsigned int) exp_u & 1) != 0) { np[3] = up[1] >> 1; - np[2] = (up[1] << (GMP_NUMB_BITS-1)) | (up[0] >> 1); + np[2] = (up[1] << (GMP_NUMB_BITS - 1)) | (up[0] >> 1); + np[1] = up[0] << (GMP_NUMB_BITS - 1); exp_u ++; } else { np[3] = up[1]; np[2] = up[0]; + np[1] = 0; } MPFR_ASSERTD (((unsigned int) exp_u & 1) == 0); exp_r = exp_u / 2; - np[1] = np[0] = 0; + np[0] = 0; sb = mpn_sqrtrem4 (rp, tp, np); sb |= tp[0] | tp[1]; rb = rp[0] & (MPFR_LIMB_ONE << (sh - 1)); @@ -713,7 +789,7 @@ mpfr_sqrt2 (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) if (exp_r > __gmpfr_emax) return mpfr_overflow (r, rnd_mode, 1); - /* See comments in mpfr_divsp1 */ + /* See comments in mpfr_div_1 */ if (exp_r < __gmpfr_emin) { if (rnd_mode == MPFR_RNDN) @@ -642,17 +642,6 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, (mpfr_eexp_t) maxexp, maxexp == MPFR_EXP_MIN ? " (MPFR_EXP_MIN)" : "")); - /* Let's copy/shift the bits [max(u,minexp),e) to the - most significant part of the destination, and zero - the least significant part (there can be one only if - u < minexp). The trailing bits of the destination may - contain garbage at this point. Then, at the same time, - take the absolute value and do an initial rounding, - zeroing the trailing bits at this point. - TODO: This may be improved by merging some operations - in particular cases. The average speed-up may not be - significant, though. To be tested... */ - sn = MPFR_PREC2LIMBS (sq); sd = (mpfr_prec_t) sn * GMP_NUMB_BITS - sq; sh = cancel % GMP_NUMB_BITS; @@ -673,7 +662,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, /* Determine the rounding bit, which is represented. */ td = tq % GMP_NUMB_BITS; - lbit = (wp[wi] >> td) & 1; + lbit = (wp[wi] >> td) & MPFR_LIMB_ONE; rbit = td >= 1 ? ((wp[wi] >> (td - 1)) & MPFR_LIMB_ONE) : (MPFR_ASSERTD (wi >= 1), wp[wi-1] >> (GMP_NUMB_BITS - 1)); MPFR_ASSERTD (rbit == 0 || rbit == 1); @@ -685,7 +674,14 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, inex = 0 if the final sum is exact, else 1, i.e. inex = rounding bit || sticky bit. In round to nearest, also determine the rounding direction: obtained from - the rounding bit possibly except in halfway cases. */ + the rounding bit possibly except in halfway cases. + Halfway cases are rounded toward -inf iff the last bit + of the truncated significand in two's complement is 0 + (in precision > 1, because the parity after rounding is + the same in two's complement and sign + magnitude; in + precision 1, one checks that the rule works for both + positive (lbit == 1) and negative (lbit == 0) numbers, + rounding halfway cases away from zero). */ if (MPFR_LIKELY (rbit == 0 || (rnd == MPFR_RNDN && lbit == 0))) { /* We need to determine the sticky bit, either to set inex @@ -824,14 +820,14 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, * not taken into account, and if it occurs, this is * necessarily on a machine number (-> tmd = 1). */ - lbit = u == minexp ? wp[0] & 1 : 0; + lbit = u == minexp ? wp[0] & MPFR_LIMB_ONE : 0; rbit = 0; inex = tmd = maxexp != MPFR_EXP_MIN; } MPFR_ASSERTD (rbit == 0 || rbit == 1); - /* neg = 0 if negative, 1 if positive. */ + /* neg = 1 if negative, 0 if positive. */ neg = wp[ws-1] >> (GMP_NUMB_BITS - 1); MPFR_ASSERTD (neg == 0 || neg == 1); @@ -947,7 +943,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, zz = zs - wi; MPFR_ASSERTD (zz >= 0 && zz < zs); if (zz > 0) - MPN_COPY_DECR (zp + zz, wp, wi); + MPN_COPY (zp + zz, wp, wi); } /* Compute minexp2 = minexp - (zs * GMP_NUMB_BITS + td) safely. */ @@ -1026,6 +1022,12 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_SIGN (sum) = sgn; + /* Let's copy/shift the bits [max(u,minexp),e) to the + most significant part of the destination, and zero + the least significant part (there can be one only if + u < minexp). The trailing bits of the destination may + contain garbage at this point. */ + if (MPFR_LIKELY (u > minexp)) { mp_size_t wi; @@ -1062,6 +1064,11 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPN_ZERO (sump, sn - en); } + /* Let's take the complement if the result is negative, and at + the same time, do the rounding and zero the trailing bits. + As this is valid only for precisions >= 2, there is special + code for precision 1 first. */ + if (MPFR_UNLIKELY (sq == 1)) /* precision 1 */ { sump[0] = MPFR_LIMB_HIGHBIT; diff --git a/src/uceil_log2.c b/src/uceil_log2.c index 7ee5166ae..8e1aa1f83 100644 --- a/src/uceil_log2.c +++ b/src/uceil_log2.c @@ -24,7 +24,8 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., /* returns ceil(log(d)/log(2)) if d > 0, -1023 if d = +0, - and floor(log(-d)/log(2))+1 if d < 0*/ + and floor(log(-d)/log(2))+1 if d < 0 +*/ long __gmpfr_ceil_log2 (double d) { @@ -36,6 +37,7 @@ __gmpfr_ceil_log2 (double d) /* The cast below is useless in theory, but let us not depend on the integer promotion rules (for instance, tcc is currently wrong). */ exp = (long) x.s.exp - 1023; + MPFR_ASSERTN (exp < 1023); /* fail on infinities */ x.s.exp = 1023; /* value for 1 <= d < 2 */ if (x.d != 1.0) /* d: not a power of two? */ exp++; diff --git a/src/ufloor_log2.c b/src/ufloor_log2.c index fb27b6e42..32cb486f0 100644 --- a/src/ufloor_log2.c +++ b/src/ufloor_log2.c @@ -26,26 +26,30 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., long __gmpfr_floor_log2 (double d) { + long exp; #if _MPFR_IEEE_FLOATS union mpfr_ieee_double_extract x; x.d = d; - return (long) x.s.exp - 1023; + /* The cast below is useless in theory, but let us not depend on the + integer promotion rules (for instance, tcc is currently wrong). */ + exp = (long) x.s.exp - 1023; + MPFR_ASSERTN (exp < 1023); /* fail on infinities */ + return exp; #else /* _MPFR_IEEE_FLOATS */ - long exp; double m; MPFR_ASSERTD (d >= 0); if (d >= 1.0) { exp = -1; - for( m= 1.0 ; m <= d ; m *=2.0 ) + for (m = 1.0; m <= d; m *= 2.0) exp++; } else { exp = 0; - for( m= 1.0 ; m > d ; m *= (1.0/2.0) ) + for (m = 1.0; m > d; m *= 0.5) exp--; } return exp; diff --git a/src/zeta.c b/src/zeta.c index e76716ce8..164a4fdf8 100644 --- a/src/zeta.c +++ b/src/zeta.c @@ -414,6 +414,10 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */ add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); prec1 = precz + add; + /* FIXME: to avoid that the working precision (prec1) depends on the + input precision, one would need to take into account the error made + when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1) + below, and also in the case y=Inf (i.e. when gamma(s1) overflows). */ prec1 = MAX (prec1, precs1) + 10; MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p); diff --git a/tests/tabs.c b/tests/tabs.c index 994f417d5..a63e2b15a 100644 --- a/tests/tabs.c +++ b/tests/tabs.c @@ -162,7 +162,6 @@ check_cmp (int argc, char *argv[]) int main (int argc, char *argv[]) { - mpfr_test_init (); tests_start_mpfr (); check_inexact (); diff --git a/tests/tdiv.c b/tests/tdiv.c index 39cab385f..5239b9153 100644 --- a/tests/tdiv.c +++ b/tests/tdiv.c @@ -645,7 +645,8 @@ check_inexact (void) mpfr_set_prec (u, 2); mpfr_set_str_binary (u, "0.1E0"); mpfr_set_prec (y, 28); - if ((inexact = test_div (y, x, u, MPFR_RNDN) >= 0)) + inexact = test_div (y, x, u, MPFR_RNDN); + if (inexact >= 0) { printf ("Wrong inexact flag (1): expected -1, got %d\n", inexact); @@ -657,7 +658,8 @@ check_inexact (void) mpfr_set_prec (u, 15); mpfr_set_str_binary (u, "0.101101000001100E-1"); mpfr_set_prec (y, 92); - if ((inexact = test_div (y, x, u, MPFR_RNDN)) <= 0) + inexact = test_div (y, x, u, MPFR_RNDN); + if (inexact <= 0) { printf ("Wrong inexact flag for rnd=MPFR_RNDN(1): expected 1, got %d\n", inexact); diff --git a/tests/tget_d.c b/tests/tget_d.c index 20d4e8672..32e4d1022 100644 --- a/tests/tget_d.c +++ b/tests/tget_d.c @@ -89,7 +89,7 @@ static void check_inf_nan (void) { /* only if nans and infs are available */ -#if _GMP_IEEE_FLOATS && !defined(MPFR_ERRDIVZERO) +#if !defined(MPFR_ERRDIVZERO) mpfr_t x; double d; @@ -139,7 +139,7 @@ check_max (void) MPFR_ASSERTN(d == -DBL_MAX); d = mpfr_get_d (u, MPFR_RNDU); MPFR_ASSERTN(d == -DBL_MAX); -#if _GMP_IEEE_FLOATS && !defined(MPFR_ERRDIVZERO) +#if !defined(MPFR_ERRDIVZERO) d = mpfr_get_d (u, MPFR_RNDN); MPFR_ASSERTN(DOUBLE_ISINF(d) && d < 0.0); d = mpfr_get_d (u, MPFR_RNDD); @@ -151,7 +151,7 @@ check_max (void) MPFR_ASSERTN(d == DBL_MAX); d = mpfr_get_d (u, MPFR_RNDD); MPFR_ASSERTN(d == DBL_MAX); -#if _GMP_IEEE_FLOATS && !defined(MPFR_ERRDIVZERO) +#if !defined(MPFR_ERRDIVZERO) d = mpfr_get_d (u, MPFR_RNDN); MPFR_ASSERTN(DOUBLE_ISINF(d) && d > 0.0); d = mpfr_get_d (u, MPFR_RNDU); diff --git a/tests/tget_d_2exp.c b/tests/tget_d_2exp.c index 34a9f404f..9e8d6957a 100644 --- a/tests/tget_d_2exp.c +++ b/tests/tget_d_2exp.c @@ -78,8 +78,7 @@ check_round (void) static void check_inf_nan (void) { - /* only if nans and infs are available */ -#if _GMP_IEEE_FLOATS +#if !defined(MPFR_ERRDIVZERO) mpfr_t x; double d; long exp; diff --git a/tests/tget_ld_2exp.c b/tests/tget_ld_2exp.c index 9e7dff838..af70af2cb 100644 --- a/tests/tget_ld_2exp.c +++ b/tests/tget_ld_2exp.c @@ -77,8 +77,7 @@ check_round (void) static void check_inf_nan (void) { - /* only if nans and infs are available */ -#if _GMP_IEEE_FLOATS +#if !defined(MPFR_ERRDIVZERO) mpfr_t x; double d; long exp; diff --git a/tests/tget_set_d64.c b/tests/tget_set_d64.c index 65c50b6f7..17c18a02f 100644 --- a/tests/tget_set_d64.c +++ b/tests/tget_set_d64.c @@ -34,11 +34,11 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., /* #define DEBUG */ -#if _GMP_IEEE_FLOATS +#if _MPFR_IEEE_FLOATS static void print_decimal64 (_Decimal64 d) { - union ieee_double_extract x; + union mpfr_ieee_double_extract x; union ieee_double_decimal64 y; unsigned int Gh, i; @@ -65,7 +65,7 @@ print_decimal64 (_Decimal64 d) { printf ("%.15Le\n", (long double) d); } -#endif /* _GMP_IEEE_FLOATS */ +#endif /* _MPFR_IEEE_FLOATS */ #define ERR_MISC(V) \ do \ diff --git a/tests/tset_ld.c b/tests/tset_ld.c index d974ba98e..fd7959f5b 100644 --- a/tests/tset_ld.c +++ b/tests/tset_ld.c @@ -58,10 +58,11 @@ Isnan_ld (long double d) /* Return the minimal number of bits to represent d exactly (0 for zero). If flag is non-zero, also print d. */ +/* FIXME: This function doesn't work if the rounding precision is reduced. */ static mpfr_prec_t print_binary (long double d, int flag) { - long double e; + long double e, f, r; long exp = 1; mpfr_prec_t prec = 0; @@ -96,7 +97,7 @@ print_binary (long double d, int flag) printf ("0.0\n"); return prec; } - /* now d > 0 */ + MPFR_ASSERTN (d > 0); e = (long double) 1.0; while (e > d) { @@ -104,26 +105,49 @@ print_binary (long double d, int flag) exp --; } if (flag == 2) printf ("1: e=%.36Le\n", e); - /* now d >= e */ - while (d >= e + e) + MPFR_ASSERTN (d >= e); + /* FIXME: There can be an overflow here, which may not be supported + on all platforms. */ + while (f = e + e, d >= f) { - e = e + e; + e = f; exp ++; } if (flag == 2) printf ("2: e=%.36Le\n", e); - /* now e <= d < 2e */ + MPFR_ASSERTN (e <= d && d < f); if (flag == 1) printf ("0."); if (flag == 2) printf ("3: d=%.36Le e=%.36Le prec=%ld\n", d, e, (long) prec); - while (d > (long double) 0.0) + /* Note: the method we use here to extract the bits of d is the following, + to deal with the case where the rounding precision is less than the + precision of d: + (1) we accumulate the upper bits of d into f + (2) when accumulating a new bit into f is not exact, we subtract + f from d and reset f to 0 + This is guaranteed to work only when the rounding precision is at least + half the precision of d, since otherwise d-f might not be exact. + This method does not work with flush-to-zero on underflow. */ + f = 0.0; /* will hold accumulated powers of 2 */ + while (1) { prec++; - if (d >= e) + r = f + e; + /* r is close to f (in particular in the cases where f+e may + not be exact), so that r - f should be exact. */ + if (r - f != e) /* f+e is not exact */ + { + d -= f; /* should be exact */ + f = 0.0; + r = e; + } + if (d >= r) { if (flag == 1) printf ("1"); - d = (long double) ((long double) d - (long double) e); + if (d == r) + break; + f = r; } else { @@ -131,6 +155,7 @@ print_binary (long double d, int flag) printf ("0"); } e *= (long double) 0.5; + MPFR_ASSERTN (e != 0); /* may fail with flush-to-zero on underflow */ if (flag == 2) printf ("4: d=%.36Le e=%.36Le prec=%ld\n", d, e, (long) prec); } @@ -139,9 +164,9 @@ print_binary (long double d, int flag) return prec; } -/* checks that a long double converted to a mpfr (with precision - MPFR_LDBL_MANT_DIG), then converted back to a long double gives the initial - value, or in other words mpfr_get_ld(mpfr_set_ld(d)) = d. +/* Checks that a long double converted exactly to a MPFR number, then + converted back to a long double gives the initial value, or in other + words, mpfr_get_ld(mpfr_set_ld(d)) = d. */ static void check_set_get (long double d) @@ -152,6 +177,7 @@ check_set_get (long double d) long double e; int inex; + /* Select a precision to ensure that the conversion of d to x be exact. */ prec = print_binary (d, 0); if (prec < MPFR_PREC_MIN) prec = MPFR_PREC_MIN; @@ -184,6 +210,7 @@ check_set_get (long double d) printf (" x = "); mpfr_dump (x); printf (" MPFR_LDBL_MANT_DIG=%u\n", MPFR_LDBL_MANT_DIG); + printf (" prec=%lu\n", prec); print_binary (d, 2); exit (1); } @@ -321,12 +348,39 @@ check_subnormal (void) printf ("d=%Le\n", d); printf ("x="); mpfr_dump (x); printf ("e=%Le\n", e); + exit (1); } d *= 0.5; } mpfr_clear (x); } +static void +check_overflow (void) +{ + long double d, e; + mpfr_t x; + int i; + + mpfr_init2 (x, MPFR_LDBL_MANT_DIG); + for (i = 0; i < 2; i++) + { + d = i == 0 ? LDBL_MAX : -LDBL_MAX; + mpfr_set_ld (x, d, MPFR_RNDN); + mpfr_mul_2ui (x, x, 1, MPFR_RNDN); + e = mpfr_get_ld (x, MPFR_RNDN); + if (! DOUBLE_ISINF (e) || (i == 0 ? (e < 0) : (e > 0))) + { + printf ("Error in check_overflow.\n"); + printf ("d=%Le\n", d); + printf ("x="); mpfr_dump (x); + printf ("e=%Le\n", e); + exit (1); + } + } + mpfr_clear (x); +} + /* issue reported by Sisyphus on powerpc */ static void test_20140212 (void) @@ -364,10 +418,72 @@ test_20140212 (void) mpfr_clear (fr2); } +/* bug reported by Walter Mascarenhas + https://sympa.inria.fr/sympa/arc/mpfr/2016-09/msg00005.html */ +static void +bug_20160907 (void) +{ +#if HAVE_LDOUBLE_IEEE_EXT_LITTLE + long double dn, ld; + mpfr_t mp; + long e; + mpfr_long_double_t x; + + /* the following is the encoding of the smallest subnormal number + for HAVE_LDOUBLE_IEEE_EXT_LITTLE */ + x.s.manl = 1; + x.s.manh = 0; + x.s.expl = 0; + x.s.exph = 0; + x.s.sign= 0; + dn = x.ld; + e = -16445; + /* dn=2^e is now the smallest subnormal. */ + + mpfr_init2 (mp, 64); + mpfr_set_ui_2exp (mp, 1, e - 1, MPFR_RNDN); + ld = mpfr_get_ld (mp, MPFR_RNDU); + /* since mp = 2^(e-1) and ld is rounded upwards, we should have + ld = 2^e */ + if (ld != dn) + { + printf ("Error, ld = %Le <> dn = %Le\n", ld, dn); + printf ("mp="); + mpfr_out_str (stdout, 10, 0, mp, MPFR_RNDN); + printf ("\n"); + printf ("mp="); mpfr_dump (mp); + exit (1); + } + + /* check a few more numbers */ + for (e = -16446; e <= -16381; e++) + { + mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); + ld = mpfr_get_ld (mp, MPFR_RNDU); + mpfr_set_ld (mp, ld, MPFR_RNDU); + /* mp is 2^e rounded up, thus should be >= 2^e */ + MPFR_ASSERTN(mpfr_cmp_ui_2exp (mp, 1, e) >= 0); + + mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); + ld = mpfr_get_ld (mp, MPFR_RNDD); + mpfr_set_ld (mp, ld, MPFR_RNDD); + /* mp is 2^e rounded down, thus should be <= 2^e */ + if (mpfr_cmp_ui_2exp (mp, 3, e) > 0) + { + printf ("Error, expected value <= 2^%ld\n", e); + printf ("got "); mpfr_dump (mp); + exit (1); + } + } + + mpfr_clear (mp); +#endif +} + int main (int argc, char *argv[]) { - volatile long double d, e; + volatile long double d, e, maxp2; mpfr_t x; int i; mpfr_exp_t emax; @@ -383,12 +499,11 @@ main (int argc, char *argv[]) #endif tests_start_mpfr (); + mpfr_test_init (); check_gcc33_bug (); test_fixed_bugs (); - mpfr_test_init (); - mpfr_init2 (x, MPFR_LDBL_MANT_DIG + 64); #if !defined(MPFR_ERRDIVZERO) @@ -408,7 +523,7 @@ main (int argc, char *argv[]) mpfr_set_ld (x, DBL_NEG_ZERO, MPFR_RNDN); if (MPFR_IS_POS (x)) { -#if _GMP_IEEE_FLOATS +#if defined(HAVE_SIGNEDZ) printf ("Error: sign of -0.0 is not set correctly\n"); exit (1); #else @@ -430,22 +545,20 @@ main (int argc, char *argv[]) #endif /* check the largest power of two */ - d = 1.0; while (d < LDBL_MAX / 2.0) d += d; - check_set_get (d); - check_set_get (-d); + maxp2 = 1.0; + while (maxp2 < LDBL_MAX / 2.0) + maxp2 *= 2.0; + check_set_get (maxp2); + check_set_get (-maxp2); - /* check LDBL_MAX; according to the C standard, LDBL_MAX must be - exactly (1-b^(-LDBL_MANT_DIG)).b^LDBL_MAX_EXP, where b is the - radix (in practice, b = 2), even though there can be larger - long double values, not regarded as being in the subset of - the floating-point values of the system. As a consequence - (assuming b = 2), LDBL_MAX must be exactly representable on - LDBL_MANT_DIG bits. GCC is currently buggy[*], but LDBL_MAX - is still representable on LDBL_MANT_DIG bits. - [*] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=61399 */ d = LDBL_MAX; - check_set_get (d); - check_set_get (-d); + e = d / 2.0; + if (e != maxp2) /* false under NetBSD/x86 */ + { + /* d = LDBL_MAX does not have excess precision. */ + check_set_get (d); + check_set_get (-d); + } /* check the smallest power of two */ d = 1.0; @@ -454,14 +567,38 @@ main (int argc, char *argv[]) check_set_get (d); check_set_get (-d); - /* check that 2^i, 2^i+1 and 2^i-1 are correctly converted */ + /* check that 2^i, 2^i+1, 2^i-1 and 2^i-2^(i-2)-1 are correctly converted */ d = 1.0; for (i = 1; i < MPFR_LDBL_MANT_DIG + 8; i++) { d = 2.0 * d; /* d = 2^i */ check_set_get (d); - check_set_get (d + 1.0); - check_set_get (d - 1.0); + if (d + 1.0 != d) + check_set_get (d + 1.0); + else + { + mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN); + mpfr_add_ui (x, x, 1, MPFR_RNDN); + e = mpfr_get_ld (x, MPFR_RNDN); + check_set_get (e); + } + if (d - 1.0 != d) + check_set_get (d - 1.0); + else + { + mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN); + mpfr_sub_ui (x, x, 1, MPFR_RNDN); + e = mpfr_get_ld (x, MPFR_RNDN); + check_set_get (e); + } + if (i < 3) + continue; + /* The following test triggers a failure in r10844 for i = 56, + with gcc -mpc64 on x86 (64-bit ABI). */ + mpfr_set_ui_2exp (x, 3, i-2, MPFR_RNDN); + mpfr_sub_ui (x, x, 1, MPFR_RNDN); + e = mpfr_get_ld (x, MPFR_RNDN); + check_set_get (e); } for (i = 0; i < 10000; i++) @@ -488,8 +625,10 @@ main (int argc, char *argv[]) test_small (); check_subnormal (); + check_overflow (); test_20140212 (); + bug_20160907 (); tests_end_mpfr (); diff --git a/tests/tsqrt.c b/tests/tsqrt.c index 533d5420d..1963236b0 100644 --- a/tests/tsqrt.c +++ b/tests/tsqrt.c @@ -597,6 +597,25 @@ bug20160120 (void) mpfr_clear(y); } +/* Bug in mpfr_sqrt2 when prec(u) = 2*GMP_NUMB_BITS and the exponent of u is + odd: the last bit of u is lost. */ +static void +bug20160908 (void) +{ + mpfr_t r, u; + int n = GMP_NUMB_BITS, ret; + + mpfr_init2 (r, 2*n - 1); + mpfr_init2 (u, 2 * n); + mpfr_set_ui_2exp (u, 1, 2*n-2, MPFR_RNDN); /* u=2^(2n-2) with exp(u)=2n-1 */ + mpfr_nextabove (u); + /* now u = 2^(2n-2) + 1/2 */ + ret = mpfr_sqrt (r, u, MPFR_RNDZ); + MPFR_ASSERTN(ret == -1 && mpfr_cmp_ui_2exp (r, 1, n-1) == 0); + mpfr_clear (r); + mpfr_clear (u); +} + static void testall_rndf (mpfr_prec_t pmax) { @@ -773,6 +792,7 @@ main (void) bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 8, -256, 255, 4, 128, 800, 50); bug20160120 (); + bug20160908 (); tests_end_mpfr (); return 0; diff --git a/tests/tzeta.c b/tests/tzeta.c index 9fb5a640e..596e98824 100644 --- a/tests/tzeta.c +++ b/tests/tzeta.c @@ -415,8 +415,9 @@ main (int argc, char *argv[]) mpfr_clear (y); mpfr_clear (z); - /* FIXME: change the last argument back to 5 once the mpfr_zeta issue - has been found (see TODO). */ + /* FIXME: change the last argument back to 5 once the working precision + in the mpfr_zeta implementation no longer depends on the precision of + the input. */ test_generic (MPFR_PREC_MIN, 70, 1); test2 (); diff --git a/tools/announce-text b/tools/announce-text index d70329592..060745fd4 100755 --- a/tools/announce-text +++ b/tools/announce-text @@ -210,10 +210,10 @@ http://www.mpfr.org/mpfr-$base/mpfr-$vers.tar.bz2.asc http://www.mpfr.org/mpfr-$base/mpfr-$vers.tar.gz.asc http://www.mpfr.org/mpfr-$base/mpfr-$vers.zip.asc -Each tarball is signed by Vincent Lefèvre. This can be verified -using the DSA key ID 98C3739D; this key can be retrieved with: +Each tarball is signed by Vincent Lefèvre. This can be verified using +the DSA key ID 980C197698C3739D; this key can be retrieved with: - gpg --recv-keys 98C3739D + gpg --recv-keys 980C197698C3739D or by downloading it from <https://www.vinc17.net/pgp.html>. The key fingerprint is: diff --git a/tools/mbench/mfv5-mpfr.cc b/tools/mbench/mfv5-mpfr.cc index c36ed029d..3d667429b 100644 --- a/tools/mbench/mfv5-mpfr.cc +++ b/tools/mbench/mfv5-mpfr.cc @@ -254,6 +254,12 @@ public: return mpfr_atanh (a,b,r); } }; +class mpfr_pow_test { +public: + int func (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t r) { + return mpfr_pow (a,b,c,r); + } +}; class mpfr_get_ld_test { long double ld; public: @@ -295,6 +301,7 @@ static mpfr_test<mpfr_tanh_test> testCh ("mpfr_tanh"); static mpfr_test<mpfr_acosh_test> testDh ("mpfr_acosh"); static mpfr_test<mpfr_asinh_test> testEh ("mpfr_asinh"); static mpfr_test<mpfr_atanh_test> testFh ("mpfr_atanh"); +static mpfr_test<mpfr_pow_test> testGh ("mpfr_pow"); static mpfr_test<mpfr_get_ld_test> testFj ("mpfr_get_ld"); diff --git a/tools/mpfrlint b/tools/mpfrlint index 3cd884770..53293f896 100755 --- a/tools/mpfrlint +++ b/tools/mpfrlint @@ -19,7 +19,7 @@ fi err-if-output() { local dir msg checkoutput=1 checkstatus=1 output st h line - while [[ $1 = -* ]] + while [[ $1 == -* ]] do case $1 in --dir=*) @@ -61,9 +61,11 @@ err-if-output() fi set -e [[ -n $dir ]] && popd + # Do not fail. + true } -if [[ $1 = test ]] then +if [[ $1 == test ]] then export LC_ALL=C err-if-output Err-e exit 17 err-if-output Err-f false @@ -267,7 +269,7 @@ grep -E 'mpfr_set_(inf|zero) *\([^,]*, *[-+]?([02-9]|1[^)])' $srctests | \ # such as MS Windows use a config.h, which is included by mpfr-impl.h). for file in src/*.c do - [[ "$file" = src/(jyn_asympt|mini-gmp|round_raw_generic).c ]] || \ + [[ "$file" == src/(jyn_asympt|mini-gmp|round_raw_generic).c ]] || \ grep -q '^# *include *"\(mpfr-impl\|fits.*\|gen_inverse\|random_deviate\)\.h"' $file || \ { echo "Missing '#include \"mpfr-impl.h\"' in $file?" && err=1 } done @@ -303,7 +305,7 @@ fdlv2="`sed -n '/GNU Free Documentation License/ { p q }' doc/mpfr.texi`" -[[ $fdlv1 = $fdlv2 ]] || { cat <<EOF && err=1 } +[[ $fdlv1 == $fdlv2 ]] || { cat <<EOF && err=1 } GFDL versions differ: fdl.texi: $fdlv1 mpfr.texi: $fdlv2 @@ -317,10 +319,13 @@ else echo "Warning! gcc is not installed. Cannot run ck-mparam." >&2 fi -texisvnd=`LC_ALL=C TZ=UTC svn info doc/mpfr.texi 2> /dev/null | sed -n 's/Last Changed Date:.*, [0-9]* \([A-Z][a-z][a-z] [0-9][0-9][0-9][0-9]\)).*/\1/p'` +# Note about the TZ value: GMT0 and UTC0 are both specified by POSIX, +# and UTC0 is the preferred value, but old systems only accept GMT0. +# The "0" is important ("GMT" alone does not work on Tru64 Unix). +texisvnd=`LC_ALL=C TZ=GMT0 svn info doc/mpfr.texi 2> /dev/null | sed -n 's/Last Changed Date:.*, [0-9]* \([A-Z][a-z][a-z] [0-9][0-9][0-9][0-9]\)).*/\1/p'` if [[ $? -eq 0 ]] && [[ -n "$texisvnd" ]] then texidate=`sed -n 's/@set UPDATED-MONTH \([A-Z][a-z][a-z]\).*\( [0-9][0-9][0-9][0-9]\)/\1\2/p' doc/mpfr.texi` - [[ $texidate = $texisvnd ]] || { cat <<EOF && err=1 } + [[ $texidate == $texisvnd ]] || { cat <<EOF && err=1 } mpfr.texi's UPDATED-MONTH seems to be incorrect: mpfr.texi's UPDATED-MONTH: $texidate Last Changed Date in WC: $texisvnd @@ -329,7 +334,7 @@ fi acv1=`sed -n 's/.*autoconf \([0-9.]\+\) (at least).*/\1/p' doc/README.dev` acv2=`sed -n 's/AC_PREREQ(\([0-9.]\+\).*/\1/p' acinclude.m4` -[[ $acv1 = $acv2 ]] || { cat <<EOF && err=1 } +[[ $acv1 == $acv2 ]] || { cat <<EOF && err=1 } autoconf minimal versions differ: README.dev: $acv1 acinclude.m4: $acv2 |