summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/pow_z.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/pow_z.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/pow_z.c373
1 files changed, 0 insertions, 373 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/pow_z.c b/Build/source/libs/mpfr/mpfr-src/src/pow_z.c
deleted file mode 100644
index 3b4e62fd684..00000000000
--- a/Build/source/libs/mpfr/mpfr-src/src/pow_z.c
+++ /dev/null
@@ -1,373 +0,0 @@
-/* mpfr_pow_z -- power function x^z with z a MPZ
-
-Copyright 2005-2020 Free Software Foundation, Inc.
-Contributed by the AriC and Caramba projects, INRIA.
-
-This file is part of the GNU MPFR Library.
-
-The GNU MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-The GNU MPFR Library is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
-License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
-https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
-51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
-
-#define MPFR_NEED_LONGLONG_H
-#include "mpfr-impl.h"
-
-/* y <- x^|z| with z != 0
- if cr=1: ensures correct rounding of y
- if cr=0: does not ensure correct rounding, but avoid spurious overflow
- or underflow, and uses the precision of y as working precision (warning,
- y and x might be the same variable). */
-static int
-mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t rnd, int cr)
-{
- mpfr_t res;
- mpfr_prec_t prec, err;
- int inexact;
- mpfr_rnd_t rnd1, rnd2;
- mpz_t absz;
- mp_size_t size_z;
- MPFR_ZIV_DECL (loop);
- MPFR_BLOCK_DECL (flags);
-
- MPFR_LOG_FUNC
- (("x[%Pu]=%.*Rg z=%Zd rnd=%d cr=%d",
- mpfr_get_prec (x), mpfr_log_prec, x, z, rnd, cr),
- ("y[%Pu]=%.*Rg inexact=%d",
- mpfr_get_prec (y), mpfr_log_prec, y, inexact));
-
- MPFR_ASSERTD (mpz_sgn (z) != 0);
-
- if (MPFR_UNLIKELY (mpz_cmpabs_ui (z, 1) == 0))
- return mpfr_set (y, x, rnd);
-
- absz[0] = z[0];
- SIZ (absz) = ABSIZ (absz); /* Hack to get abs(z) */
- MPFR_MPZ_SIZEINBASE2 (size_z, z);
-
- /* round toward 1 (or -1) to avoid spurious overflow or underflow,
- i.e. if an overflow or underflow occurs, it is a real exception
- and is not just due to the rounding error. */
- rnd1 = (MPFR_EXP(x) >= 1) ? MPFR_RNDZ
- : (MPFR_IS_POS(x) ? MPFR_RNDU : MPFR_RNDD);
- rnd2 = (MPFR_EXP(x) >= 1) ? MPFR_RNDD : MPFR_RNDU;
-
- if (cr != 0)
- prec = MPFR_PREC (y) + 3 + size_z + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
- else
- prec = MPFR_PREC (y);
- mpfr_init2 (res, prec);
-
- MPFR_ZIV_INIT (loop, prec);
- for (;;)
- {
- unsigned int inexmul; /* will be non-zero if res may be inexact */
- mp_size_t i = size_z;
-
- /* now 2^(i-1) <= z < 2^i */
- /* see below (case z < 0) for the error analysis, which is identical,
- except if z=n, the maximal relative error is here 2(n-1)2^(-prec)
- instead of 2(2n-1)2^(-prec) for z<0. */
- MPFR_ASSERTD (prec > (mpfr_prec_t) i);
- err = prec - 1 - (mpfr_prec_t) i;
-
- MPFR_BLOCK (flags,
- inexmul = mpfr_sqr (res, x, rnd2);
- MPFR_ASSERTD (i >= 2);
- if (mpz_tstbit (absz, i - 2))
- inexmul |= mpfr_mul (res, res, x, rnd1);
- for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
- {
- inexmul |= mpfr_sqr (res, res, rnd2);
- if (mpz_tstbit (absz, i))
- inexmul |= mpfr_mul (res, res, x, rnd1);
- });
- if (MPFR_LIKELY (inexmul == 0 || cr == 0
- || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
- || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
- break;
- /* Can't decide correct rounding, increase the precision */
- MPFR_ZIV_NEXT (loop, prec);
- mpfr_set_prec (res, prec);
- }
- MPFR_ZIV_FREE (loop);
-
- /* Check Overflow */
- if (MPFR_OVERFLOW (flags))
- {
- MPFR_LOG_MSG (("overflow\n", 0));
- inexact = mpfr_overflow (y, rnd, mpz_odd_p (absz) ?
- MPFR_SIGN (x) : MPFR_SIGN_POS);
- }
- /* Check Underflow */
- else if (MPFR_UNDERFLOW (flags))
- {
- MPFR_LOG_MSG (("underflow\n", 0));
- if (rnd == MPFR_RNDN)
- {
- mpfr_t y2, zz;
-
- /* We cannot decide now whether the result should be rounded
- toward zero or +Inf. So, let's use the general case of
- mpfr_pow, which can do that. But the problem is that the
- result can be exact! However, it is sufficient to try to
- round on 2 bits (the precision does not matter in case of
- underflow, since MPFR does not have subnormals), in which
- case, the result cannot be exact due to previous filtering
- of trivial cases. */
- MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
- MPFR_EXP (x) - 1) != 0);
- mpfr_init2 (y2, 2);
- mpfr_init2 (zz, ABSIZ (z) * GMP_NUMB_BITS);
- inexact = mpfr_set_z (zz, z, MPFR_RNDN);
- MPFR_ASSERTN (inexact == 0);
- inexact = mpfr_pow_general (y2, x, zz, rnd, 1,
- (mpfr_save_expo_t *) NULL);
- mpfr_clear (zz);
- mpfr_set (y, y2, MPFR_RNDN);
- mpfr_clear (y2);
- __gmpfr_flags = MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
- }
- else
- {
- inexact = mpfr_underflow (y, rnd, mpz_odd_p (absz) ?
- MPFR_SIGN (x) : MPFR_SIGN_POS);
- }
- }
- else
- inexact = mpfr_set (y, res, rnd);
-
- mpfr_clear (res);
- return inexact;
-}
-
-/* The computation of y = pow(x,z) is done by
- * y = set_ui(1) if z = 0
- * y = pow_ui(x,z) if z > 0
- * y = pow_ui(1/x,-z) if z < 0
- *
- * Note: in case z < 0, we could also compute 1/pow_ui(x,-z). However, in
- * case MAX < 1/MIN, where MAX is the largest positive value, i.e.,
- * MAX = nextbelow(+Inf), and MIN is the smallest positive value, i.e.,
- * MIN = nextabove(+0), then x^(-z) might produce an overflow, whereas
- * x^z is representable.
- */
-
-int
-mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t rnd)
-{
- int inexact;
- mpz_t tmp;
- MPFR_SAVE_EXPO_DECL (expo);
-
- MPFR_LOG_FUNC
- (("x[%Pu]=%.*Rg z=%Zd rnd=%d",
- mpfr_get_prec (x), mpfr_log_prec, x, z, rnd),
- ("y[%Pu]=%.*Rg inexact=%d",
- mpfr_get_prec (y), mpfr_log_prec, y, inexact));
-
- /* x^0 = 1 for any x, even a NaN */
- if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
- return mpfr_set_ui (y, 1, rnd);
-
- if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
- {
- if (MPFR_IS_NAN (x))
- {
- MPFR_SET_NAN (y);
- MPFR_RET_NAN;
- }
- else if (MPFR_IS_INF (x))
- {
- /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
- /* Inf ^(-n) = 0, sign = + if x>0 or z even */
- if (mpz_sgn (z) > 0)
- MPFR_SET_INF (y);
- else
- MPFR_SET_ZERO (y);
- if (MPFR_UNLIKELY (MPFR_IS_NEG (x) && mpz_odd_p (z)))
- MPFR_SET_NEG (y);
- else
- MPFR_SET_POS (y);
- MPFR_RET (0);
- }
- else /* x is zero */
- {
- MPFR_ASSERTD (MPFR_IS_ZERO(x));
- if (mpz_sgn (z) > 0)
- /* 0^n = +/-0 for any n */
- MPFR_SET_ZERO (y);
- else
- {
- /* 0^(-n) if +/- INF */
- MPFR_SET_INF (y);
- MPFR_SET_DIVBY0 ();
- }
- if (MPFR_LIKELY (MPFR_IS_POS (x) || mpz_even_p (z)))
- MPFR_SET_POS (y);
- else
- MPFR_SET_NEG (y);
- MPFR_RET(0);
- }
- }
-
- MPFR_SAVE_EXPO_MARK (expo);
-
- /* detect exact powers: x^-n is exact iff x is a power of 2
- Do it if n > 0 too as this is faster and this filtering is
- needed in case of underflow. */
- if (MPFR_UNLIKELY (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
- MPFR_EXP (x) - 1) == 0))
- {
- mpfr_exp_t expx = MPFR_EXP (x); /* warning: x and y may be the same
- variable */
-
- MPFR_LOG_MSG (("x^n with x power of two\n", 0));
- mpfr_set_si (y, mpz_odd_p (z) ? MPFR_INT_SIGN(x) : 1, rnd);
- MPFR_ASSERTD (MPFR_IS_FP (y));
- mpz_init (tmp);
- mpz_mul_si (tmp, z, expx - 1);
- MPFR_ASSERTD (MPFR_GET_EXP (y) == 1);
- mpz_add_ui (tmp, tmp, 1);
- inexact = 0;
- if (MPFR_UNLIKELY (mpz_cmp_si (tmp, __gmpfr_emin) < 0))
- {
- MPFR_LOG_MSG (("underflow\n", 0));
- /* |y| is a power of two, thus |y| <= 2^(emin-2), and in
- rounding to nearest, the value must be rounded to 0. */
- if (rnd == MPFR_RNDN)
- rnd = MPFR_RNDZ;
- inexact = mpfr_underflow (y, rnd, MPFR_SIGN (y));
- }
- else if (MPFR_UNLIKELY (mpz_cmp_si (tmp, __gmpfr_emax) > 0))
- {
- MPFR_LOG_MSG (("overflow\n", 0));
- inexact = mpfr_overflow (y, rnd, MPFR_SIGN (y));
- }
- else
- MPFR_SET_EXP (y, mpz_get_si (tmp));
- mpz_clear (tmp);
- MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
- }
- else if (mpz_sgn (z) > 0)
- {
- inexact = mpfr_pow_pos_z (y, x, z, rnd, 1);
- MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
- }
- else
- {
- /* Declaration of the intermediary variable */
- mpfr_t t;
- mpfr_prec_t Nt; /* Precision of the intermediary variable */
- mpfr_rnd_t rnd1;
- mp_size_t size_z;
- MPFR_ZIV_DECL (loop);
-
- MPFR_MPZ_SIZEINBASE2 (size_z, z);
-
- /* initial working precision */
- Nt = MPFR_PREC (y);
- Nt = Nt + size_z + 3 + MPFR_INT_CEIL_LOG2 (Nt);
- /* ensures Nt >= bits(z)+2 */
-
- /* initialize of intermediary variable */
- mpfr_init2 (t, Nt);
-
- /* We will compute rnd(rnd1(1/x) ^ (-z)), where rnd1 is the rounding
- toward sign(x), to avoid spurious overflow or underflow. */
- rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ :
- (MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD);
-
- MPFR_ZIV_INIT (loop, Nt);
- for (;;)
- {
- MPFR_BLOCK_DECL (flags);
-
- /* compute (1/x)^(-z), -z>0 */
- /* As emin = -emax, an underflow cannot occur in the division.
- And if an overflow occurs, then this means that x^z overflows
- too (since we have rounded toward 1 or -1). */
- MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
- MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
- /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
- if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
- goto overflow;
- MPFR_BLOCK (flags, mpfr_pow_pos_z (t, t, z, rnd, 0));
- /* Now if z=-n, t = x^z*(1+theta)^(2n-1) where |theta| <= 2^(-Nt),
- with theta maybe different from above. If (2n-1)*2^(-Nt) <= 1/2,
- which is satisfied as soon as Nt >= bits(z)+2, then we can use
- Lemma \ref{lemma_graillat} from algorithms.tex, which yields
- t = x^z*(1+theta) with |theta| <= 2(2n-1)*2^(-Nt), thus the
- error is bounded by 2(2n-1) ulps <= 2^(bits(z)+2) ulps. */
- if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
- {
- overflow:
- MPFR_ZIV_FREE (loop);
- mpfr_clear (t);
- MPFR_SAVE_EXPO_FREE (expo);
- MPFR_LOG_MSG (("overflow\n", 0));
- return mpfr_overflow (y, rnd,
- mpz_odd_p (z) ? MPFR_SIGN (x) :
- MPFR_SIGN_POS);
- }
- if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
- {
- MPFR_ZIV_FREE (loop);
- mpfr_clear (t);
- MPFR_LOG_MSG (("underflow\n", 0));
- if (rnd == MPFR_RNDN)
- {
- mpfr_t y2, zz;
-
- /* We cannot decide now whether the result should be
- rounded toward zero or away from zero. So, like
- in mpfr_pow_pos_z, let's use the general case of
- mpfr_pow in precision 2. */
- MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
- MPFR_EXP (x) - 1) != 0);
- mpfr_init2 (y2, 2);
- mpfr_init2 (zz, ABSIZ (z) * GMP_NUMB_BITS);
- inexact = mpfr_set_z (zz, z, MPFR_RNDN);
- MPFR_ASSERTN (inexact == 0);
- inexact = mpfr_pow_general (y2, x, zz, rnd, 1,
- (mpfr_save_expo_t *) NULL);
- mpfr_clear (zz);
- mpfr_set (y, y2, MPFR_RNDN);
- mpfr_clear (y2);
- MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
- goto end;
- }
- else
- {
- MPFR_SAVE_EXPO_FREE (expo);
- return mpfr_underflow (y, rnd, mpz_odd_p (z) ?
- MPFR_SIGN (x) : MPFR_SIGN_POS);
- }
- }
- if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_z - 2, MPFR_PREC (y),
- rnd)))
- break;
- /* actualization of the precision */
- MPFR_ZIV_NEXT (loop, Nt);
- mpfr_set_prec (t, Nt);
- }
- MPFR_ZIV_FREE (loop);
-
- inexact = mpfr_set (y, t, rnd);
- mpfr_clear (t);
- }
-
- end:
- MPFR_SAVE_EXPO_FREE (expo);
- return mpfr_check_range (y, inexact, rnd);
-}