summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/ai.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/ai.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/ai.c163
1 files changed, 97 insertions, 66 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/ai.c b/Build/source/libs/mpfr/mpfr-src/src/ai.c
index 58d4addce30..66650dddd16 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/ai.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/ai.c
@@ -1,6 +1,6 @@
/* mpfr_ai -- Airy function Ai
-Copyright 2010-2019 Free Software Foundation, Inc.
+Copyright 2010-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -106,6 +106,9 @@ mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
return mpfr_check_range (y, r, rnd);
}
+ /* now x is not zero */
+ MPFR_ASSERTD(!MPFR_IS_ZERO(x));
+
/* FIXME: underflow for large values of |x| ? */
@@ -133,43 +136,56 @@ mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
mpfr_init2 (tmp2_sp, MPFR_SMALL_PRECISION);
mpfr_abs (tmp_sp, x, MPFR_RNDU);
mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU);
- mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ x^3/2 */
+ mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ |x|^(3/2) */
/* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */
mpfr_set_str (tmp2_sp, "0.96179669392597567", 10, MPFR_RNDU);
mpfr_mul (tmp2_sp, tmp_sp, tmp2_sp, MPFR_RNDU);
/* cond represents the number of lost bits in the evaluation of the sum */
- if ( (MPFR_IS_ZERO (x)) || (MPFR_GET_EXP (x) <= 0) )
+ if (MPFR_GET_EXP (x) <= 0)
cond = 0;
else
- cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU) - (MPFR_GET_EXP (x)-1)/4 - 1;
+ {
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
+ cond -= (MPFR_GET_EXP (x) - 1) / 4 + 1;
+ }
/* The variable assumed_exponent is used to store the maximal assumed */
/* exponent of Ai(x). More precisely, we assume that |Ai(x)| will be */
/* greater than 2^{-assumed_exponent}. */
- if (MPFR_IS_ZERO (x))
- assumed_exponent = 2;
- else
+ if (MPFR_IS_POS (x))
{
- if (MPFR_IS_POS (x))
+ if (MPFR_GET_EXP (x) <= 0)
+ assumed_exponent = 3;
+ else
{
- if (MPFR_GET_EXP (x) <= 0)
- assumed_exponent = 3;
- else
- assumed_exponent = (2 + (MPFR_GET_EXP (x)/4 + 1)
- + mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ unsigned long int t;
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
+ assumed_exponent = t + 2 + (MPFR_GET_EXP (x) / 4 + 1);
+ MPFR_ASSERTN (assumed_exponent > t);
}
- /* We do not know Ai (x) yet */
- /* We cover the case when EXP (Ai (x))>=-10 */
- else
- assumed_exponent = 10;
}
+ /* We do not know Ai (x) yet */
+ /* We cover the case when EXP (Ai (x))>=-10 */
+ else
+ assumed_exponent = 10;
{
- mpfr_prec_t incr =
- MPFR_INT_CEIL_LOG2 (prec) + 5 + cond + assumed_exponent;
- wprec = MPFR_ADD_PREC (prec, incr);
+ unsigned long int t, u;
+
+ t = assumed_exponent + cond;
+ MPFR_ASSERTN (t >= cond);
+ u = MPFR_INT_CEIL_LOG2 (prec) + 5;
+ t += u;
+ MPFR_ASSERTN (t >= u);
+ wprec = MPFR_ADD_PREC (prec, t);
}
mpfr_init (ti);
@@ -244,15 +260,13 @@ mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
MPFR_LOG_MSG (("Roundoff error: %Pu\n", err));
MPFR_LOG_MSG (("Approxim error: %Pu\n", wprec-prec-1));
- if (wprec < err+1)
- correct_bits=0;
+ if (wprec < err + 1)
+ correct_bits = 0;
+ else if (wprec < err + prec +1)
+ correct_bits = wprec - err - 1; /* since wprec > err + 1,
+ correct_bits > 0 */
else
- {
- if (wprec < err+prec+1)
- correct_bits = wprec - err - 1;
- else
- correct_bits = prec;
- }
+ correct_bits = prec;
if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd)))
break;
@@ -262,26 +276,24 @@ mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
assumed_exponent *= 2;
MPFR_LOG_MSG (("Not a single bit correct (assumed_exponent=%lu)\n",
assumed_exponent));
- wprec = prec + 5 + MPFR_INT_CEIL_LOG2 (k) + cond + assumed_exponent;
+ wprec = prec + 5 + MPFR_INT_CEIL_LOG2 (prec) + cond +
+ assumed_exponent;
+ }
+ else if (correct_bits < prec)
+ { /* The precision was badly chosen */
+ MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x)"
+ " (E=%" MPFR_EXP_FSPEC "d)\n",
+ (mpfr_eexp_t) MPFR_GET_EXP (s)));
+ wprec = prec + err + 1;
}
else
- {
- if (correct_bits < prec)
- { /* The precision was badly chosen */
- MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x)"
- " (E=%" MPFR_EXP_FSPEC "d)\n",
- (mpfr_eexp_t) MPFR_GET_EXP (s)));
- wprec = prec + err + 1;
- }
- else
- { /* We are really in a bad case of the TMD */
- MPFR_ZIV_NEXT (loop, prec);
+ { /* We are really in a bad case of the TMD */
+ MPFR_ZIV_NEXT (loop, prec);
- /* We update wprec */
- /* We assume that K will not be multiplied by more than 4 */
- wprec = prec + (MPFR_INT_CEIL_LOG2 (k)+2) + 5 + cond
- - MPFR_GET_EXP (s);
- }
+ /* We update wprec */
+ /* We assume that K will not be multiplied by more than 4 */
+ wprec = prec + (MPFR_INT_CEIL_LOG2 (k) + 2) + 5 + cond
+ - MPFR_GET_EXP (s);
}
} /* End of ZIV loop */
@@ -305,7 +317,7 @@ mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
/* Airy function Ai evaluated by Smith algorithm.
- Assume that x is a finite number. */
+ Assume that x is a finite non-zero number. */
static int
mpfr_ai2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
{
@@ -338,7 +350,6 @@ mpfr_ai2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
/* FIXME: underflow for large values of |x| */
-
/* Set initial precision */
/* See the analysis for the naive evaluation */
@@ -350,40 +361,57 @@ mpfr_ai2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
mpfr_init2 (tmp2_sp, MPFR_SMALL_PRECISION);
mpfr_abs (tmp_sp, x, MPFR_RNDU);
mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU);
- mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ x^3/2 */
+ mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ |x|^(3/2) */
/* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */
mpfr_set_str (tmp2_sp, "0.96179669392597567", 10, MPFR_RNDU);
mpfr_mul (tmp2_sp, tmp_sp, tmp2_sp, MPFR_RNDU);
/* cond represents the number of lost bits in the evaluation of the sum */
- if ( (MPFR_IS_ZERO (x)) || (MPFR_GET_EXP (x) <= 0) )
+ if (MPFR_GET_EXP (x) <= 0)
cond = 0;
else
- cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU) - (MPFR_GET_EXP (x) - 1)/4 - 1;
+ {
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
+ cond -= (MPFR_GET_EXP (x) - 1) / 4 + 1;
+ }
/* This variable is used to store the maximal assumed exponent of */
- /* Ai (x). More precisely, we assume that |Ai (x)| will be greater than */
- /* 2^{-assumedExp}. */
- if (MPFR_IS_ZERO (x))
- assumed_exponent = 2;
- else
+ /* Ai(x). More precisely, we assume that |Ai(x)| will be greater than */
+ /* 2^{-assumed_exponent}. */
+ if (MPFR_IS_POS (x))
{
- if (MPFR_IS_POS (x))
+ if (MPFR_GET_EXP (x) <= 0)
+ assumed_exponent = 3;
+ else
{
- if (MPFR_GET_EXP (x) <= 0)
- assumed_exponent = 3;
- else
- assumed_exponent = (2 + (MPFR_GET_EXP (x)/4 + 1)
- + mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ unsigned long int t;
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
+ MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
+ assumed_exponent = t + 2 + (MPFR_GET_EXP (x) / 4 + 1);
+ MPFR_ASSERTN (assumed_exponent > t);
}
- /* We do not know Ai (x) yet */
- /* We cover the case when EXP (Ai (x))>=-10 */
- else
- assumed_exponent = 10;
}
+ /* We do not know Ai(x) yet */
+ /* We cover the case when EXP(Ai(x))>=-10 */
+ else
+ assumed_exponent = 10;
- wprec = prec + MPFR_INT_CEIL_LOG2 (prec) + 6 + cond + assumed_exponent;
+ {
+ unsigned long int t, u;
+
+ t = assumed_exponent + cond;
+ MPFR_ASSERTN (t >= cond);
+ u = MPFR_INT_CEIL_LOG2 (prec) + 6;
+ t += u;
+ MPFR_ASSERTN (t >= u);
+ wprec = MPFR_ADD_PREC (prec, t);
+ }
/* We assume that the truncation rank will be ~ prec */
L = __gmpfr_isqrt (prec);
@@ -629,6 +657,7 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
}
else if (MPFR_IS_INF (x))
return mpfr_set_ui (y, 0, rnd);
+ /* the cases x = +0 or -0 will be treated below */
}
/* The exponent range must be large enough for the computation of temp1. */
@@ -655,5 +684,7 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
MPFR_SAVE_EXPO_FREE (expo); /* Ignore all previous exceptions. */
+ /* we use ai2 if |x|*AI_THRESHOLD1/3 + PREC(y)*AI_THRESHOLD2 > AI_SCALE,
+ which means x cannot be zero in mpfr_ai2 */
return use_ai2 ? mpfr_ai2 (y, x, rnd) : mpfr_ai1 (y, x, rnd);
}